With over half a century of development, laser technology has matured significantly^[1]. High-power lasers, characterized by their high pointing accuracy and long propagation distances, are capable of delivering energy swiftly and precisely to distant targets^[2-3]. However, during the propagation of high-power lasers , the small beam diameter and high energy density often cause localized heating of the propagation medium, leading to thermal blooming^[4-5]. This phenomenon results in beam phase distortion, significantly degrading the output beam quality and the far-field beam energy concentration of laser systems^[6]. Wu^[7] analyzed the deterioration of the beam’s β factor in a closed region with different gas media under the effect of airflow alone using simulation analysis. By comparing the analysis results with experimental data, it was shown that both exhibited the same variation trend, with the difference between the two being an order of magnitude smaller than the actual result. Therefore, it is crucial to suppress the thermal blooming effect during the propagation of high-power lasers in internal channels. Researchers have explored and confirmed the positive effects of introducing gas curtains or low-absorption gases in internal channels to suppress the thermal blooming^[8-9]. Nevertheless, in practical engineering, it remains essential to accurately and rapidly assess the impact of the thermal blooming effect to provide data support for selecting appropriate thermal blooming suppression strategies.

Regarding the numerical simulation of the thermal blooming effect in high-power lasers, researchers have already achieved numerous significant results. In the 1990s, Johnson B^[10] developed a computer simulation program for laser propagation in absorbing media, which effectively analyzed steady-state thermal blooming effects but was limited in analyzing transient thermal effects of high-power lasers. In 2011, Sun Q^[11] from the National University of Defense Technology utilized the custom subroutine function of FLUENT to couple the optical field and flow field. By dividing the propagation area into multiple sub-regions for iterative solving, this method effectively determined the optical field at the laser system’s exit. In 2017, Zhu F Y^[12] from the Institute of Optics and Electronics, Chinese Academy of Sciences, employed a similar approach to analyze the impact of thermal blooming on phase differences at the exit under both non-blowing and blowing conditions, and experimentally validated the method's accuracy. However, these methods require complex secondary development of simulation software and separate meshing of the optical field for calculations, making them unsuitable for rapid assessment needs during the early stages of engineering design. In 2022, Hu P^[13] and colleagues from the Beijing Institute of Applied Physics and Computational Mathematics studied the effects of various factors on thermal effects in horizontal enclosed light-conducting pipes and introduced internal optical path analysis software. This software can perform the thermal blooming simulation task but requires substantial computational resources and the ability to rapidly assess the thermal blooming effect in complex three-dimensional propagation paths can still be improved. In 2024, Zhou R^[14] introduced a method for analyzing thermal blooming effect using grid-wise integration approach in their study. This approach, which leverages finite element analysis, significantly simplifies the process of analyzing thermal blooming. However, grid-wise integration approach itself may encounter mesh mismatch issues in complex regions. Therefore, designing a practical and straightforward analysis model that meets the needs for rapid assessment of the thermal blooming effect of high-power laser gas in engineering remains a crucial issue in the current engineering field.

To provide a reference solution to the aforementioned issues, an equivalent analysis method for thermal blooming effect in the optical path is proposed. This method utilizes finite element analysis for flow field calculations. By constructing an equivalent analysis model and dividing the regions with identical thermal blooming effect, the impact of thermal blooming effect in the internal channel on the beam phase can be determined. Compared to the existing methods, this model offers a simpler equivalent process, a more flexible data handling method, and eliminates the need for complex software development, making it highly practical and accurate. The analytical results of this method provide references for engineering design.

The thermal blooming effect in the internal channel light path occurs when some of the gas in the propagation path absorbs laser energy upon irradiation. The irradiated gas undergoes a localized temperature increase, altering its density. In the absence of external lateral wind interference within the internal channel, the gas flows under the influence of gravity once its density changes. Furthermore, the non-uniform density distribution leads to a non-uniform refractive index during beam propagation, resulting in a phase difference^[15].

Based on an analysis of the mechanisms behind thermal blooming, the process includes two main components: the analysis of density distribution resulting from the interaction between the fluid and high-power laser, and the optical path difference analysis as the beam traverses a non-uniform medium.

The density of a gas is influenced by pressure and temperature, and their relationship is described by the ideal gas law, expressed as follows:

$ PV = nRT \quad,$ (1)

where $ P $ represents the gas pressure; $ V $ is the gas volume; $ T $ denotes the gas temperature; $ n $ is the amount of substance; $ R $ is the molar gas constant, typically $ R $=8.31 J/(mol·K).

The gas density within the infinitesimal element is:

$ \rho = \frac{m}{V} = \frac{{M \cdot n}}{V} = M \cdot \frac{P}{{RT}}\quad, $ (2)

where $ \rho $ is the gas density, $ m $ is the infinitesimal element mass, $ M $ is the molar mass.

The gas flow inside the closed pipe is driven by natural convection, caused by non-uniform density distribution and gravity. Due to the relatively small temperature rise, the impact of density changes on volume variation is negligible. Therefore, the model is simplified to an incompressible laminar flow model^[16-17]. The laser irradiation energy is considered as a heat source in the spot area. Under these conditions, the density distribution in the airflow can be determined using the conservation equations of mass, momentum, and energy^[18].

The mass conservation equation states that the mass entering an infinitesimal element per unit time is equal to the increase in mass within that infinitesimal element:

$ \frac{{\partial \rho }}{{\partial t}} + \nabla (\rho U) = 0\quad, $ (3)

where $ U $ is the velocity vector; in incompressible flow, $ \nabla (\rho U) = 0 $.

The momentum conservation equation states that the rate of change of momentum is equal to the net force acting on the infinitesimal volume element:

$ \frac{{\partial U}}{{\partial t}} + U \cdot \nabla U = - \frac{1}{\rho }\nabla p + G \quad,$ (4)

where $ G $ is the gravitational force, $\nabla p $ is the pressure Gradient of the fluid.

The energy conservation equation states that the increase in energy of an infinitesimal volume element is equal to the work done on it by the energy entering the element.

$ \frac{{\partial (\rho {T_{}})}}{{\partial t}} + \nabla \cdot (\rho U{T_{}}) = \nabla \cdot \left(\frac{k}{{\rho c}}\nabla T\right) + \frac{Q}{{\rho c}} \quad,$ (5)

where $ k $ is the thermal conductivity of the gas; $ c $ is the specific heat capacity of the gas; $ Q $ is the heat entering the infinitesimal element.

In the analysis of the heating effect of laser on gas, the energy Q entering the gas infinitesimal volume element is related to the gas absorption rate and the laser intensity distribution. The absorbed energy $ Q(x,y,\textit{z}) $ of the infinitesimal element at point $ (x,y,\textit{z}) $ in space can be expressed as:

$ Q(x,y,\textit{z}) = t \cdot P(x,y,\textit{z}) = t \cdot \eta \cdot I(x,y,\textit{z})\quad, $ (6)

where $ P(x,y,\textit{z}) $ is the power of the energy absorbed by the infinitesimal volume element at the point $ (x,y,\textit{z}) $ in space; $ t $ is the duration; $ \eta $ is the absorption rate of the gas for the specific wavelength of the laser; $ I(x,y,\textit{z}) $ is the intensity distribution at that point.

Using the finite element analysis method, the medium region is divided into multiple infinitesimal elements. By applying the aforementioned equations, the temperature change of the gas within each infinitesimal volume element after heating can be determined. Consequently, the density change within each infinitesimal volume element can be calculated.

After determining the density within each infinitesimal element, it is necessary to evaluate the specific impact of thermal blooming on beam quality deterioration through the beam’s phase difference. For the analysis of thermal blooming effect in complex internal channels, this study employs an equivalent substitution method to simplify the calculation process. The equivalent modeling procedure is as follows.

Assume the L-type light path within a spatial coordinate system, containing an equivalent mirror, incident light, and reflected light. In the L-type light path, the spatial coordinates of the center of any infinitesimal volume element are $ ({x_i},{y_i},{\textit{z}_i}) $, the gas density within this infinitesimal volume element is $ {\rho _h} $, and the coordinates of any three points on the equivalent plane of the mirror are $ {P_1}({x_1},{y_1},{\textit{z}_1}) $, $ {P_2}({x_2},{y_2},{\textit{z}_2}) $, and $ {P_3}({x_3},{y_3},{\textit{z}_3}) $. The positions are illustrated in Fig. 1.

The two vectors $ {\vec v_1} $ and $ {\vec v_2} $ on the equivalent plane of the mirror can be expressed as:

$ \left\{ \begin{gathered} {{\vec v}_1} = {P_1} - {P_2} = ({x_1} - {x_2},{y_1} - {y_2},{\textit{z}_1} - {\textit{z}_2}) \\ {{\vec v}_2} = {P_3} - {P_2} = ({x_3} - {x_2},{y_3} - {y_2},{\textit{z}_3} - {\textit{z}_2}) \\ \end{gathered} \right.\quad. $ (7)

The normal vector $ \vec n $ of the equivalent plane of the mirror is expressed as:

$ \vec n = {\vec v_1} \times {\vec v_2} \quad.$ (8)

The coordinates of the selected point are $ h(x,y,\textit{z}) $. The vector from $ h(x,y,\textit{z}) $ to $ {P_1}({x_1},{y_1},{\textit{z}_1}) $, which can be determined based on the coordinates of the two points, is denoted as $ \vec m = (x - {x_1},y - {y_1}, \textit{z} - {\textit{z}_1}) $. According to the spatial relationship between vectors, the mirror image point $ h'(x',y',\textit{z} ') $ of the selected point $ h(x,y,\textit{z}) $, formed by reflection across the mirror plane, can be calculated. The relationship between $h(x,y,\textit{z}) $ and $h'(x',y',\textit{z} ') $ is as follows:

$ (x',y',\textit{z} ') = (x,y,\textit{z}) - \frac{{2(\vec m \cdot \vec n)}}{{\vec n \cdot \vec n}}\vec n\quad.$ (9)

After the equivalent mirroring, the density at the equivalent point $ h' $ of $ h $ is $ {\rho _{h'}} = {\rho _h} $. The original L-type light path can be equivalently represented as a straight light path, as shown in Fig. 2.

The perturbation method is the simplest and most efficient approach commonly used for calculating the thermal blooming effect. This method treats the distortion caused by thermal blooming as a perturbation term added to the original undistorted beam. As a geometric optics method, it ignores the effect of diffraction^[19]. Given the short propagation distance and minimal external influences in the internal light path, this method been as the theoretical basis for evaluating the phase distortion induced by thermal blooming effect within the internal light path.

Based on the equivalent straight light path, the equivalent phase distribution of the beam affected by the thermal blooming effect as it traverses the system is calculated. The coordinate system is oriented such that the Z-axis aligns with the direction of light propagation. In the cross-section perpendicular to the light propagation direction, the area is divided into m units along both the X-axis and Y-axis, creating m² regions within the equivalent flow field. When these regions are sufficiently small, the thermal blooming effect is assumed to be uniformly distributed within each region, referred to as a region with identical thermal blooming effect. The side lengths of each region with identical thermal blooming effect in the X-axis and Y-axis directions are Δx and Δy, respectively, while the unit thickness corresponds to the light propagation distance. This subdivision process is depicted in Fig. 3.

Each region with identical thermal blooming effects contains H infinitesimal volume element . The average density $ {\rho _{{\text{average}}}} $ within each region with identical thermal blooming effect is calculated as follows:

$ {\rho _{{\text{average}}}} = \frac{{\displaystyle\sum\limits_{i = 1}^H {{\rho _i}} }}{H}\quad, $ (10)

where $ {\rho _i} $ is the density of an infinitesimal volume element within the region with identical thermal blooming effect.

According to the Gladstone-Dale relation ^[7,20] the average refractive index $ {n_{{\text{average}}}} $ within each region with identical thermal blooming effect can be calculated as follows:

$ {n_{{\text{average}}}} = K \cdot {\rho _{{\text{average}}}} + 1 \quad,$ (11)

$ K\approx \left({n}_{\text{s}}-1\right)\left(1+\frac{\alpha T}{P}\right)\quad, $ (12)

where $ K $ is the Gladstone-Dale constant, $ \alpha $ = 1/273.15 K, the refractive index $ {n_{\text{s}}} $ can be obtained through measurements under standard conditions ($ T $ = 293.15 K, $ P $ = 1 atm).

The phase difference $ \Delta {\varphi _i} $ induced by the beam passing through the regions with identical thermal blooming effect is expressed as:

$ \Delta {\varphi _i} = \frac{{2{\text{π}} }}{\lambda }OPD = \frac{{2{\text{π}} }}{\lambda }({n_{{\text{average}}}} - {n_s})L\quad, $ (13)

where OPD is the optical path difference, $ L $ is the beam propagation distance, and $ \lambda $ is the wavelength of the beam, $ {n_s} $ is the refractive index of the gas under standard conditions.

After determining the phase difference for each region with identical thermal blooming effect, the spatial position distribution of these regions can be used to obtain the equivalent analysis results of the phase difference distribution caused by thermal blooming effect during beam propagation.

The laser propagation channel analyzed in this study is a serpentine path, with the simplified model of the optical path shown in Fig. 4 (color online). The optical path includes three 45° mirrors and two 22.5° mirrors. Gravity acts in the negative Y-axis direction on the medium within the optical path. Due to the influence of gravity, the light path can be divided into two horizontal light-transmitting regions and two vertical light-transmitting regions.

The beam is incident perpendicular to surface A and exits perpendicular to surface B after five reflections. The beam diameter is 150 mm, and the light propagation channel diameter is 250 mm. The type of beam used in this study is a Gaussian flat-top beam, with the cross-sectional intensity distribution given by the following expression:

$ I(x,y) = {I_o}\exp \left[ - \left(\frac{{{x^2} + {y^2}}}{{w_{o}^2}} \right)^N\right] ,$ (14)

where, $ N $ is the order of expansion, $ {w_0} $ is the beam waist radius of the flat-top beam spot, and $ {I_o} $ is the peak light intensity.

In the study, the wavelength of the light beam is 1064 nm. The distribution of the gas absorption power in the cross-section, obtained by multiplying the gas absorption coefficient $ \eta $ with the light intensity $ {I_{}} $, is shown in Fig. 5 (color online). The power density at the flat-top region is 400 W/m³. The beam propagation time in this study is 3 seconds.

In this study, the finite element method (FEM) was employed for flow field calculations. COMSOL Multiphysics was used to discretize and analyze the flow field in detail.

After simplifying the flow field region, the structure becomes straight forward, and the flow field mesh is divided using a structured grid. The structured mesh can effectively enhance the speed of data storage and access. The boundary regions, including the fluid area in contact with the pipe wall and the lens, exhibited significant variations in medium properties. Therefore, eight layers of boundary layer mesh were added, with each layer having a thickness of D/85, where D represents the diameter of the light propagation region. Given the substantial temperature variations in the irradiated area at the center of the light propagation path, this region's mesh was refined to improve the accuracy of medium temperature rise calculations. The overall mesh element growth rate was maintained below 1.1. High-order elements were used in the calculations to improve the accuracy of the sampled node data in representing the overall simulation data. The total number of mesh elements is 1,851,075. The mesh element division model is shown in Fig. 6.

Due to the short laser irradiation time, the air density near the tube wall in the flow field changes minimally during laser loading, and there is almost no gas flow, the contact regions between medium and the pipe wall are set as no-slip boundaries. Since the medium in the light propagation path is influenced only by gravity and the gas flow speed is slow, the flow field simulation uses an incompressible flow. Considering that the gas under study is at normal temperature and pressure with low viscosity, the $ k - \varepsilon $ model is employed for simulation. Given that the analyzed environment is static, with no air movement in the flow field, the flow velocities at surfaces A and B of the simplified light propagation pipeline model are set to 0 m/s.

In this study, the propagation medium used is dry air at standard atmospheric pressure. The material properties are shown in Tab.1. The initial temperature of the environment is set to 300 K. After 3 seconds of light propagation, the temperature distribution at the center cross-section under no-wind condition is shown in Fig. 7 (color online). The temperature rise occurs in the regions where the light paths intersect, with a maximum temperature increase of 2.163 K. In areas without overlapping light path, the temperature rise is approximately 1 K.

The density, which directly influences the refractive index, is initially set to 1.17 kg/m³ for the propagation medium within the optical path,as shown in Fig. 8(a) (color online). The density values at each node after 3 seconds of light propagation are calculated using finite element analysis. The density distribution on the central cross-section is shown in Fig. 8(b) (color online). The density significantly decreases in the light path coverage area due to the temperature effect, with the non-overlapping region of the light path experiencing a density reduction of approximately 0.005 kg/m³ and the overlapping region experiencing a reduction of 0.01 kg/m³. Additionally, it is observed that under the influence of gravity, the propagation medium with lower density exhibits a clear buoyancy effect.

The node density distribution obtained from the finite element analysis software is extracted and shown in Fig. 8(c) (color online). Based on Equation (9), the equivalent model of the refracted light path was established, as shown in Fig. 8(d) (color online).

In the equivalent light propagation region model, the division of regions with identical thermal blooming effect and the calculation of the average density and optical path difference within these regions were completed.

The specific process is as follows: regions with identical thermal blooming effect are divided in the cross-section perpendicular to the light propagation direction, with m axial units along both the X-axis and Y-axis. As the number of axial units m increases, the calculation accuracy improves, but the computation time also increases. The relationships between the number of axial units m and the calculated phase difference, as well as between the number of axial units m and the computation time, are analyzed. The results are shown in Fig. 9. As the number of axial units varies from 64 to 1064, the computation time exhibits a clear upward trend. When the number of axial units exceeds 256, the calculation results of phase difference stabilize, and the data show significant convergence.

In the final analysis, the number of axial elements along the X-axis and Y-axis was set to $ m = 256 $.

Based on Equation (11) and equation (13), the phase distribution at the output port is obtained, as shown in Fig. 10 (color online). From the phase distribution image, it can be observed that under the influence of gravity, the calculation result of the equivalent model proposed in this study exhibit a phase distribution with distinct double-fisheye characteristics. Due to the energy distribution characteristics，at the cross-section of the flat-top beam, the phase variation rapidly decreases at the edges of the beam. After 3 seconds of laser irradiation, the results indicate a phase distortion of 4.7568 μm and a phase difference peak-to-valley (PV) of 2.9504 μm.

Using a self-developed program validated by experiments^[12], the simplified model of the same light propagation path was analyzed. After analyzing the thermal blooming effect following 3 seconds of laser irradiation in a no-wind environment, the results indicate a phase distortion of 4.9364 μm and a phase difference PV of 3.1949 μm.

By comparing the calculation results of the two methods, the deviation in phase distortion magnitude is less than 3.6%, and the deviation in phase PV is less than 7.7%. The obtained normalized beam wavefront aberration is shown in Fig. 11 (color online).

The comparison of wavefront aberration clearly demonstrates that the results calculated by the analytical model proposed in this study are highly consistent with those from existing programs.

Zernike polynomials are a set of orthogonal polynomials that can accurately describe wavefront errors and aberrations on a unit circle^[21-22]. The phase distortion calculated from the equivalent model and the self-developed program were fitted using Zernike polynomials. Standard Zernike polynomials were selected for the fitting process. The comparative fitting results are shown in Fig. 12 (color online).

According to the comparative results in Fig. 12, the calculations from the equivalent model and the self-developed program show a high degree of similarity. The differences are mainly reflected in the coefficients of the first, fifth, and thirteenth terms. The first, fifth, and thirteenth terms of the standard Zernike polynomials correspond to the tilt, the defocus, and the spherical aberration, respectively^[23].

Except for the first coefficient, which shows a significant difference, the variation trends of the other coefficients are generally consistent. The coefficients that have a substantial impact on phase distortion are the defocus and spherical aberration terms. The first term of the standard Zernike polynomials primarily represents the overall phase shift, which does not affect the actual beam quality.

Since this study uses a method of calculating phase distortion through regions with identical thermal blooming effect for equivalence, the input data for Zernike fitting is smoother compared to the input data from program calculations. As a result, the Zernike coefficients are slightly smaller than the program calculation results. For the fifth term (defocus), the coefficients from the equivalent model and the program calculations are −2.95 $ \lambda $ and −3.84 $ \lambda $, respectively. For the thirteenth term (spherical aberration), the coefficients are 3.5 $ \lambda $ and 2.9 $ \lambda $, respectively.

The phase distortion types calculated by both methods are similar. Combined with the comparison results of phase distortion magnitude and PV, the calculation results of this method have a certain reference value.

The L-type unit is a common basic structure in optical paths. Analyzing the thermal blooming effect within L-type units aids in optimizing the overall performance of optical systems. Therefore, in addition to the comprehensive optical path validation mentioned earlier, this study also includes validation of the accuracy of the thermal blooming effect analysis results for L-type units.

The L-type unit used in this study is shown in the Fig. 13. The ratio of the beam diameter to the width of the optical path is 3/8. The beam is a Gaussian beam with an absorption power density of 25 W/m³. The light entering from A, is reflected, and exits from B, with the gravitational direction aligned with the direction of light incidence.

During the analysis, the finite element mesh discretization method remains consistent with that described earlier. Helium is used as the propagation medium, with its material properties listed in the Tab. 2.

During the analysis, the initial temperature is set to 293.15 K. After 5 seconds of beam irradiation, the temperature distribution at the cross-section is shown in the Fig. 14 (color online). The temperature rise is highest in the overlapping region of the optical path near the mirror, with gas flow phenomena occurring due to the influence of gravity. The temperature rise in the non-overlapping region of the optical path is approximately 0.4 K.

The phase distribution at the output aperture is analyzed using the method described earlier. In this study, standard Zernike polynomials were used for fitting, with the first 15 terms selected for evaluation. The first three Zernike polynomials, which do not affect beam quality, were excluded from the analysis. The resulting Zernike polynomial coefficients at the output aperture are shown in the Fig. 15.

In the standard Zernike polynomial, the primary terms affecting wavefront aberrations are the fourth term (astigmatism-A), the fifth term (defocus), the sixth term (astigmatism-B), the eighth term (coma-A), the ninth term (coma-B), and the thirteenth term (spherical). The variation of the main Zernike terms after 5 seconds of beam irradiation is shown in Fig. 16 (color online).

It is evident that, after 5 seconds of beam irradiation, the primary aberrations at the output aperture are defocus and spherical aberration.

An analysis of the thermal blooming effect on the same L-type region has also been conducted in the reference [14]. The magnitude of aberrations or overall wavefront aberration is commonly assessed using the absolute values of Zernike coefficients, a common practice in wavefront analysis and optical design. Therefore, the variations in the absolute values of defocus and spherical aberration over 5 seconds of beam irradiation were compared between the results of this study and those reported in the reference, as shown in the Fig. 17 (color online).

From the Fig. 17, it is evident that, over 5 seconds of beam irradiation, the variation patterns of defocus and spherical aberration from the simplified model proposed in this study are consistent with those in the reference. While specific numerical values differ due to variations in optical path distances, the results are sufficiently accurate to validate the simplified model's effectiveness in assessing thermal blooming effect in the L-type unit.

Gravity, as an input term in the momentum conservation equation, is a significant factor contributing to gas turbulence within internal channels^[24,25]. Different gravity directions result in varying thermal blooming effects. This study employs the aforementioned analytical methods to investigate the thermal blooming effects of gas in a straight circular internal channel under the influence of gravity in different directions.

The straight circular internal channel has a diameter of 250 mm. The distribution of beam energy density, the properties of the propagation medium within the internal channel, and the discretization method for the flow field region are the same as described previously.

The angle between the direction of gravity and the beam direction is denoted as $ \theta $, as shown in Fig. 18. The range of $ \theta $ varies from 0° to 180°. The angle $ \theta $ is taken at observation points of 0°, 30°, 60°, 90°, 120°, 150°, and 180°. The variation process of phase distortion with angle under the influence of thermal blooming effect is documented, as shown in Fig. 19.

Based on the analysis results shown in Fig. 19, it is evident that under thermal blooming, the phase distortion within the straight circular channel changes with the angle $ \theta $. Specifically, when 0°<$ \theta $<90°, as $ \theta $ increases, the phase distortion shows a displacement, with larger angle changes resulting in greater displacements. When 90°<$ \theta $<180°, as $ \theta $ increases, the displacement of the phase distortion gradually decreases until the phase distortion returns to the central position at $ \theta $=180°. The phase distortion distribution exhibits symmetry around $ \theta $=90° under different gravity directions. In the ranges 0°<$ \theta $<60° and 120°<$ \theta $<180°, the displacement of the phase distortion changes significantly with $ \theta $. Within these angle ranges, the phase distortion distribution is sensitive to the direction of gravity.

As the angle $ \theta $ changes, the component of gravity perpendicular to the optical axis undergoes significant variation. The relationship between the gravity component perpendicular to the optical axis and the phase distortion with changing $ \theta $ is shown in Fig. 20. The analysis results show that as the angle $ \theta $ changes, the magnitude of phase distortion increases from 3.5498 μm at $ \theta $=0° to 4.7568 μm at $ \theta $=90°, and then decreases back to 3.5121 μm at $ \theta $=180°. The phase distortion trend strongly correlates with the variation trend of the gravity component perpendicular to the optical axis as $ \theta $ changes. In the ranges 0°<$ \theta $<60°and 120°<$ \theta $<180°, both the phase distortion and the gravity component perpendicular to the optical axis exhibit high rates of change. In the range 60°<$ \theta $<120°, the phase distortion is more gradual, and the rate of change of the gravity component perpendicular to the optical axis is also lower.

In this work, a rapid simulation model based on the finite element method (FEM) is developed. This model is used to evaluate phase distortion caused by thermal blooming effect during the propagation of high-power lasers within an internal channel. Flow field analysis and micro-element division were first performed using FEM principles. Subsequently, an equivalent model was constructed to calculate the phase distortion induced by thermal blooming effect after 3 seconds of light propagation by a high-power laser with a beam diameter of 150 mm. The results show a distinct 'double fisheye' pattern in the phase distortion distribution, with a phase distortion magnitude of 4.7568 μm and a PV of 2.9604 μm. The primary contributors to phase distortion were identified as the defocus and the spherical aberration. The method was validated by comparing it with existing methods, finding a high degree of similarity in the calculation results between the two approaches. The types of phase distortion were consistent, and the deviation in phase distortion magnitude was less than 3.6%. For L-type units, both methods identified the same primary factors influencing aberrations and exhibited consistent trends in their variation. Additionally, this model was used to evaluate the impact of thermal blooming under different gravity directions. The results show a strong correlation between the magnitude of phase distortion and the component of gravity perpendicular to the optical axis. As the perpendicular component of gravity increases, the phase distortion exhibits centroid shift, with its magnitude increasing from 3.5498 μm to 4.7568 μm, and then decreasing back to 3.5121 μm. This model eliminates the need for developing and debugging complex analysis code, offering significant versatility and flexibility. It is particularly valuable in the early stages of engineering design for the rapid assessment of thermal blooming effect caused by laser propagation in light-propagation paths and for selecting appropriate control strategies.