Fast Calculation Method for Predicting the Morphology of Steady-State Ablation

: Predicting the surface morphology of materials during steady-state ablation is important in rocket motor nozzles and the heat shields of vehicles performing atmospheric re-entry. When designing ablative materials, a high number of calculations is required for analyzing surface morphology. To effectively design these materials and reduce the number of experiments, a fast, effective, and simple calculation method is required. Although a fundamental theory for ablation has been established, quick and effective prediction of the morphology of the composites remains a challenge. In this study, we propose a fast, effective, and simple numerical calculation method to predict the surface morphology of steady-state ablation based on the geometric characteristics of the materials. The results obtained in this study were consistent with the experimental observations. The calculation time was signiﬁcantly reduced. In addition, our method was found to be useful for analyzing the physical and chemical properties and surface roughness of ablative materials.


Introduction
The thermal protection system (TPS) is considered essential technology in the development of hypersonic vehicles [1][2][3][4][5][6][7][8].The ideal materials for thermal protection systems can be divided into three categories according to different ablation mechanisms: molten, pyrolysis, and gasification types.Among these, the excellent high-temperature properties of gasification-type ablative materials, including high-temperature stability and good mechanical properties, have been extensively studied [9][10][11][12].The received heat flux of the materials is converted into an outward mass flux through physical or chemical processes, which cause surface recession and destruction of the materials.Surface recession leads to a change in the structural shape, which is an inevitable phenomenon with various hazardous repercussions [13,14].The change in the structural shape affects the flight trajectory and attitude-control characteristics of the aircraft.Meanwhile, the surface recession leads to regular morphology of the material surface, making the material surface progressively rougher.A rough surface not only induces laminar-to-turbulent transition in the boundary layer, but also increases the aerodynamic heating rate of the interaction between the materials and airflow [15].Therefore, ablation morphology must be considered when designing ablative materials.Owing to the large number of calculations required for analyzing the material surface morphology, a fast, effective, and simple calculation method is required.
Analyzing the surface morphology of materials during steady-state ablation is important for the efficient design of thermal structure protection.Previously, thermochemical ablation behavior and data have been obtained by high-temperature flame experiments and observed by scanning electron microscopy (SEM) [16][17][18][19][20][21][22][23][24].The results showed that the fibers and matrix were exposed to the atmospheric environment, and there was a thin reactive interphase between the fiber and matrix [25][26][27][28][29][30].Fan et al. [31,32] discovered the needle shapes and sheath package structures of composites, following ablation using a hypersonic flowing propane torch.The reasons for the microstructural evolution of the C/C and C/C-SiC composites were discussed and analyzed using ANSYS Fluent software.The results showed that the temperature of the fibers did not change significantly, and that the distribution of oxygen concentration at the edges of the fibers was possibly the main reason for the formation of needle shapes.Wang et al. [33] examined laser ablation behaviors of C/SiC composites in a hypersonic wind tunnel.The results showed significant differences in ablation morphology between micro and macro scales.Under high supersonic flow, the needle-shaped microstructure of carbon fiber became rough and tattered under static air.In addition, Nikan et al. [34][35][36][37][38] proposed numerical approaches for modeling heat conduction and gas diffusion.Numerical results confirmed the accuracy and efficiency of the numerical approaches.
Based on analysis of the regular ablation morphology of composites [39], Lachaud et al. [40] assumed that regular ablation morphology reached a stable state.They then established a theoretical analysical model to predict the steady-state ablation morphology of composite components.The results showed that the dimensionless relative reactivity and Sherwood number were two of the parameters necessary to describe the steady-state ablation morphology characteristics of the bundle.They discussed the morphological characteristics of the reaction-limited, transition, and diffusion-limited regimes.Lachaud et al. [41] also established an analytical model for transient ablation morphology.However, the research samples in these studies were simple circular fibers and rectangular fiber bundles.Li et al. [42] studied the microstructural evolution of carbon fibers during ablation in air, as indicated in Figure 1a.The fiber cross-section was not a perfect circle but was slightly undulated around a mean circular line.To simulate more accurately the surface ablation morphology of realistic fibers, a Brownian motion simulation method was implemented [43].The morphological evolution of the model from the initial state to the steady state was complex, and included molecular Brownian motion and 3D volume element deletion.In addition, the cross section of the fiber bundles showed elliptic and rectangular shapes, as in Figure 1b,c [44,45].It is also necessary to study the steady-state ablation morphology of different cross sections.
Coatings 2022, 12, x FOR PEER REVIEW 3 of 16 each point of the complex cross section are defined in Sections 2.3 and 2.4, respectively.A fast calculation method using a theoretical steady-state ablation morphology model is proposed, which can directly and quickly predict the steady-state ablation morphology of composites with arbitrary sections.The calculated results were compared with the experimental and simulation results.The errors are stated, confirming the accuracy of the model.The results provide a basis for analyzing the physical/ and chemical properties and surface roughness of ablative materials.The surface roughness features may provide a basis for better understanding the effect of surface morphology on aerodynamic analysis.

Problem Statement
Gasification-type ablative materials usually sublimate or react with air at high temperatures.In these reactions, the products are all gases, causing mass loss and surface recession.Therefore, the residue of reaction products is not considered when assessing the ablation morphology of the material.It is difficult for gas products to protect gasification-type ablative materials, and they cannot effectively weaken the ablation behavior.However, divergent ablation performance in thermal protection materials is not appropriate.Thus, a qualified gasification-type ablative material can achieve steady-state abla- Due to the complexity of previous models, the speed and effectiveness of morphology prediction for composites have become concerning issues.For this purpose, a series of numerical simulation methods have been used to study the ablation morphology evolution of composites [46], such as the level-set method [47,48], volume-of-fluid (VOF) method [49][50][51], and lattice Boltzmann method [52].It is worth noting that the compu-tational process for the steady-state ablation morphology of the realistic fiber simulated using the above methods remains complex and time-consuming.There are two reasons for this; (1) the steady-state ablation morphology needs to be calculated from the initial state to the transient state to the steady state, and (2) the current numerical simulation methods include gas diffusion, chemical reaction (recession criterion), and deformation or deletion of 3D volume elements to realize boundary movement.However, the steady-state ablation of composites meant that their ablation morphology did not change.The ablation morphology was independent of the ablation process under specified ablation conditions.The computation did not need to track the entire evolution process from transient to steadystate ablation.Based on these steady-state characteristics, an effective and simple method can be used to rapidly calculate morphology.
As reviewed above, no method that can quickly predict the steady-state ablation morphology of composites.Therefore, this study investigates the characteristics of steady-state ablation morphology from the perspective of mathematics and geometry.The characteristic length and direction of the steady-state ablation morphology corresponding to each point of the complex cross section are defined in Sections 2.3 and 2.4, respectively.A fast calculation method using a theoretical steady-state ablation morphology model is proposed, which can directly and quickly predict the steady-state ablation morphology of composites with arbitrary sections.The calculated results were compared with the experimental and simulation results.The errors are stated, confirming the accuracy of the model.The results provide a basis for analyzing the physical/ and chemical properties and surface roughness of ablative materials.The surface roughness features may provide a basis for better understanding the effect of surface morphology on aerodynamic analysis.

Problem Statement
Gasification-type ablative materials usually sublimate or react with air at high temperatures.In these reactions, the products are all gases, causing mass loss and surface recession.Therefore, the residue of reaction products is not considered when assessing the ablation morphology of the material.It is difficult for gas products to protect gasification-type ablative materials, and they cannot effectively weaken the ablation behavior.However, divergent ablation performance in thermal protection materials is not appropriate.Thus, a qualified gasification-type ablative material can achieve steady-state ablation in a certain environment.
Carbon-fiber-reinforced-carbon-matrix composites (C/C composites), a gasificationtype ablative material, are taken as an example in the following description.Within a certain temperature range, the gas in the boundary layer adheres to the material surface and reacts with the material.C/C composites easily react with oxygen at lower temperatures (approximately 370 • C), and the carbon-nitrogen reactions and sublimation occur at higher temperatures.The gas products escape from the material surface into the boundary layer.The thin interphases lying between the fibers and matrix are highly reactive, and the ablation recession velocity of the interphases tends to be greater than that of the fibers and matrix.Therefore, the interphases were oxidized first, and both the fibers and the matrix were exposed to a gaseous environment, as shown in Figure 2a,b.The surface of the material changed from smooth to rough.Finally, the surface recession velocity and morphology of the C/C composites were close to a steady state with the continuous progress of various reactions, as in Figure 2c,d.

Ablation Morphology Theoretical Model
Considering the position of the fiber, i.e., vertical, inclined, or rotated around the axis, the steady-state ablation morphology was determined, which only changed with the coordinates.In steady-state ablation, the calculation can be significantly simplified using a morphology-control equation.To simplify the description of the model, we assumed that the fiber was an ideal pillar, vertically consistent, and added only the case in which the fiber was perpendicular to the material surface.This implies that the steady-state ablation morphology of the material only receded in the negative direction along the Z axis.The steady-state ablation morphology was not rotated around the Z axis.Because the fiber boundary PAPB (shown in Figure 3a) in the initial state was parallel to the Z axis, the morphology evolution PA2PB2 (shown in Figure 3b) of PAPB was also in this plane parallel to the Z axis.The evolution of the steady-state ablation morphology from PA1PB1 of t1 to PA2PB2 of t1 + Δt is shown in detail in Figure 3c.It is evident that the local real ablation velocity ve (m⋅s -1 ) of point PE1 was perpendicular to the fiber surface.The global ablation velocity va (m⋅s -1 ) of point PE1 was parallel to the Z axis and perpendicular to the material macrosurface, and equal to the ablation velocity vi (m⋅s -1 ) of the interphase.In addition, it is acceptable to assume that the fiber surface was flat within a sufficiently small region.Thus, the geometry formed by connecting ve and va can be represented by a triangle.The angle between ve and va is θ (°) and θ < 90°.According to the geometric relationship of the two velocities, the steady-state ablation condition can be written as: where x, y, and z are the lengths in the X, Y, and Z axis direction (m), respectively.When the steady-state ablation morphology is in the LZ plane, (Δl) 2 = (Δx) 2 + (Δy) 2 , Equation (2) can be replaced by:

Ablation Morphology Theoretical Model
Considering the position of the fiber, i.e., vertical, inclined, or rotated around the axis, the steady-state ablation morphology was determined, which only changed with the coordinates.In steady-state ablation, the calculation can be significantly simplified using a morphology-control equation.To simplify the description of the model, we assumed that the fiber was an ideal pillar, vertically consistent, and added only the case in which the fiber was perpendicular to the material surface.This implies that the steady-state ablation morphology of the material only receded in the negative direction along the Z axis.The steady-state ablation morphology was not rotated around the Z axis.Because the fiber boundary P A P B (shown in Figure 3a) in the initial state was parallel to the Z axis, the morphology evolution P A2 P B2 (shown in Figure 3b) of P A P B was also in this plane parallel to the Z axis.The evolution of the steady-state ablation morphology from P A1 P B1 of t 1 to P A2 P B2 of t 1 + ∆t is shown in detail in Figure 3c.It is evident that the local real ablation velocity v e (m•s −1 ) of point P E1 was perpendicular to the fiber surface.The global ablation velocity v a (m•s −1 ) of point P E1 was parallel to the Z axis and perpendicular to the material macro-surface, and equal to the ablation velocity v i (m•s −1 ) of the interphase.In addition, it is acceptable to assume that the fiber surface was flat within a sufficiently small region.Thus, the geometry formed by connecting v e and v a can be represented by a triangle.The angle between v e and v a is θ ( • ) and θ < 90 • .According to the geometric relationship of the two velocities, the steady-state ablation condition can be written as: where x, y, and z are the lengths in the X, Y, and Z axis direction (m), respectively.When the steady-state ablation morphology is in the LZ plane, (∆l) 2 = (∆x) 2 + (∆y) 2 , Equation ( 2) can be replaced by: where cosθ > 0 and ∆l/∆z > 0. The L axis is in the same direction as the inner normal direction n of point P A2 .
where cosθ > 0 and Δl Δz ⁄ > 0. The L axis is in the same direction as the inner normal direction n of point PA2.In addition, the ablation velocity v is denoted [53] by: where Ω is the solid molar volume of the phase (m 3 ⋅mol -1 ), J is the molar rate of ablation (mol⋅m -2 ⋅s -1 ), and nv is the velocity direction normal to the surface and pointing towards the interior.
Assuming that the chemical reactions of C/C composites in the ablation process satisfy the first-order kinetic law, and the temperature is uniform at the microscopic scale [40], then the molar rate of ablation J is related only to the material reactivity and gas concentration.Thus, the molar rate of ablation J can be calculated as follows: where k is the reactivity of the phase (m⋅s -1 ) and C is the molar concentration (mol⋅m -3 ), which is related to diffusion of molecules.According to the mass conservation equation for gas diffusion, ∂C/∂t + ∇⋅(-D∇C) = 0, and the boundary conditions of Ci = C(z = 0) and C0 = C(z = hf + δc), the molar concentration C(z) can be written as: where Ci and C0 are the molar concentrations of the interphase and boundary layer (mol⋅m -3 ), respectively.Ci = C0/(1 + Dai), where Dai = ki(hf + δc)/D is the Damköhler number, hf is the characteristic height of the fiber (m), δc is the thickness of the concentration boundary layer (m), and D is the diffusion coefficient (m 2 ⋅s -1 ).
Incorporating Equations ( 3)-( 6) into Equation ( 1), Δl Δz ⁄ can be expressed [40] as: In addition, the ablation velocity v is denoted [53] by: where Ω is the solid molar volume of the phase (m 3 •mol −1 ), J is the molar rate of ablation (mol•m −2 •s −1 ), and n v is the velocity direction normal to the surface and pointing towards the interior.
Assuming that the chemical reactions of C/C composites in the ablation process satisfy the first-order kinetic law, and the temperature is uniform at the microscopic scale [40], then the molar rate of ablation J is related only to the material reactivity and gas concentration.Thus, the molar rate of ablation J can be calculated as follows: where k is the reactivity of the phase (m•s −1 ) and C is the molar concentration (mol•m −3 ), which is related to diffusion of molecules.According to the mass conservation equation for gas diffusion, ∂C/∂t + ∇•(−D∇C) = 0, and the boundary conditions of , the molar concentration C(z) can be written as: where C i and C 0 are the molar concentrations of the interphase and boundary layer (mol•m −3 ), respectively.C i = C 0 /(1 + Da i ), where Da i = k i (h f + δ c )/D is the Damköhler number, h f is the characteristic height of the fiber (m), δ c is the thickness of the concentration boundary layer (m), and D is the diffusion coefficient (m 2 •s −1 ).Incorporating Equations ( 3)-( 6) into Equation ( 1), ∆l/∆z can be expressed [40] as: where η = (k i Ω i )/(k f Ω f ) is a dimensionless number that represents the relative reactivity of the interphase and fiber, and χ = D/k i is the relative parameter of diffusion and reactivity (m).
The steady-state ablation morphology, S(l, z) = 0, is determined by l and z.The relationship between length l(z) and height z at each point of the fibers is given by the integration of Equation (7) [40]: where l c is the characteristic length of the point on the fiber cross-section.Clearly, the variables which affect the relationship between length l(z) and height z are η, χ and l c .η and χ are related to the material properties and environmental parameters, respectively, and l c is related to the cross-section of the carbon fibers.
The contour arithmetic mean deviation Sa of the cross-sectional area S is one of the parameters that can be used to describe the roughness, which can be defined as:

Definition of the Characteristic Length l c
When studying steady-state ablation morphology, it is very important to define the characteristic length l c of each point on the fiber cross-section.For an ideal cylindrical fiber, the characteristic length l c of each point on the fiber's cross section can be defined as the fiber radius.However, the cross section is far from a perfect circle, but has an irregular shape.It can be difficult to determine the fiber radius to describe the characteristic length l c .
According to Equation (8), we can draw a steady-state ablation morphology, as shown in Figure 4a.If the characteristic length l c of each point on the fiber's cross section is defined as the curvature radius, it is easy to encounter a situation in which the curvature radius is infinite.Without the influence of other factors, the morphology could extend unboundedly and form a platform at the top.However, the fiber cross section is not infinite.The black dotted line in Figure 4a shows the morphology truncated by another surface; the vertex of the morphology corresponds to P C on the fiber cross-section.Line P A P C is the characteristic length l c of point P A .Therefore, an infinite radius of curvature is unreasonable.
To determine the value of the line P A P C (in Figure 4a), a triangular fiber cross section in the XY plane was considered.As shown in Figure 4b, point P A is the midpoint of the side of the triangular fiber.Three circles pass through point P A .The centers of these circles are points P B , P C , and P D .The inner normal direction n of point P A on the fiber cross-section is in the same direction as the normal of point P A to the circle centers.Because of the particularity of point P A , here we only interpret the effect on the morphology of point P A resulting from the points on the left side of the triangle.Lines P B P E , P C P F , and P D P G are perpendicular to the side of the triangular fiber.This means that points P E , P F , and P G may affect the morphology of point P A at points P B , P C , and P D , respectively.However, when r 1 > r B , the morphology height z BE (according to Equation ( 8)) corresponding to points P B and P E is larger than z BA at point P B .Because r 3 < r D , height z DG is smaller than z DA at point P D .Clearly, neither r B nor r D is a suitable characteristic length l c for point P A .The red circle is the largest circle passing through point P A in the domain of the triangular fiber cross-section.As shown in Figure 4b, P C is a special center, where r 2 = r C and height z CF = z CA at point P C .This indicates that the steady-state ablation morphology of point P A was truncated by the steady-state ablation morphology of point P F above point P C .Therefore, r C is the optimum characteristic length l c for point P A .The largest circle (red circle) tangent to point P A in the domain of the triangular fiber cross-section is the characteristic circle of point P A .
There are two special points, P H and P I , in Figure 4c.Points P H and P I are the convex and concave singularities, respectively.It is easy to determine that the characteristic length l c for point P H is zero.However, the characteristic length and inner normal direction of point P I are not unique; the steady-state ablation morphology of point P I has been truncated by the morphology of the points on lines P J P K and P L .The green lines and curves in Figure 4c indicate the morphology vertices and characteristic centers of each point on the fiber cross-section.As shown in Figure 4d, the definition of the characteristic circle and length is also suitable for a fiber with a concave complex cross section.
the optimum characteristic length lc for point PA.The largest circle (red circle) tangent to point PA in the domain of the triangular fiber cross-section is the characteristic circle of point PA.
There are two special points, PH and PI, in Figure 4c.Points PH and PI are the convex and concave singularities, respectively.It is easy to determine that the characteristic length lc for point PH is zero.However, the characteristic length and inner normal direction of point PI are not unique; the steady-state ablation morphology of point PI has been truncated by the morphology of the points on lines PJPK and PL.The green lines and curves in Figure 4c indicate the morphology vertices and characteristic centers of each point on the fiber cross-section.As shown in Figure 4d, the definition of the characteristic circle and length is also suitable for a fiber with a concave complex cross section.

Calculation Process
As mentioned in Equation ( 8), the steady-state ablation morphology of the fiber can be determined by three parameters: the relative reactivity η, relative parameter of diffusion and reactivity χ, and characteristic length lc.In other words, the steady-state ablation morphology of an infinite extension can be directly determined according to the environment and material properties.The accurate steady-state ablation morphology only needs to be truncated according to the characteristic length lc of each point and placed along the normal direction n on the coordinates of each point.Thus, the calculation process for the steady-state ablation morphology is simple, as shown in Figure 5.
The first step of the experimental method was to obtain data on the cross-sectional profile shape.The fibers and matrix were distinguished using different colors, and we obtained and selected the data from the microscopic picture of the composite material using MATLAB.Additional data points were supplemented by fitting the profile shapes.Thus, the predicted steady-state ablation morphology was more detailed.It is also possible to write the profile shape functions directly into the MATLAB program.

Calculation Process
As mentioned in Equation ( 8), the steady-state ablation morphology of the fiber can be determined by three parameters: the relative reactivity η, relative parameter of diffusion and reactivity χ, and characteristic length l c .In other words, the steady-state ablation morphology of an infinite extension can be directly determined according to the environment and material properties.The accurate steady-state ablation morphology only needs to be truncated according to the characteristic length l c of each point and placed along the normal direction n on the coordinates of each point.Thus, the calculation process for the steady-state ablation morphology is simple, as shown in Figure 5.The second step was to determine the inner normal direction of each point on the cross-sectional profile shape.Owing to the sufficiently dense data points, the inner normal direction n of point Pi can be computed using the surrounding data points Pi-1 and Pi+1.As illustrated in Figure 6, the signs of nx and ny must be determined based on the sizes of xi-1 and xi+1, as well as the sizes of yi-1 and yi+1.As the data are stored clockwise, when xi-1 < xi+1 and yi-1 < yi+1, the sign of n is (+, -).When xi-1 < xi+1 and yi-1 > yi+1, the sign of n is (-, -).When xi-1 > xi+1 and yi-1 > yi+1, the sign of n is (-, +).When xi-1 > xi+1 and yi-1 < yi+1, the sign of n is (+, +).When xi-1 < xi+1 and yi-1 = yi+1, the value of n is (0, -1).When xi-1 = xi+1 and yi-1 > yi+1, the value of n is (-1, 0).When xi-1 > xi+1 and yi-1 = yi+1, the value of n is (0, 1).When xi-1 = xi+1 and yi-1 < yi+1, the value of n is (1, 0).The unit normal direction (nx, ny) of point Pi can be expressed as: , when x i-1 ≠ x i+1 and y i-1 ≠ y i+1 , If the profile shape is written directly into the MATLAB program with functions, concave singularities may exist and need to be determined.The two extreme inner normal directions of the concave singularity can be calculated according to lines Pi-1Pi and PiPi+1.The remaining directions are selected every 1° between the two extreme inner normal directions.The first step of the experimental method was to obtain data on the cross-sectional profile shape.The fibers and matrix were distinguished using different colors, and we obtained and selected the data from the microscopic picture of the composite material using MATLAB.Additional data points were supplemented by fitting the profile shapes.Thus, the predicted steady-state ablation morphology was more detailed.It is also possible to write the profile shape functions directly into the MATLAB program.
The second step was to determine the inner normal direction of each point on the cross-sectional profile shape.Owing to the sufficiently dense data points, the inner normal direction n of point P i can be computed using the surrounding data points P i−1 and P i+1 .As illustrated in Figure 6, the signs of n x and n y must be determined based on the sizes of x i−1 and x i+1 , as well as the sizes of y i−1 and y i+1 .As the data are stored clockwise, when x i−1 < x i+1 and y i−1 < y i+1 , the sign of n is (+, −).When x i−1 < x i+1 and y i−1 > y i+1 , the sign of n is (−, −).When x i−1 > x i+1 and y i−1 > y i+1 , the sign of n is (−, +).When x i−1 > x i+1 and y i−1 < y i+1 , the sign of n is (+, +).When x i−1 < x i+1 and y i−1 = y i+1 , the value of n is (0, −1).When x i−1 = x i+1 and y i−1 > y i+1 , the value of n is (−1, 0).When x i−1 > x i+1 and y i−1 = y i+1 , the value of n is (0, 1).When x i−1 = x i+1 and y i−1 < y i+1 , the value of n is (1, 0).The unit normal direction (n x , n y ) of point P i can be expressed as: If the profile shape is written directly into the MATLAB program with functions, concave singularities may exist and need to be determined.The two extreme inner normal directions of the concave singularity can be calculated according to lines P i−1 P i and P i P i+1 .The remaining directions are selected every 1 • between the two extreme inner normal directions.The third step was to calculate the characteristic length lc of each point.The center of the characteristic circle is the designated turning point.When the distance between point Pi and the center of the circle ld < lc, then ld is the smallest distance from the center of this circle to each point of the cross-section.This dichotomy is a suitable method for solving the characteristic length lc of each point.The characteristic length lc must be in the range [ld,min, ld,max], where ld,min = 0, ld,max = Lmax, and Lmax is the maximum geometric length of the cross-section.Clearly, Lmax is larger than the distances from the center of this largest circle to the other points of the cross-section, and it is necessary to determine whether (ld,min + ld,max)/2 is the smallest of the distances between the center of this middle-sized circle and the points of the cross-section.When (ld,min + ld,max)/2 is the smallest length, ld,min = (ld,min + ld,max)/2; otherwise, ld,max = (ld,min + ld,max)/2.This processshould be repeated until ld,min -ld,max ≤ 10 -12 m.Then, lc = (ld,min + ld,max)/2.
The fourth step was to reconstruct the steady-state ablation morphology using a theoretical model.According to Equation (8), the relationship between lc -l(z) and z, which can also represent a typical steady-state ablation morphology, can be obtained from the relative reactivity η and the relative parameter of diffusion and reactivity χ.Accurate coordinates (xz, yz, zz) of the steady-state ablation morphology can be obtained by truncating the typical steady-state ablation morphology according to the characteristic length lc of each point, and placing it along the normal direction n on the coordinates (x0, y0) of each point.That is, the accurate coordinates (xz, yz, zz) of the steady-state ablation morphology can be obtained using the coordinate transformation: The final step was to plot the steady-state ablation morphology.

Method Validation
To verify the correctness of our model, we analyzed fiber morphology after oxidation at 1200 °C for 15 min by scanning electron microscope (SEM), as shown in Figure 7a.When the material and environmental parameters were η = 600 and χ = 3.5 × 10 -7 m, the calculation results for the fiber were in good agreement with the experimental results, as shown in Figure 7b.The fiber heights of the experimental and calculated results were 29.83 and 29.97 μm, respectively.The error of fiber height was approximately 0.46%.The fiber contour areas of the experimental and calculated results were 95.11 and 87.47 μm 2 , respectively.The error of fiber contour areas was approximately 8%, representing a verification The third step was to calculate the characteristic length l c of each point.The center of the characteristic circle is the designated turning point.When the distance between point P i and the center of the circle l d < l c , then l d is the smallest distance from the center of this circle to each point of the cross-section.This dichotomy is a suitable method for solving the characteristic length l c of each point.The characteristic length l c must be in the range [l d,min , l d,max ], where l d,min = 0, l d,max = L max , and L max is the maximum geometric length of the cross-section.Clearly, L max is larger than the distances from the center of this largest circle to the other points of the cross-section, and it is necessary to determine whether (l d,min + l d,max )/2 is the smallest of the distances between the center of this middle-sized circle and the points of the cross-section.When (l d,min + l d,max )/2 is the smallest length, l d,min = (l d,min + l d,max )/2; otherwise, l d,max = (l d,min + l d,max )/2.This processshould be repeated until l d,min − l d,max ≤ 10 −12 m.Then, l c = (l d,min + l d,max )/2.
The fourth step was to reconstruct the steady-state ablation morphology using a theoretical model.According to Equation (8), the relationship between l c − l(z) and z, which can also represent a typical steady-state ablation morphology, can be obtained from the relative reactivity η and the relative parameter of diffusion and reactivity χ.Accurate coordinates (x z , y z , z z ) of the steady-state ablation morphology can be obtained by truncating the typical steady-state ablation morphology according to the characteristic length l c of each point, and placing it along the normal direction n on the coordinates (x 0 , y 0 ) of each point.That is, the accurate coordinates (x z , y z , z z ) of the steady-state ablation morphology can be obtained using the coordinate transformation: The final step was to plot the steady-state ablation morphology.

Method Validation
To verify the correctness of our model, we analyzed fiber morphology after oxidation at 1200 • C for 15 min by scanning electron microscope (SEM), as shown in Figure 7a.When the material and environmental parameters were η = 600 and χ = 3.5 × 10 −7 m, the calculation results for the fiber were in good agreement with the experimental results, as shown in Figure 7b.The fiber heights of the experimental and calculated results were 29.83 and 29.97 µm, respectively.The error of fiber height was approximately 0.46%.The fiber contour areas of the experimental and calculated results were 95.11 and 87.47 µm 2 , respectively.The error of fiber contour areas was approximately 8%, representing a verification of the 2D results.There were several potential causes of error: (1) It is impossible for the fibers to be completely consistent from top to bottom; (2) SEM results can only extract two morphologies, which may be composed of the vertices of other point morphologies, and the left and right morphologies may have been truncated; (3) there were errors in fiber contour extraction; (4) There was deviation in parameter selection; etc.
of the 2D results.There were several potential causes of error: (1) It is impossible for the fibers to be completely consistent from top to bottom; (2) SEM results can only extract two morphologies, which may be composed of the vertices of other point morphologies, and the left and right morphologies may have been truncated; (3) there were errors in fiber contour extraction; (4) There was deviation in parameter selection; etc.When the material and environmental parameters were η = 20 and χ = 5 × 10 -3 m, the calculation results for the concave complex-fiber cross section are shown in Figure 8b.When compared with the calculation results of the Brownian motion simulation, as shown in Figure 8a, it can be seen that the calculation method in this study produced accurate results.It is worth pointing out that the number of fiber data points in the calculation result was 3321 × 251.The computational time using the calculation method proposed in this study was only approximately 58.55 s, approximately 1475 times faster than the Brownian motion simulation method proposed by Lachaud (approximately 24 h [43]).This is because the calculation method proposed in our study only needs to calculate one line and place it by coordinate transformation, without gradually evolving from the 3D initial state to the steady state over time.The computational complexity of the calculation method proposed in this study is approximately N, much smaller than that of the Brownian motion simulation method (approximately N 4 ).When the material and environmental parameters were η = 20 and χ = 5 × 10 −3 m, the calculation results for the concave complex-fiber cross section are shown in Figure 8b.When compared with the calculation results of the Brownian motion simulation, as shown in Figure 8a, it can be seen that the calculation method in this study produced accurate results.It is worth pointing out that the number of fiber data points in the calculation result was 3321 × 251.The computational time using the calculation method proposed in this study was only approximately 58.55 s, approximately 1475 times faster than the Brownian motion simulation method proposed by Lachaud (approximately 24 h [43]).This is because the calculation method proposed in our study only needs to calculate one line and place it by coordinate transformation, without gradually evolving from the 3D initial state to the steady state over time.The computational complexity of the calculation method proposed in this study is approximately N, much smaller than that of the Brownian motion simulation method (approximately N 4 ).
of the 2D results.There were several potential causes of error: (1) It is impossible for the fibers to be completely consistent from top to bottom; (2) SEM results can only extract two morphologies, which may be composed of the vertices of other point morphologies, and the left and right morphologies may have been truncated; (3) there were errors in fiber contour extraction; (4) There was deviation in parameter selection; etc.When the material and environmental parameters were η = 20 and χ = 5 × 10 -3 m, the calculation results for the concave complex-fiber cross section are shown in Figure 8b.When compared with the calculation results of the Brownian motion simulation, as shown in Figure 8a, it can be seen that the calculation method in this study produced accurate results.It is worth pointing out that the number of fiber data points in the calculation result was 3321 × 251.The computational time using the calculation method proposed in this study was only approximately 58.55 s, approximately 1475 times faster than the Brownian motion simulation method proposed by Lachaud (approximately 24 h [43]).This is because the calculation method proposed in our study only needs to calculate one line and place it by coordinate transformation, without gradually evolving from the 3D initial state to the steady state over time.The computational complexity of the calculation method proposed in this study is approximately N, much smaller than that of the Brownian motion simulation method (approximately N 4 ).Table 1 shows that the number of data points had little effect on the calculation results.The maximum height of the fiber was only affected by the number of data points in the XY plane.The reason is that the maximum characteristic length l c,max is independent of the number of data points on the Z axis, while the maximum height is calculated according to the maximum characteristic length and Equation (8).Roughness was less affected by the number of XY plane and Z axis data points, at approximately 0.24%.

Analysis of Steady-State Ablation Morphology
The temperature field can be considered uniform at the microscopic scale.The effect of temperature on the steady-state ablation morphology can be considered through the influence of the relative parameter of diffusion and reactivity χ, and the relative reactivity of the material η, which are defined in Section 2.2.As shown in Figure 9a-f, the steady-state ablation morphology of the fiber changed from a needle shape to a transition shape and then a platform shape.The height of the fiber gradually decreased when χ or η decreased.When the relative parameter of diffusion and reactivity χ was large, the gas molecules reached deeper below the surface and reacted with the interphase, thus forming a needle shape.When the relative reactivity of the material η was small, the reactivity of the fiber and interphase was similar, thus forming a platform shape.In summary, the needle shape was found to be controlled by the reaction rate, the shape of the platform was controlled by gas diffusion, while the transition shape was simultaneously controlled by gas diffusion and reaction rate.In addition, the characteristic length l c affected the steady-state ablation morphology.Steady-state ablation morphologies with different radii are shown in Figure 9g.The results showed that different characteristic lengths l c of the point corresponded to the morphology of different segments, consistent with the assumptions of previous modelling.As the characteristic length l c increased, the steady-state ablation morphology of the fiber also changed from a needle shape to a transition shape then a platform shape.This result explains why the fibers and bundles had needle and platform shapes, respectively.
As shown in Figure 10, the effect of different parameters on the steady-state ablation morphology of rectangular, hexagonal, and elliptical cross sections was similar to that for the circular fiber cross sections.Based on the analysis of Figures 4c and 10c,f,i, it is reasonable to define the characteristic length l c .The maximum characteristic length l c,max was 2 × 10 −4 m, indicated in Figure 10.The maximum heights of the different crosssections were the same: z max = 9.83 × 10 −4 m in Figure 10a,d,g, and z max = 7.50 × 10 −5 m in Figure 10b,e,h.Thus, the maximum characteristic length l c,max of the cross section determined the morphology height when the material and environmental parameters were the same.Fiber bundles of braided composites are likely to present rectangular and elliptical cross sections.The platform shapes (Figure 10b,e,h) were common to the steady-state ablation morphologies of the observed fiber bundles.

Steady-State Ablation Morphology of C/C Composites
The reconstruction of the steady-state ablation morphology of carbon fibers has been described previously.Carbon matrix can also be reconstructed in the same way as carbon fibers.Figure 11a shows the SEM morphology of the C/C composites [31].A simplified section of the fibers and matrix was established by observing the red rectangular area of the SEM morphology, as shown in Figure 11b.Figure 11c,d shows the C/C composite morphology and its equipotential map, according to calculation.In this model, the material and environmental parameters were the same as those described in Section 3.1: ηf = ηm = 20, χf = χm = 5 × 10 -3 m.The calculation results were similar to the SEM morphology of C/C composites.The calculation time was 9.23 s.The characteristic lengths lc of the fibers and matrix were in the ranges of [2.60 × 10 -6 , 2.95 × 10 -6 ] m and [4.78 × 10 -7 , 2.80 × 10 -6 ] m, respectively.The numbers of fibers and matrix data reflected in the calculation results were 2527 × 31 and 3658 × 31, respectively.The maximum heights of the fibers and matrix were 5.86 × 10 -5 m and 5.57 × 10 -5 m, respectively.The morphology of the matrix showed obvious undulations because the characteristic lengths lc of the points in the matrix were different.Owing to the limited depth of field of SEM, only the red and yellow parts may be observed.However, according to the analysis of the steady-state ablation morphology (the maximum characteristic lengths and heights of the fibers and matrix are approximate), the reactivity of the fiber and matrix was approximated by observing the SEM morphology.Because the fibers were needle-shaped, the steady-state ablation morphology was limited by the reaction.In addition, according to Equation ( 9), the roughness Sa = 2.00 × 10 -5 m.However, roughness is greatly affected by the material and environmental parameters.Thus, the roughness results can only be used as a reference.It is worth noting that the prediction of steady-state ablation morphology can guide the determination of ablation surface roughness to a certain extent when the material and environmental parameters are sufficiently accurate.The calculation results can have a guiding effect on the aerodynamic analysis of aircraft.

Steady-State Ablation Morphology of C/C Composites
The reconstruction of the steady-state ablation morphology of carbon fibers has been described previously.Carbon matrix can also be reconstructed in the same way as carbon fibers.Figure 11a shows the SEM morphology of the C/C composites [31].A simplified section of the fibers and matrix was established by observing the red rectangular area of the SEM morphology, as shown in Figure 11b.Figure 11c,d shows the C/C composite morphology and its equipotential map, according to calculation.In this model, the material and environmental parameters were the same as those described in Section 3.1: η f = η m = 20, χ f = χ m = 5 × 10 −3 m.The calculation results were similar to the SEM morphology of C/C composites.The calculation time was 9.23 s.The characteristic lengths l c of the fibers and matrix were in the ranges of [2.60 × 10 −6 , 2.95 × 10 −6 ] m and [4.78 × 10 −7 , 2.80 × 10 −6 ] m, respectively.The numbers of fibers and matrix data reflected in the calculation results were 2527 × 31 and 3658 × 31, respectively.The maximum heights of the fibers and matrix were 5.86 × 10 −5 m and 5.57 × 10 −5 m, respectively.The morphology of the matrix showed obvious undulations because the characteristic lengths l c of the points in the matrix were different.Owing to the limited depth of field of SEM, only the red and yellow parts may be observed.However, according to the analysis of the steady-state ablation morphology (the maximum characteristic lengths and heights of the fibers and matrix are approximate), the reactivity of the fiber and matrix was approximated by observing the SEM morphology.Because the fibers were needle-shaped, the steady-state ablation morphology was limited by the reaction.In addition, according to Equation ( 9), the roughness Sa = 2.00 × 10 −5 m.However, roughness is greatly affected by the material and environmental parameters.Thus, the roughness results can only be used as a reference.It is worth noting that the prediction of steady-state ablation morphology can guide the determination of ablation surface roughness to a certain extent when the material and environmental parameters are sufficiently accurate.The calculation results can have a guiding effect on the aerodynamic analysis of aircraft.

Advantages and Limitations
The advantages of our method are as follows: (1) Our method is fast, directly calculating the steady-state ablation morphology according to the 2D cross section and theoretical model.Namely, our method does not need to calculate the evolution process from the initial state to the transient state to the steady state.(2) Our method is effective.Both the gas diffusion equation and the chemical reaction equation are considered in the steady-state ablation morphology equation.Errors of fiber height and fiber contour areas were found to be small.(3) Our method is simple, only needing to calculate one line, and placing it by coordinate transformation, without gradually evolving from the 3D initial state to the steady state over time.There is no boundary movement in the modeling, and the physical and chemical processes are considered in the steady-state ablation morphology equation.
The limitations of our method are as follows: (1) Certain more detailed physical and chemical processes have not been considered, such as the effect of gas flow erosion on the steady-state ablation morphology.(2) The material properties of the interphase are difficult to obtain.(3) It is difficult for fibers to maintain the same cross-section from top to bottom.
The innovations of this paper are as follows: (1) Lachaud's theoretical model was extended to a 3D morphology with an irregular cross section.(2) The characteristic length lc of the theoretical model has been innovatively defined in this study.(3) The errors of our method have been quantified by fiber morphology.
The prospects for the method and the scope of further research include the following: (1) Our method is suitable for gasification-type ablative materials, such as carbon-carbon composites and resin matrix composites.(2) Other ablation materials (such as C/SiC, C/ZrC, SiC-SiC composites) also form needle shapes.Our method can be extended to these

Advantages and Limitations
The advantages of our method are as follows: (1) Our method is fast, directly calculating the steady-state ablation morphology according to the 2D cross section and theoretical model.Namely, our method does not need to calculate the evolution process from the initial state to the transient state to the steady state.(2) Our method is effective.Both the gas diffusion equation and the chemical reaction equation are considered in the steady-state ablation morphology equation.Errors of fiber height and fiber contour areas were found to be small.(3) Our method is simple, only needing to calculate one line, and placing it by coordinate transformation, without gradually evolving from the 3D initial state to the steady state over time.There is no boundary movement in the modeling, and the physical and chemical processes are considered in the steady-state ablation morphology equation.
The limitations of our method are as follows: (1) Certain more detailed physical and chemical processes have not been considered, such as the effect of gas flow erosion on the steady-state ablation morphology.(2) The material properties of the interphase are difficult to obtain.(3) It is difficult for fibers to maintain the same cross-section from top to bottom.
The innovations of this paper are as follows: (1) Lachaud's theoretical model was extended to a 3D morphology with an irregular cross section.(2) The characteristic length l c of the theoretical model has been innovatively defined in this study.(3) The errors of our method have been quantified by fiber morphology.
The prospects for the method and the scope of further research include the following: (1) Our method is suitable for gasification-type ablative materials, such as carbon-carbon composites and resin matrix composites.(2) Other ablation materials (such as C/SiC, C/ZrC, SiC-SiC composites) also form needle shapes.Our method can be extended to these ceramic-based ablative materials.(3) More detailed physical and chemical processes can be considered in the ablation morphology equation, to reduce the calculation errors when determining the steady-state ablation morphology.(4) The material properties of the interphase need to be measured in detail.

Conclusions
Based on the steady-state morphology equation, a fast, effective, and simple calculation method for steady-state ablation morphology was proposed and verified.The main conclusions are as follows: (1) Lachaud's theoretical model was extended to a 3D morphology with an irregular cross section.The characteristic length l c of the theoretical model was innovatively defined.The calculation method was verified and the correctness of the characteristic length was proved.(2) The calculation method proposed in this paper was found to be much faster and more efficient than the previous methods.This is because this method directly calculates the steady-state ablation morphology according to the 2D cross section and the theoretical model.There is no need to evolve from a 3D initial state to a steady state.(3) In addition, this method can accurately describe the morphology from the needle shape to the platform shape, consistent with the experimental results for fibers and bundles.As the characteristic length increased, the steady-state ablation morphology of the fiber changed from a needle shape to a transition shape then a platform shape.(4) For composite application, we reconstructed the steady-state ablation morphology of C/C composites using the calculation method proposed in this paper.Our results were consistent with the experimental observations, useful for analyzing physical and chemical properties and surface roughness features.(5) Furthermore, this method can improve the understanding of the interaction between composites and their surrounding environments.The material property relationship of the composite components can be affected by ablation morphology.The surface roughness features can provide a basis for better understanding the effect of surface morphology on aerodynamic analysis.

Figure 2 .
Figure 2. SEM images of C/C composites (a) before oxidation; (b) after oxidation at 1200 °C in air for 5 min; (c,d) after oxidation at 1200 °C in air for 15 min.

Figure 2 .
Figure 2. SEM images of C/C composites (a) before oxidation; (b) after oxidation at 1200 • C in air for 5 min; (c,d) after oxidation at 1200 • C in air for 15 min.

Figure 3 .
Figure 3. Schematic diagram of composites' surface morphology evolution, (a) initial microstructure of the material; (b) steady-state ablation morphology; (c) velocity relationship of the steady-state ablation morphology.

Figure 3 .
Figure 3. Schematic diagram of composites' surface morphology evolution, (a) initial microstructure of the material; (b) steady-state ablation morphology; (c) velocity relationship of the steady-state ablation morphology.

Figure 4 .
Figure 4. Schematic diagram for the characteristic lengths and circles of each point: (a) factor limiting fiber morphology; (b) selection of the suitable characteristic lengths and circles; (c) discussion on two singular points; (d) characteristic circles and centers of a concave complex fiber cross section.

Figure 4 .
Figure 4. Schematic diagram for the characteristic lengths and circles of each point: (a) factor limiting fiber morphology; (b) selection of the suitable characteristic lengths and circles; (c) discussion on two singular points; (d) characteristic circles and centers of a concave complex fiber cross section.

Figure 5 .
Figure 5. Flow chart of the calculation process of the steady-state ablation morphology.

Figure 5 .
Figure 5. Flow chart of the calculation process of the steady-state ablation morphology.

Figure 7 .
Figure 7.Comparison between experimental results and calculation results: (a) SEM of C/C composites oxidized at 1200 °C in air for 15 min, (b) experimental results and calculation results.

Figure 8 .
Figure 8. Calculation results for the irregular fiber cross-section: (a) Concave complex-fiber cross section and calculation result by Brownian motion simulation method (reprinted/adapted with permission from Elsevier, 2009 [43]); (b) the result in this study.

Figure 7 .
Figure 7.Comparison between experimental results and calculation results: (a) SEM of C/C composites oxidized at 1200 • C in air for 15 min, (b) experimental results and calculation results.

Figure 7 .
Figure 7.Comparison between experimental results and calculation results: (a) SEM of C/C composites oxidized at 1200 °C in air for 15 min, (b) experimental results and calculation results.

Figure 8 .
Figure 8. Calculation results for the irregular fiber cross-section: (a) Concave complex-fiber cross section and calculation result by Brownian motion simulation method (reprinted/adapted with permission from Elsevier, 2009 [43]); (b) the result in this study.

Figure 8 .
Figure 8. Calculation results for the irregular fiber cross-section: (a) Concave complex-fiber cross section and calculation result by Brownian motion simulation method (reprinted/adapted with permission from Elsevier, 2009 [43]); (b) the result in this study.

Figure 10 .
Figure 10.Different cross-sections and their steady-state ablation morphology: (a-c) rectangular cross section; (d-f) hexagonal cross section; (g-i) elliptical cross section.

Figure 10 .
Figure 10.Different cross-sections and their steady-state ablation morphology: (a-c) rectangular cross section; (d-f) hexagonal cross section; (g-i) elliptical cross section.

Figure 10 .
Figure 10.Different cross-sections and their steady-state ablation morphology: (a-c) rectangular cross section; (d-f) hexagonal cross section; (g-i) elliptical cross section.

Figure 11 .
Figure 11.Comparison of experimental observation with the steady-state ablation morphology of C/C composites: (a) SEM morphology of C/C composite (reprinted/adapted with permission from Hindawi, 2017 [31]); (b) C/C composite section; (c) C/C composite morphology by calculation; (d) equipotential map of the C/C composite morphology by calculation.

Figure 11 .
Figure 11.Comparison of experimental observation with the steady-state ablation morphology of C/C composites: (a) SEM morphology of C/C composite (reprinted/adapted with permission from Hindawi, 2017 [31]); (b) C/C composite section; (c) C/C composite morphology by calculation; (d) equipotential map of the C/C composite morphology by calculation.

Table 1 .
Calculation results with different numbers of data points.