Heat-fluid-solid coupling analysis method based on direct simulation Monte Carlo method

Abstract

The invention discloses a heat-fluid-solid coupling analysis method based on a direct simulation Monte Carlo method, which adopts an improved time loose coupling strategy combining large-step thermodynamic elastic bidirectional coupling and small-step aeroelastic bidirectional coupling, and introduces a DSMC method into an analysis process to realize data exchange between a fluid domain and a solid domain. And finally, the change of the temperature and the displacement of the pneumatic control surface of the reentry vehicle in a certain space along with the time and the transient response at some moments are given. The method provided by the invention has the advantages of relatively small calculation amount and flexible calculation point selection.

Description

Technical Field

[0001] This invention belongs to the field of aircraft technology, specifically relating to a thermal-fluid-structure interaction analysis method based on the direct simulation Monte Carlo method. Background Technology

[0002] Space reentry vehicles generally refer to aircraft that fly at a Mach number greater than 5, escape the atmosphere, and then re-enter the atmosphere. For example... Figure 1 Since most of the outer surfaces of spacecraft use thin-walled structures, the aeroelasticity problem caused by the coupling effect between aerodynamic forces and thin-walled structures during high-speed operation has always been one of the challenges in this field. In supersonic or hypersonic environments, the thermal stress caused by aerodynamic heating makes this problem even more complicated, and it is necessary to consider the influence of the coupling effect of thermal stress, aerodynamic forces and elastic forces on thin-walled structures.

[0003] The coupled solutions to thermo-aeroelastic problems include both spatial and temporal coupling. Spatial coupling includes surface coupling and volume coupling. In thermo-aeroelastic problems, spatial coupling exhibits both surface and volume coupling. For example, the coupling between structural deformation and aerodynamic / thermal forces is a typical two-way surface coupling, while the coupling between thermal response and thermal stress / deformation is a typical two-way volume coupling. Spatial coupling depends on the relationships between the various physical fields and is generally deterministic.

[0004] Temporal coupling strategies can be broadly categorized into two types: fully coupled and loosely coupled. The fully coupled temporal analysis scheme best reflects the coupling characteristics of real physical conditions. During the calculation progresses over time, the flow field and structure are processed simultaneously. For each time step, the flow field calculation is accompanied by corresponding structural heat transfer-stress-strain calculations, and the selected time step is the minimum value among all physical field characteristic times. Figure 2 As shown, this coupling scheme involves a computational burden comparable to solving simultaneous equations of different physical fields, and yields the most accurate results. However, the computational cost and time required for this coupling strategy are extremely high, and its practical application is still immature.

[0005] To address the computational burden of the fully coupled time strategy, the loosely coupled time strategy makes the following assumption: the flow field is considered transiently stable compared to the structural calculation time step. Calculations show that even with large boundary heat fluxes, the structural temperature field does not change significantly over a short time interval, and this change has minimal impact on the flow field. Based on this assumption, the coupling strategy can be simplified. Instead of using the flow field calculation time step as the coupling time step, the temperature field calculation time step is used as the outer time step of the coupling calculation, nesting small-time-step aeroelastic calculations within it, such as... Figure 3 As shown.

[0006] The rarefied gas effect refers to the phenomenon where, with increasing altitude, the atmosphere becomes increasingly thinner, its density gradually decreases, the molecular number density decreases, molecular collisions decrease, the continuous flow assumption fails, the tangential forces generated by molecular collisions with the spacecraft surface cannot be ignored, the proportion of viscous terms in various aerodynamic parameters increases, and the viscous effect of the airflow intensifies. It is generally believed that the rarefied gas effect occurs above an altitude of 60 km, at which point the wall shear stress and wall heat flux given by the NSF (Navier-Stokes-Fourier) equations fail. In the simulation of rarefied flows, the Direct Simulation Monte Carlo (DSMC) method proposed by GA Bird is one of the most widely used methods. For over half a century, this method has been widely applied in various fields and validated by experimental results.

[0007] With the rapid improvement of computing power, thermo-aeroelastic methods using CFD / CSD (computational fluid dynamics / computational structural dynamics) for structural analysis have been increasingly studied and used. Currently, research on applying the DSMC method to (thermal)aeroelasticity is not very extensive, and it mainly focuses on the vibration of microstructures. Summary of the Invention

[0008] To overcome the shortcomings of existing technologies, this invention provides a thermo-fluid-structure interaction (TFI) analysis method based on the direct simulation Monte Carlo method. It employs an improved time-loose coupling strategy combining large-step thermo-aeroelastic bidirectional coupling and small-step aeroelastic bidirectional coupling. Simultaneously, the DSMC method is introduced into the analysis process to achieve data exchange between the fluid and solid domains, enabling TFI bidirectional coupling in rarefied fluid domains. Finally, the temperature and displacement variations of the aerodynamic control surfaces of a space reentry vehicle over time, as well as the transient responses at certain moments, are presented. This invention offers the advantages of relatively low computational cost and flexible selection of computation points.

[0009] The technical solution adopted by this invention to solve its technical problem is as follows: Step 1: Initialize the finite element program calculation; Step 2: SPARTA calculates and outputs; Step 3: Iterative calculation; Step 4: Calculation of structural transient response; Step 5: TPS panel thermo-aeroelasticity simulation.

[0010] Preferably, step 1 specifically comprises: In the SPARTA program, the surface of a three-dimensional watertight solid is composed of multiple two-dimensional triangles. In the finite element simulation, the element used is uniformly set as a tetrahedron. After setting the material properties and meshing, perform the initial calculation, i.e. the first load step, and output the surface triangle coordinates and surface temperature to SPARTA. The output file name is to_spa.dat.

[0011] Preferably, step 2 specifically comprises: The file to_spa.dat has n rows and 12 columns, where n is the total number of triangles, and the 12 columns are the number, type, xyz coordinates of the three points, and temperature of each triangle. The type is a user-defined positive real number. SPARTA performs DSMC calculations using the three-dimensional watertight solid enclosed by the triangles and the temperature of each triangle as boundary conditions. The calculated surface pressure and heat flux density are written to the file to_apdl.dat, which also has n rows and 12 columns.

[0012] Preferably, step 3 specifically comprises: The finite element program uses the pressure and heat flux density in the SPARTA output file as loads to perform the second load step of the thermo-solid coupling calculation and outputs the to_spa.dat file. This process is repeated to achieve the two-way thermo-solid coupling calculation.

[0013] Preferably, step 4 specifically comprises: Assume the time step from step 1 to step 3 is in seconds. After the calculation is completed, all results are saved in the relevant files of the finite element program. The file names are file.db and file.rst, which represent the settings and result data, respectively. Repeat steps 1 to 3, but this time replace the finite element program initialization calculation of step 1 with reading the result of a certain load step, and then output to_spa.dat to SPARTA for DSMC calculation. After that, read the load output by SPARTA in the finite element program, restart the calculation from the next load step, and apply an additional load at this time. Then repeat steps 1 to 3, and so on, setting the time step to the millisecond level to realize transient response calculation.

[0014] Preferably, step 5 specifically comprises: Step 5-1: Select the TPS panel geometry model; Step 5-2: Parameter equivalence treatment of dimensionless dynamic pressure; The formula for calculating the dimensionless dynamic pressure on TPS wall panels is: (1) In the formula For local dynamic pressure, For the length of the board, The local Mach number, The formula for calculation is: (2) In the formula This represents the Young's modulus of the metal sheet at 300K. For the thickness of the metal plate, Poisson's ratio; The dimensionless dynamic pressure calculated from the incoming flow parameters The dimensionless dynamic pressure calculated from the local inflow parameters at the head of the TPS plate Equal, that is: (3) Assuming only oxygen and nitrogen molecules are involved in the calculation, with oxygen molecules accounting for 22% and nitrogen molecules accounting for 78%, and the temperature remaining unchanged before and after the conversion, the Young's modulus of the metal plate after the conversion becomes 0.0232558 times, or 1 / 43 times, while the Poisson's ratio remains unchanged at 0.32. In order to perform three-dimensional finite element simulation, the radiation layer and the isolation layer are given extremely small stiffness, and the Young's modulus of both are set to one-thousandth of the Young's modulus of the metal plate before the conversion, and the Poisson's ratio is set to 0.32. Natural frequency Determined by the system's stiffness and inertia, it is positively correlated with the square root of stiffness k and negatively correlated with the square root of mass m, i.e., equation (4); the density of the three-layer material is simultaneously reduced to 0.112 times the original, while the thickness remains unchanged; since heat capacity is positively correlated with density, the specific heat capacity of the three-layer material is increased to 8.955 times the original. (4) Step 5-3: Three-dimensional conversion of the two-dimensional geometric model; The two-dimensional plate is stretched by 10mm along the lateral direction (z-axis) to form a three-dimensional wall panel. The metal plate, i.e., the bottom layer, is divided into two parts along the neutral plane. In the finite element calculation, the two line segments at the front and rear ends of the neutral plane of the metal plate are set as simply supported constraints. The z-direction displacement of the eight planes of the TPS wall panel in the lateral direction (z-axis) is set to 0. In SPARTA, only the upper surface of the TPS wall panel is set to participate in the particle collision calculation. The other surfaces directly "disappear" when they come into contact with particles, without generating force or heat. Step 5-4: Calculation parameter settings; The DSMC calculation uses the VHS variable hard sphere model. The grid step size and time step size used in the simulation follow the physical requirements of DSMC implementation. The grid width is less than 0.3 times the local molecular mean free path, and the time step size is less than 0.2 times the molecular mean collision time. Adaptive grid refinement is used. When the calculation reaches steady state, there are no more than 400 simulated particles in each grid. The macroscopic variables of the flow field are obtained by statistically analyzing 1000 steps of microscopic information. SPARTA through The parameter controls the number of simulated particles; this parameter is the ratio of the number of actual particles to the number of simulated particles. (Setting...) In the formula The actual particle number density, i.e., the simulated particle number density used in the calculation, is... per cubic meter; The CLL (Collision-Limited Linear Collision) model is adopted, which has two fitness parameters: the energy fitness coefficient of the normal velocity component. Energy adaptation coefficient of tangential velocity component The two parameters are set as follows: (5) In large-step calculations, the time step is set to... The total calculation time is 200 seconds. A pressure load equivalent to the pressure exerted on the lower surface of the plate when it is undeformed is applied. In small-step calculations, set the time step. Select respectively Transient calculations of the structure with small step size were performed at three time points. At the beginning of the small step size calculation, a transient disturbance was applied to the plate. Specifically, an additional pressure of 150 Pa was applied to the lower surface for 0.005 s.

[0015] An electronic device includes: a processor and a memory; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to enable the electronic device to perform the above-described thermal-fluid-structure interaction analysis method.

[0016] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described thermal-fluid-structure interaction analysis method.

[0017] A chip includes a processor for calling and running a computer program from a memory, causing a device on which the chip is mounted to perform the above-described thermal-fluid-structure interaction analysis method.

[0018] A computer program product includes a computer storage medium storing a computer program, the computer program including instructions executable by at least one processor, which, when executed by the at least one processor, implement the above-described thermal-fluid-structure interaction analysis method.

[0019] The beneficial effects of this invention are as follows: 1) The computational cost of wing / rudder thermo-aeroelastic analysis of hypersonic vehicles based on CFD / CSD is large. The improved time loose coupling strategy proposed in this invention has the advantages of relatively small computational cost and flexible selection of computation points. 2) This invention introduces the DSMC method into thermo-aeroelastic simulation analysis, which solves the problem of inaccurate CFD calculation in thin water domains and realizes thermo-aeroelastic simulation in thin water domains. 3) This invention extends the 2D TPS thermal protection wall panel to 3D and performs equivalent processing, and calculates displacement amplitude response results similar to those in the literature, proving the effectiveness of this invention. Attached Figure Description

[0020] Figure 1 A schematic diagram of the leeward side of the control surfaces (left) and the side of the fuselage (right) during Starship's fifth test flight reentry; Figure 2 Flowchart of the time-coupling strategy; Figure 3 Flowchart of the time-loose coupling strategy; Figure 4 Flowchart for the improved time-loose coupling strategy; Figure 5 This is a schematic diagram of a spatial coupling strategy; Figure 6 Here are the geometric schematics for SOLID227 and SURF252 elements; Figure 7 Here is a flowchart of the program's cyclic calculation process; Figure 8 A schematic diagram of the TPS thermal protection panel structure on the inclined surface of a scramjet engine. Figure 9 shows the material properties as a function of temperature: (a) specific heat capacity of the TPS structure; (b) thermal conductivity of the TPS structure; (c) Young's modulus and coefficient of thermal expansion of the wall panel. Figure 10 A schematic diagram of a 3D TPS panel (left) and mesh generation (right); Figure 11 This is a schematic diagram showing the upper surface temperature of the three-layer plate at x = 0.5, 1, and 1.5. Figure 12 This is a schematic diagram showing the y-direction displacement at 3 / 4 of the wall panel. Figure 12 (a) represents a large step size. Figure 12 (b) is t=10s. Figure 12 (c) is t = 80s. Figure 12 (d) is t=190s; Figure 13 This is a schematic diagram of the y-direction displacement results at 3 / 4 of the TPS wall panel in the literature. Figure 13 (a) is for the overall view. Figure 13 (b) is t = 10s. Figure 13 (c) is t=50s. Figure 13 (d) is t=100s; Figure 14Screenshot of the coupled computation program running; Figure 15 This is a flow field mass density contour map at time t = 0s; Figure 16 The deformation of the plate and a molecular snapshot at t = 190s; Figure 17 This is a geometric schematic diagram of the computational model for the aerodynamic rudder of the Quasar. Figure 18 This is a schematic diagram showing the main dimensions of the aerodynamic rudder calculation model; Figure 19 shows the temperature cloud map of the aerodynamic rudder at t = 500s. Figure 19(a) shows the windward side and Figure 19(b) shows the leeward side. Figure 20 The aerodynamic rudder displacement contour plot at t = 500s; Figure 21 This is a schematic diagram of the trailing edge displacement response of the aerodynamic rudder tip. Figure 21 (a) shows the total displacement result for the large step size. Figure 21 (b) shows the transient calculation results at t=500s. Detailed Implementation

[0021] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0022] To address the failure of CFD methods in high-altitude sparse watersheds, the DSMC method is introduced into thermo-aeroelastic simulation analysis. At the same time, an improved time-loose coupling strategy is proposed, which effectively reduces the computational load and time of thermo-aeroelastic analysis while ensuring computational accuracy, thereby improving the flexibility of the analysis.

[0023] To address the computationally intensive nature of current CFD / CSD-based thermo-aeroelastic analysis of hypersonic vehicle wings / rudders, an improved time-loose coupling strategy combining large-step-length and small-step-length thermo-aeroelastic bidirectional coupling is proposed. Simultaneously, the DSMC method is introduced into the analysis flow to enable data exchange between the fluid and solid domains, facilitating thermo-fluid-structure interaction in the rarefied domain. Based on this invention, the temperature and displacement variations of the aerodynamic control surfaces of a space reentry vehicle over time, as well as their transient responses at certain moments, are presented.

[0024] The following are the main contents of the invention: The improved time-loosely coupled strategy makes two assumptions: 1) Similar to the assumptions of the loosely coupled time strategy mentioned above, it is assumed that the time scale of changes in the structural temperature field and thermal strain is large, much larger than the coupling effect (flutter) between structural vibration and flow field.

[0025] 2) The structure temperature continues to rise, but it remains stable without being disturbed and will not flutter. The structure deforms slowly as the temperature rises.

[0026] Based on the above assumptions, an improved time-loosened thermo-aeroelastic strategy is as follows: Figure 4 As shown. The first part of the improved time-loosely coupled strategy follows a similar process to the time-fully coupled strategy, employing a thermo-fluid-structure interaction (TFISE) bidirectional coupling strategy. The difference lies in the absence of nested small-step aeroelastic calculations. In the second part, after saving the structural state, including displacement and temperature, at each analysis moment, a small-step (millisecond-level) aeroelastic or thermo-aeroelastic bidirectional coupling analysis is performed using the temperature and displacement fields at a given moment as initial conditions. This yields the transient response of the wing to a small disturbance at that moment. Since the timescale of structural vibration is extremely small, the effects of aerodynamic heating and structural heat transfer can be disregarded. Furthermore, because these moments can be manually selected based on the situation, the computational load can be significantly reduced, increasing the flexibility of the analysis.

[0027] This invention is mainly achieved by coupling a general-purpose finite element calculation program and the SPARTA open-source DSMC calculation program using shell scripts under the Linux system. The overall spatial coupling strategy is as follows: Figure 5 As shown.

[0028] Step 1: Initialize the finite element program calculation; In the SPARTA program, the surface of a three-dimensional watertight solid is composed of multiple two-dimensional triangles. Therefore, in order to facilitate the transfer of loads between the two programs, the elements used in finite element simulation are uniformly set to tetrahedrons, such as the SOLID227 element in MAPDL. This element can perform thermo-solid bidirectional coupling calculations and consider thermal effects such as thermal expansion, thermal stress, and thermal strain. It can also generate SURF252 elements on the outer surface of the model for thermal radiation simulation.

[0029] After setting the material properties and meshing, perform the initial calculation, which is the first load step calculation. At this time, there are generally no loads on the model, only constraints. The purpose is simply to output the surface triangle coordinates and surface temperature to SPARTA. The output file name is to_spa.dat.

[0030] Step 2: SPARTA calculates and outputs; According to the requirements of the SPARTA program, the main part of the to_spa.dat file has n rows and 12 columns, where n is the total number of triangles, and the 12 columns are the number, type, xyz coordinates of the three points, and temperature of each triangle. The type is a user-defined positive real number, used to distinguish different surface types in SPARTA. SPARTA performs DSMC calculations using the three-dimensional watertight solid enclosed by triangles and the temperature of each triangle as boundary conditions. The calculated surface pressure and heat flux density are written to the to_apdl.dat file, which also has n rows and 12 columns, with the 12 columns being the number, xyz coordinates of the three points, pressure, and heat flux density of each triangle.

[0031] Step 3: Iterative calculation; The finite element program uses the pressure and heat flux density from the SPARTA output file as loads to perform the second load step of the thermo-mechanical coupling calculation and outputs the to_spa.dat file. This process is repeated to achieve the two-way thermo-mechanical coupling calculation. The specific implementation process of steps one to three is as follows: Figure 7 As shown in the diagram, steps ②, ③, and ④ are the loop parts. The control.sh file is the overall control file for automatic calculations and belongs to the Linux Shell script category. The xxx.mac file is the finite element program input script file, and in.xxx is the SPARTA input file.

[0032] Step 4: Calculation of structural transient response; Assuming the time step for the first three steps is in the second range (e.g., 1 second), after the calculation is complete, all results are saved in the relevant files of the finite element program, with filenames such as file.db and file.rst, representing the settings and results data, respectively. The calculation of the first three steps is then repeated, but this time the finite element program initialization calculation of the first step is replaced by reading the results of a specific load step (e.g., step 50), and then outputting to_spa.dat to SPARTA for DSMC calculation. Afterwards, the finite element program reads the load output from SPARTA and restarts the calculation from the next load step (e.g., step 51 at 51 seconds), applying an additional load (perturbation) at this time. Steps one through three are then repeated in a loop, with the time step set to the millisecond range (e.g., 0.001 seconds), thus realizing the transient response calculation of the model at a certain moment.

[0033] Step 5: TPS thermo-aeroelastic simulation of the wall panel and comparison with literature; 1) Selection of the geometric model for the TPS panel; The literature “Xie Dan, Ji Chunxiu, Jing Xingjian. Thermo-aeroelastic dynamics analysis of a wall panel under typical hypersonic trajectories [J]. Acta Aeronautica Sinica, 2021(11): 375-390” presents a two-dimensional thermo-aeroelastic coupled modeling and analysis of a two-dimensional wall panel with a thermal protection system (TPS) in hypersonic flow. The aerodynamic force is calculated using the third-order piston theory, the aerodynamic heat flow is obtained through the reference enthalpy method, and the structural heat conduction is calculated based on the finite difference method.

[0034] Figure 9 shows the temperature-dependent changes in the material parameters of each layer of the TPS wall panel, representing the specific heat capacity, thermal conductivity, Young's modulus, and coefficient of thermal expansion. In the figure, the blue line represents the metal wall panel material, the yellow line represents the insulation layer material, and the red line represents the radiating layer material. The literature simplifies the insulation and radiating layer materials, considering only their mass and ignoring stiffness; therefore, Young's modulus and coefficient of thermal expansion are not included for these two layers. Table 1 shows the parameters of each layer of the TPS wall panel. The Poisson's ratio for the metal plate is 3.2, and the surface of the radiating layer emits non-blackbody radiation with an emissivity of... .

[0035] Table 1 Parameters of each layer of the two-dimensional TPS wall panel

[0036] 2) Parameter equivalence treatment for dimensionless dynamic pressure; To compare with the literature, the Mach number from the literature was selected. ,high The calculations are performed under the following operating conditions. Since the inflow density under these conditions is too high for DSMC calculations, it is necessary to perform an equivalent transformation of the inflow parameters and wall panel material based on the dimensionless dynamic pressure exerted on the wall panel. The formula for calculating the dimensionless dynamic pressure exerted on the TPS wall panel is: (1) In the formula For local dynamic pressure, For the length of the board, The local Mach number, The formula for calculation is: (2) In the formula This represents the Young's modulus of the metal sheet at 300K. For the thickness of the metal plate, Let be Poisson's ratio. Let the dimensionless dynamic pressure calculated from the incoming flow parameters in this paper be equal to the dimensionless dynamic pressure calculated from the local incoming flow parameters at the head of the TPS plate in the literature, that is: (3) Assuming only oxygen and nitrogen molecules are involved in the calculation, with oxygen molecules accounting for 22% and nitrogen molecules accounting for 78%, and the temperature remaining unchanged before and after the conversion, the incoming flow parameters after the conversion are shown in Table 6. After the conversion, the Young's modulus of the metal plate becomes 0.0232558 times its original value, or 1 / 43 times, while the Poisson's ratio remains unchanged at 0.32. For the purpose of three-dimensional finite element simulation, the radiating layer and the insulating layer are given extremely low stiffness; their Young's modulus is set to one-thousandth of the original Young's modulus of the metal plate, and their Poisson's ratio is set to 0.32.

[0037] Table 2. Transformed Incoming Flow Parameters

[0038] Generally, the natural frequency is determined by the system's stiffness and inertia (mass or density), and is positively correlated with the square root of the stiffness and negatively correlated with the square root of the mass, i.e., equation (4). After the above treatment, the natural frequency of the wall panel decreases because the stiffness of the metal plate decreases while the mass remains unchanged. In order to match the natural frequency of the original wall panel in the literature, the density of the three layers of materials is simultaneously reduced to 0.112 times the original, while the thickness remains unchanged. Since heat capacity is positively correlated with density (mass), the specific heat capacity of the three layers of materials also needs to be increased to 8.955 times the original. The parameters of the wall panel after treatment are shown in Table 3, and the comparison of the natural frequencies of the wall panel before and after treatment at a temperature of 300K is shown in Table 4.

[0039] (4) Table 3 Parameters of each layer of the treated wall panel

[0040] Table 4 Comparison of natural frequencies of the wall panel before and after material density treatment

[0041] 3) Three-dimensional conversion of two-dimensional geometric models; To perform three-dimensional finite element calculations, a two-dimensional plate was stretched 10 mm laterally (along the z-axis) to form a three-dimensional wall panel with a very small width, such as... Figure 10 As shown, the metal plate (bottom layer) is divided in two along the neutral plane to set up simply supported constraints. In the finite element method (FEM) calculation, the two line segments at the front and rear ends (x=0 and x=1.5) of the metal plate's neutral plane are set as simply supported constraints, and the z-direction displacement of the eight lateral (z-axis) planes of the TPS panel is set to 0. In SPARTA, only the upper surface of the TPS panel is configured to participate in particle collision calculations; the other surfaces simply "disappear" upon contact with particles, generating no force or heat. This allows for the simulation of two-dimensional situations under three-dimensional conditions.

[0042] 4) Calculation parameter settings; DSMC calculations employ the VHS variable hard sphere model, which does not involve chemical reactions. The grid step size and time step size used in the simulation follow the physical requirements of DSMC implementation. The grid width is less than 0.3 times the local molecular mean free path, and the time step size is less than 0.2 times the molecular mean collision time. Adaptive mesh refinement is used; once the calculation reaches steady state, each grid contains no more than 400 simulated particles. Macroscopic variables of the flow field are obtained by statistically analyzing 1000 steps of microscopic information. SPARTA primarily utilizes... The parameter controls the number of simulated particles; it is the ratio of actual particles to simulated particles. A smaller value indicates more simulated particles are involved in the calculation. (Settings...) In the formula The actual particle number density, i.e., the simulated particle number density used in the calculation, is... Each cubic meter. The CLL surface collision model is adopted, which has two adaptation parameters: the energy adaptation coefficient of the normal velocity component. Energy adaptation coefficient of tangential velocity component The two parameters are set as follows: (5) In large-step calculations, the time step is set to... The calculation lasts for a total of 200 seconds. A pressure load equivalent to the pressure experienced by the undeformed upper surface is applied to the lower surface of the plate. In the small-step calculation, a time step is set. Select respectively Transient calculations of the structure with small step size were performed at three time points. At the beginning of the small step size calculation, a transient disturbance was applied to the plate, specifically by applying an additional pressure of 150 Pa on the lower surface for 0.005 s.

[0043] 5) Calculation results; Figure 11 The three-layer board surface is shown on x = Temperature changes at 0.5, 1, and 1.5. Figure 12 For the 3 / 4 section of the wall panel (x / a = 0.75) y The displacement results are shown in the figure, where (a) is the result calculated with a large step size, and (b), (c), and (d) are the transient responses at three time points obtained from the calculation with a small step size. w / h For dimensionless displacement, i.e., removal of the metal plate thickness. Figure 13 The TPS wall panel is given at 3 / 4 of its length. yDisplacement results. Comparison with literature results shows that the present invention can calculate the heating process of the TPS wall panel over time, and simulate the phenomenon that the amplitude of the displacement response of the wall panel under small disturbance changes from small to medium and then to large over time. The correspondence with the literature results is shown in Table 5, which to a certain extent proves the correctness and effectiveness of the present invention.

[0044] Table 5. Correspondence between displacement response results and literature results

[0045] Example: Figure 17 The aerodynamic control system consists of a hollow, framed rudder and a forebody baffle at the wing root. Its shape is modeled after the aft aerodynamic control system of a Starship, but approximately 1 / 10 the size. The square holes in the frame are for SPARTA calculations and to avoid creating cavities. The rudder has a wing root thickness of 60mm and a wingtip thickness of 30mm, with a linear transition from root to wingtip. The skin thickness is 3mm, and the frame thickness is 5mm. TPS thermal protection structures are installed on the windward surfaces of both the control surface and the forebody. The radiant layer is 2mm thick, and the insulating layer is 8mm thick. The materials for the radiant and insulating layers are treated the same as described above, but the metal plates are left untreated, using the original material properties from the literature. The control surface has a 30° dihedral angle, and the forebody baffle has a 45° dihedral angle. The main dimensions and incoming airflow direction are as follows. Figure 18 As shown in Table 3, the inflow parameters at an altitude of 65 km were used for the calculations, and the inflow velocity was referenced from the Starship's 10th test flight reentry data.

[0046] Table 6. Flow parameters at an altitude of 65 km

[0047] First, the time step is set to 1 second for large-step calculations. The calculation takes a total of 500 steps, or 500 seconds. Figure 19 shows... t Temperature cloud maps of the windward and leeward sides of the aerodynamic rudder at time 500s show that the highest temperature occurs at the very front of the baffle, at 1776K. The highest temperature on the rudder surface occurs immediately next to the baffle, at 1465K. Figure 20 Showing t =500s total displacement cloud map of the aerodynamic rudder, the maximum displacement is 17.7mm, which appears at the trailing edge of the wingtip.

[0048] An additional disturbance load was applied to the control surface at time t=500s. Specifically, an additional pressure of 500Pa was applied to the windward side for 0.005s. The displacement response at the wingtip trailing edge was then selected. To conduct analysis, Figure 21 (a) and (b) respectively show the total displacement results for the large step size and tTransient calculation results at time 500s. Figure 21 It can be concluded that the aerodynamic rudder undergoes aerodynamic heating and slow deformation over 500 seconds. t If disturbed at time 500s, a small displacement vibration will occur, and the vibration will show a convergence trend.

Claims

1. A thermo-fluid-structure interaction analysis method based on direct simulation Monte Carlo method, characterized in that, Includes the following steps: Step 1: Initialize the finite element program calculation; Step 2: SPARTA calculates and outputs; Step 3: Iterative calculation; Step 4: Calculation of structural transient response; Step 5: TPS panel thermo-aeroelasticity simulation.

2. The thermal-fluid-structure interaction analysis method based on direct simulation Monte Carlo method according to claim 1, characterized in that, Step 1 specifically involves: In the SPARTA program, the surface of a three-dimensional watertight solid is composed of multiple two-dimensional triangles. In the finite element simulation, the element used is uniformly set as a tetrahedron. After setting the material properties and meshing, perform the initial calculation, i.e. the first load step, and output the surface triangle coordinates and surface temperature to SPARTA. The output file name is to_spa.dat.

3. The thermal-fluid-structure interaction analysis method based on direct simulation Monte Carlo method according to claim 2, characterized in that, Step 2 specifically involves: The file to_spa.dat has n rows and 12 columns, where n is the total number of triangles, and the 12 columns are the number, type, xyz coordinates of the three points, and temperature of each triangle. The type is a user-defined positive real number. SPARTA performs DSMC calculations using the three-dimensional watertight solid enclosed by the triangles and the temperature of each triangle as boundary conditions. The calculated surface pressure and heat flux density are written to the file to_apdl.dat, which also has n rows and 12 columns.

4. The thermal-fluid-structure interaction analysis method based on direct simulation Monte Carlo method according to claim 3, characterized in that, Step 3 specifically involves: The finite element program uses the pressure and heat flux density in the SPARTA output file as loads to perform the second load step of thermo-mechanical coupling calculation and outputs the to_spa.dat file. This process is repeated to achieve bidirectional thermo-mechanical coupling calculation.

5. A thermo-fluid-structure interaction analysis method based on direct simulation Monte Carlo method according to claim 4, characterized in that, Step 4 specifically involves: Assume the time step from step 1 to step 3 is in seconds. After the calculation is completed, all results are saved in the relevant files of the finite element program. The file names are file.db and file.rst, which represent the settings and result data, respectively. Repeat steps 1 to 3, but this time replace the finite element program initialization calculation of step 1 with reading the result of a certain load step, and then output to_spa.dat to SPARTA for DSMC calculation. After that, read the load output by SPARTA in the finite element program, restart the calculation from the next load step, and apply an additional load at this time. Then repeat steps 1 to 3, and so on, setting the time step to the millisecond level to realize transient response calculation.

6. A thermo-fluid-structure interaction analysis method based on direct simulation Monte Carlo method according to claim 5, characterized in that, Step 5 specifically involves: Step 5-1: Select the TPS panel geometry model; Step 5-2: Parameter equivalence treatment of dimensionless dynamic pressure; The formula for calculating the dimensionless dynamic pressure on TPS wall panels is: （1） In the formula For local dynamic pressure, For the length of the board, The local Mach number, The formula for calculation is: （2） In the formula This represents the Young's modulus of the metal sheet at 300K. For the thickness of the metal plate, Poisson's ratio; The dimensionless dynamic pressure calculated from the incoming flow parameters The dimensionless dynamic pressure calculated from the local inflow parameters at the head of the TPS plate Equal, that is: （3） Assuming only oxygen and nitrogen molecules are involved in the calculation, with oxygen molecules accounting for 22% and nitrogen molecules accounting for 78%, and the temperature remaining unchanged before and after the conversion, the Young's modulus of the metal plate after the conversion becomes 0.0232558 times, or 1 / 43 times, while the Poisson's ratio remains unchanged at 0.

32. In order to perform three-dimensional finite element simulation, the radiation layer and the isolation layer are given extremely small stiffness, and the Young's modulus of both are set to one-thousandth of the Young's modulus of the metal plate before the conversion, and the Poisson's ratio is set to 0.

32. Natural frequency Determined by the system's stiffness and inertia, it is positively correlated with the square root of stiffness k and negatively correlated with the square root of mass m, i.e., equation (4); the density of the three-layer material is simultaneously reduced to 0.112 times the original, while the thickness remains unchanged; since heat capacity is positively correlated with density, the specific heat capacity of the three-layer material is increased to 8.955 times the original. （4） Step 5-3: Three-dimensional conversion of the two-dimensional geometric model; The two-dimensional plate is stretched by 10mm along the lateral direction (z-axis) to form a three-dimensional wall panel. The metal plate, i.e., the bottom layer, is divided into two parts along the neutral plane. In the finite element calculation, the two line segments at the front and rear ends of the neutral plane of the metal plate are set as simply supported constraints. The z-direction displacement of the eight planes of the TPS wall panel in the lateral direction (z-axis) is set to 0. In SPARTA, only the upper surface of the TPS wall panel is set to participate in the particle collision calculation. The other surfaces directly "disappear" when they come into contact with particles, without generating force or heat. Step 5-4: Calculation parameter settings; The DSMC calculation uses the VHS variable hard sphere model. The grid step size and time step size used in the simulation follow the physical requirements of DSMC implementation. The grid width is less than 0.3 times the local molecular mean free path, and the time step size is less than 0.2 times the molecular mean collision time. Adaptive mesh refinement is used, and once the calculation reaches steady state, there are no more than 400 simulated particles in each mesh. The macroscopic variables of the flow field are obtained by statistically analyzing 1000 steps of microscopic information. SPARTA through The parameter controls the number of simulated particles; this parameter is the ratio of the number of actual particles to the number of simulated particles. (Setting...) In the formula The actual particle number density, i.e., the simulated particle number density used in the calculation, is... per cubic meter; The CLL (Collision-Limited Linear Collision) model is adopted, which has two fitness parameters: the energy fitness coefficient of the normal velocity component. Energy adaptation coefficient of tangential velocity component The two parameters are set as follows: （5） In large-step calculations, the time step is set to... The total calculation time is 200 seconds. A pressure load equivalent to the pressure exerted on the lower surface of the plate when it is undeformed is applied. In small-step calculations, set the time step. Select respectively Transient calculations of the structure with small step size were performed at three time points. At the beginning of the small step size calculation, a transient disturbance was applied to the plate. Specifically, an additional pressure of 150 Pa was applied to the lower surface for 0.005 s.

7. An electronic device, characterized in that, include: Processor and memory; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to cause the electronic device to perform the method as described in any one of claims 1 to 6.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method as described in any one of claims 1 to 6.

9. A chip, characterized in that, include: A processor for retrieving and running a computer program from memory, causing a device on which the chip is mounted to perform the method as described in any one of claims 1 to 6.

10. A computer program product, characterized in that, The computer program product includes a computer storage medium storing a computer program, the computer program including instructions executable by at least one processor, which, when executed by the at least one processor, implement the method as described in any one of claims 1 to 6.

Images