Hi,
Your procedure is totally correct to do what you want to do.
About the B*, the SGP4 specific ballistic coefficient, the ratio from 1 to 3 (input vs output) is not so bad. Theoretically, as related to drag, it depends on Cd, A, m and also from the atmospheric model in your case, when converting from TLEPropagator to NumericalPropagator and then back to TLE. When building the forces, you set the Cd and the drag area and the HarrisPriester as the atmospheric model, the mass is set by default to 1000 kg by Orekit. All these parameters affect the double fitting process and changing any of them will conduct to a different B*. Recovering the drag often conducts to indetermination.
Setting for example the mass to 3000 kg:
propagatorBuilder.setMass(3000.);
gives to the following results:
INPUT TLE
1 31135U 07013A 11003.00000000 .00000816 00000+0 47577-4 0 11
2 31135 2.4656 183.9084 0021119 236.4164 60.4567 15.10546832 15
OUTPUT TLE
1 31135U 07013A 11003.00000000 .00000000 00000-0 48942-4 0 14
2 31135 2.4656 183.9084 0021119 236.4164 60.4567 15.10546809 08
The output B* is now closer from the input and the orbital elements are also better.
At last, you set a duration of one day, 86400 s, for fitting that is more than 15 orbital periods and is a bit long in this case, try with 2 orbital periods as a duration, it will greatly reduce the processing time with comparable results.