Nothing Special   »   [go: up one dir, main page]

CN114547790A - Calculation method for evaluating heat insulation performance of complex multi-layer thermal protection structure - Google Patents

Calculation method for evaluating heat insulation performance of complex multi-layer thermal protection structure Download PDF

Info

Publication number
CN114547790A
CN114547790A CN202210049435.XA CN202210049435A CN114547790A CN 114547790 A CN114547790 A CN 114547790A CN 202210049435 A CN202210049435 A CN 202210049435A CN 114547790 A CN114547790 A CN 114547790A
Authority
CN
China
Prior art keywords
heat
protection structure
thermal protection
matrix
layer
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210049435.XA
Other languages
Chinese (zh)
Other versions
CN114547790B (en
Inventor
白晓辉
王玉清
张先龙
刘存良
刘冰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Beijing Power Machinery Institute
Original Assignee
Northwestern Polytechnical University
Beijing Power Machinery Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University, Beijing Power Machinery Institute filed Critical Northwestern Polytechnical University
Priority to CN202210049435.XA priority Critical patent/CN114547790B/en
Publication of CN114547790A publication Critical patent/CN114547790A/en
Application granted granted Critical
Publication of CN114547790B publication Critical patent/CN114547790B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

The invention relates to a calculation method for evaluating the heat insulation performance of a complex multilayer heat protection structure, which is characterized in that aiming at the multilayer heat protection structure, the method determines the initial parameters of the multilayer heat protection structure, including the number of layers, the thickness of each layer, the equation of each layer of thermophysical parameters such as density, heat conductivity coefficient, specific heat capacity, heat absorption heat source, air gap thickness and the change relational expression of the thermophysical parameters, establishes a differential equation set of the multilayer heat protection structure, gives a heat boundary condition and an interlayer continuous condition, divides units based on a finite unit method, disperses the differential equation set, calculates a heat capacity matrix, a heat conductivity coefficient matrix and a heat load matrix, and establishes a matrix equation set; and solving the matrix equation set by adopting a Gauss-Seidel method to obtain the node temperature and the outer wall surface temperature distribution of the multilayer thermal protection structure. The method can establish a heat transfer calculation model of the multilayer thermal protection structure, and is suitable for evaluating the heat insulation performance of the multilayer thermal protection structure in different complex thermal environments.

Description

Calculation method for evaluating heat insulation performance of complex multi-layer thermal protection structure
Technical Field
The invention belongs to the field of engine auxiliary test, and particularly relates to an infection calculation method for a multilayer thermal protection structure of a combustion chamber.
Background
With the increasing requirements on the performance of aero-engines, the temperature of fuel gas in an engine combustion chamber rises continuously, the requirement for efficient thermal protection of the metal outer wall surface of the engine is more urgent, and particularly in solid fuel ramjet engines, the thermal insulation effect of a thermal protection structure is particularly important. The engine combustion chamber thermal protection can be divided into active thermal protection and passive thermal protection, and the passive thermal protection is effective thermal protection on an engine metal shell by adopting a high-temperature-resistant and ablation-resistant thermal protection structure so as to prevent the engine metal shell from deforming or reducing the reliability under the action of high-temperature gas.
Materials with high temperature resistance and ablation resistance are often adopted in the passive thermal protection structure as an ablation layer, such as carbon/phenolic aldehyde materials, however, due to pyrolysis of the materials under the action of high temperature, pyrolysis gas is generated, a part of heat is absorbed in the process, and when the thermophysical properties of the materials are changed, the mass injection and the immersion of the pyrolysis gas into other layers can affect the thermophysical property parameters of other layers; in addition, due to the different thermal expansion coefficients between the metal structure and the heat insulation solid material, air gaps are generated between layers at a certain high temperature to play a certain heat insulation role, but the generation of the air gaps and the development of the thickness are also related to time and temperature. The above factors bring difficulty to the accurate calculation of the heat insulation performance of the multilayer heat protection structure. The existing calculation method generally adopts a finite difference method to calculate the one-dimensional transient heat conduction, and extra heat sources which change along with time can be edited in some commercial software, however, a calculation model which considers the influence of the thermal physical property parameters of the materials along with the time change or along with the temperature change, the influence of the change of an internal heat source along with the time, the influence of air gaps between layers which generate dynamic change and can quickly and accurately predict is rare. The establishment of the heat transfer calculation method of the multilayer thermal protection structure considering the complex internal thermal environment has important significance for accurately predicting the temperature change of the metal outer wall surface of the engine and designing the multilayer thermal protection structure.
Disclosure of Invention
The invention aims to avoid the defects of the prior art and provides a calculation method for evaluating the heat insulation performance of a complex multi-layer thermal protection structure, which can establish a heat transfer calculation model of the multi-layer thermal protection structure containing an internal heat source and is suitable for evaluating the heat insulation performance of the multi-layer thermal protection structure under different complex thermal environments.
In order to achieve the purpose, the invention adopts the technical scheme that: a calculation method for evaluating the thermal insulation performance of a complex multi-layer thermal protection structure, comprising the steps of:
(1) establishing a heat transfer model control differential equation set of the multilayer thermal protection structure:
establishing a one-dimensional transient multilayer heat conduction differential equation set according to the initial parameters of the multilayer heat protection structure:
layer 1
Figure BDA0003473371310000021
Layer 2
Figure BDA0003473371310000022
……
The n-th layer
Figure BDA0003473371310000023
Wherein: rhon、cn、λnThe density, specific heat and heat conductivity of the nth layer of thermal protection structure material are respectively shown, T is temperature, T is time, and x is a coordinate in the thickness direction;
the initial condition and the two boundary conditions of the multilayer thermal protection structure are expressed as:
the initial condition is that T is 0T (x,0) Tf；
The boundary condition isx＝x1: a first type or a second type or a third type of thermal boundary condition;
Figure BDA0003473371310000031
wherein, TwgIs the high temperature side wall surface temperature, h is the convective heat transfer coefficient, TwfIs the low temperature side wall surface temperature, TfIs the ambient temperature;
the boundary surface continuous condition of the multilayer thermal protection structure is as follows:
Figure BDA0003473371310000032
Figure BDA0003473371310000033
……
Figure BDA0003473371310000034
(2) adopting a finite element method to disperse a heat transfer model control differential equation set of the multilayer thermal protection structure:
dividing a multilayer thermal protection structure with the thickness of L into N one-dimensional heat conduction discrete units and N +1 nodes, wherein one discrete unit consists of 2 nodes, one discrete unit is arbitrarily selected, the nodes are i and j respectively, the length of the discrete unit is L, if the temperature on the nodes is Ti and Tj respectively, the distance node i is the temperature T (x) on the point x, and the distance node i is expressed by a shape function interpolation formula as follows:
T(x)＝NiTi+NjTj
wherein Ni and Nj are shape functions expressed as:
Figure BDA0003473371310000041
(3) obtaining a heat conduction matrix and a heat capacity matrix of each discrete unit according to the discrete units in the step (2):
the heat capacity matrix [ C ] and the heat conduction matrix [ K ] of each discrete unit are both 2 x 2 matrixes, the temperature { T } and the heat load { P } are both vector matrixes, and the matrixes are respectively expressed as follows:
Figure BDA0003473371310000042
wherein,
Figure BDA0003473371310000043
a first derivative matrix of the node temperature with respect to time;
(4) based on a variational principle, calculating the multilayer thermal protection structure area of the N one-dimensional heat conduction discrete units to obtain a transient heat conduction finite element equation as follows:
Figure BDA0003473371310000044
wherein [ C ]](N+1)×(N+1)For a matrix of heat capacities superposed by a matrix of discrete cells, [ K ]](N+1)×(N+1)A heat conduction matrix which is a superposition of matrices of discrete elements, { T }(N+1)×1For temperature arrays of superimposed nodes of each matrix of discrete elements, { P }(N+1)×1A superimposed array of temperature loads for each matrix of discrete elements,
Figure BDA0003473371310000045
a first derivative array of the superimposed node temperatures of each discrete element matrix with respect to time;
(5) and (3) solving the finite element equation of the transient heat conduction obtained in the step (4) by adopting a Gauss-Seidel method to obtain the node temperature change and the outer wall surface temperature change of the multilayer thermal protection structure, and evaluating the heat insulation performance of the multilayer thermal protection structure.
Further, among the multilayer thermal protection structure, including at least one deck pyrolysis layer that can produce pyrolysis gas under the high temperature effect, just pyrolysis layer not set up at multilayer thermal protection structure's innermost or outermost layer, pyrolysis layer corresponding one-dimensional transient state heat conduction differential equation system be:
Figure BDA0003473371310000051
in the formula:
Figure BDA0003473371310000052
the endothermic heat of the pyrolysis gas is generally calculated by the mass change rate of the pyrolysis gas and the enthalpy of pyrolysis, and the calculation formula is:
Figure BDA0003473371310000053
the interface continuous condition corresponding to the pyrolysis layer is represented as:
Figure BDA0003473371310000054
Figure BDA0003473371310000055
wherein q (t) is the heat quantity of thermal blockage formed between layers after pyrolysis gas is released, and psi is the thermal blockage coefficient.
Further, in the heat transfer model control differential equation set of the multilayer thermal protection structure, the density ρ, specific heat c and thermal conductivity λ of the multilayer thermal protection material are temperature functions, i.e. ρ (T), c (T) and λ (T), however, when the surface temperature of the multilayer thermal protection material becomes 1800K within ten seconds, the density ρ, specific heat c and thermal conductivity λ of the multilayer thermal protection material are time functions, i.e. ρ (T), c (T) and λ (T).
Further, when the metal shell of the multi-layer thermal protection structure is affected by thermal expansion, an air gap layer is added between the metal shell and the multi-layer thermal protection structure, and the thickness of the air gap layer changes with temperature, the calculation method further comprises the following steps:
(a) adopting a dynamic grid to simulate an air gap, and assuming that the coordinate of the last node of the multilayer thermal protection structure is x according to the thermophysical parameters of mixed gas in the air gap0Meanwhile, the thickness of the metal shell layer is assumed to have I units, and the length of each unit is leThe coordinate of each node is x1，x2，...，xI. When air gaps are generated among the layers, the movable grids divide the grids of the I metal layers into the air gaps, the thickness of the metal shell layer is kept unchanged, and the number of the grids of the metal layers is I-I;
(b) according to the grid number of the metal layer obtained in the step a, a one-dimensional transient heat conduction differential equation of the air gap layer is obtained as follows:
Figure BDA0003473371310000061
wherein I is the number of metal layers, I is the number of gas layers, leThe original unit length of the metal layer;
(c) the differential equation of the interlayer interface continuity condition of adding the air gap layer between the outermost layer and the adjacent layer of the outermost layer is as follows:
xn＝xn-1+δ(T)(n＞1)
in the formula: δ (T) is the thickness of the air gap.
Further, the step (4) is based on a variational principle, and the specific steps of calculating the multilayer thermal protection structure region of the N one-dimensional thermal conduction discrete units are as follows:
(41) calculating the heat capacity matrix [ C ] in the step (3) according to the shape function interpolation formula in the step (2), and concretely, the following steps are carried out:
(1) calculating heat capacity matrix
Figure BDA0003473371310000071
Figure BDA0003473371310000072
Wherein N isiAnd NjFor the shape function, there are:
Figure BDA0003473371310000073
Figure BDA0003473371310000074
thus, there are:
Figure BDA0003473371310000075
in the formula: a is the cross-sectional area of the discrete unit, and for one-dimensional calculation, the A is 1; l is the length of the unit, i.e. the distance between two nodes; rho is the density kg/m of the heat insulation layer material in which the unit is arranged3(ii) a c is the specific heat capacity J/(kg.K) of the heat insulation layer material where the unit is located;
(42) calculating the heat conduction matrix [ K ] in the step (3) according to the shape function interpolation formula in the step (2), and concretely, the following steps are carried out:
heat conduction matrix
Figure BDA0003473371310000076
Figure BDA0003473371310000077
Namely: the heat-conducting matrix is composed of two parts,
Figure BDA0003473371310000078
contributes to the heat transfer matrix for the cell interior domains;
Figure BDA0003473371310000079
contribution to the thermal conduction matrix for a third type of boundary;
then there are:
Figure BDA0003473371310000081
Figure BDA0003473371310000082
thus, there are:
Figure BDA0003473371310000083
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1; l is the length of the unit, i.e. the distance between two nodes; lambda is the thermal conductivity coefficient of the heat insulation layer material where the unit is located, W/(mK);
for the
Figure BDA0003473371310000084
On the third class of boundaries:
Figure BDA0003473371310000085
the same can be obtained:
Figure BDA0003473371310000086
thus, there are:
Figure BDA0003473371310000087
in the formula, A is a surface of a third type boundary, and for a one-dimensional unit, the calculation is that A is 1; h is the heat transfer coefficient on the third type boundary, W/(m2. K);
integrating the contributions within the discrete elements of the multilayer thermal protection structure and at the boundary of the third type yields:
Figure BDA0003473371310000088
(43) calculating the heat load { P } in the step (3) according to the shape function interpolation formula in the step (2), specifically as follows:
Pethe medicine consists of three parts:
Figure BDA0003473371310000091
in the formula:
Figure BDA0003473371310000092
is the heat load generated by the heat source in the cell,
Figure BDA0003473371310000093
the thermal load generated for the heat flow on the second type boundary,
Figure BDA0003473371310000094
heat load generated for heat exchange on the third type boundary;
wherein
Figure BDA0003473371310000095
For the node i, the following are:
Figure BDA0003473371310000096
for the node j, the following are:
Figure BDA0003473371310000097
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1; l is the length of the unit, rho is the density of the heat insulation layer material where the unit is located, Q is the heat source density W/kg of the heat insulation layer where the unit is located, heat release is positive, and heat absorption is negative.
Considering the heat load generated by the second type of boundary, if there is a density of flow q at node iiThe density of heat flow q existing at j nodejThen, there are:
Figure BDA0003473371310000098
Figure BDA0003473371310000099
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1;
considering the heat load generated by the third type boundary, if the heat exchange coefficient at the i node is hiHeat transfer coefficient h at j nodejAt an ambient temperature of TfThen, there are:
Figure BDA0003473371310000101
Figure BDA0003473371310000102
in conclusion, the following results are obtained:
Figure BDA0003473371310000103
further, the heat capacity matrix [ C ], the heat conduction matrix [ K ] and the heat load { P } are formed by superposing the N discrete unit matrices or vector matrices, and are shown as follows:
Figure BDA0003473371310000104
Figure BDA0003473371310000111
Figure BDA0003473371310000112
where ρ isi、ci、λiDensity, specific heat and thermal conductivity of each node respectively; liIs the length of the ith cell; h is the equivalent heat exchange coefficient of the outer wall surface;
Figure BDA0003473371310000113
qiand
Figure BDA0003473371310000114
internal heat source, thermal obstruction and surface heat flow density influence of each node respectively, and the unit nodes without the internal heat source, the thermal obstruction and the surface heat flow density are subjected to the internal heat source, the thermal obstruction and the surface heat flow density influence
Figure BDA0003473371310000115
qiAnd
Figure BDA0003473371310000116
are all 0.
Further, the specific steps of the step (5) are as follows:
according to the transient heat conduction finite element equation obtained in the step (4), for the first derivative array
Figure BDA0003473371310000117
The derivation can be in differential format as follows:
supposing that after m iterations, the node temperature vector { T ] at the mth step is calculatedmStep (1-1) { T (T + θ Δ T) } is performed during the iteration from step (m) to step (m +1)m}+θ{Tm+1And then, there are:
Figure BDA0003473371310000121
while letting { P (t + θ Δ t) } ═ 1- θ { P }m}+θ{Pm+1And (5) simplifying to obtain:
Figure BDA0003473371310000122
the temperature array of each node (the outer wall surface is the node n +1) of the (m +1) th step is expressed as follows:
Figure BDA0003473371310000123
in the formula:
Figure BDA0003473371310000124
then the following results are obtained:
Figure BDA0003473371310000125
wherein [ C ]]Is a heat capacity matrix, [ K ]]Is a heat conduction matrix, Δ t is a time step; { TmIs the temperature array of the mth step, { P }mThe m-th step of the thermal load array, { P }m+1The heat load array of step m +1, Tm+1，N+1The temperature of the outer wall surface in the step (m + 1).
The beneficial effects of the invention are:
(1) the invention aims to provide a heat transfer calculation method for a multilayer thermal protection structure containing a complex internal thermal environment, which can evaluate the heat insulation performance of the multilayer thermal protection structure, has a certain reference value for the design work of the multilayer thermal protection structure, and further can perform inverse calculation on the range of internal thermal parameters according to the comparison of a calculation result and an experiment result.
(2) The invention adopts a finite element method to disperse the differential equation set of the multilayer thermal protection structure, simultaneously considers the influence of complex thermal environments in the multilayer thermal protection structure, such as an internal heat source, a thermal plug, a dynamic interlayer air gap and the like, can simply and conveniently adjust geometric structure parameters, material physical property parameters, boundary conditions and internal thermal boundary conditions, adds a moving grid considering the air gap, establishes a thermal transmission calculation model of the thermal protection structure, and can more accurately predict the rule and the characteristic of the temperature change of the outer wall surface of the multilayer thermal protection structure;
(3) the invention has reliable calculation result and short calculation time, can reduce the times of the complete machine experiment on the ground of the engine, saves manpower and material resources and reduces time cost and research cost.
Drawings
FIG. 1 is a model diagram of a calculation method of a multi-layer thermal protection structure according to the present invention;
FIG. 2 is a flow chart of the calculation of the present invention using the Gaussian-Seidel method;
FIG. 3 is a one-dimensional heat conduction calculation model of a multilayer thermal protection structure of thickness L according to the present invention;
FIG. 4 is a schematic diagram of the composition of the equation set of the step m +1 of the multi-layer thermal protection structure of the present invention;
FIG. 5 is a schematic diagram of a heat transfer structure of the present invention with a five-layer structure for the multi-layer thermal protect structure;
FIG. 6 is a graph illustrating the temperature of the outer wall of the thermal protection structure according to the present invention;
fig. 7 is a three-dimensional variation of thickness versus time versus temperature within the thermal protection structure of the present invention.
Detailed Description
The principles and features of this invention are described below in conjunction with the following drawings, which are set forth by way of illustration only and are not intended to limit the scope of the invention.
Because the requirement of the performance of the aircraft engine is continuously improved, the passive thermal protection structure of the thermal protection of the combustion chamber of the engine is often made of materials with high temperature resistance and ablation resistance, in addition, because the thermal expansion coefficients of the metal structure and the thermal insulation solid materials are different, air gaps can be generated between layers at a certain high temperature to play a certain thermal insulation role, the difficulty is brought to the accurate calculation of the thermal insulation performance of the multilayer thermal protection structure, and the accurate estimation of the thermal insulation performance of the multilayer thermal protection structure becomes a technical problem to be solved, so the invention provides a calculation method for solving the problem, which comprises the following steps:
1) establishing a heat transfer model control differential equation set: the known initial parameters of the multilayer thermal protection structure comprise n layers, rho density of each layer, lambda heat conductivity coefficient, c specific heat capacity, an equation of a heat absorption heat source, air gap thickness and a change relational expression of thermophysical parameters; meanwhile, establishing a heat transfer model control differential equation set of the multilayer thermal protection structure according to the constant heat boundary condition and the interlayer interface continuous condition;
2) establishing a finite element equation set: dividing units based on a finite unit method, and dispersing the heat transfer model control differential equation set obtained in the step one, so as to calculate and obtain a heat capacity matrix, a heat conductivity coefficient matrix and a heat load matrix of the multilayer thermal protection structure; further obtaining a finite element equation set of the multilayer thermal protection structure according to the overall matrix and the vector integration;
3) solving the node temperature and the outer wall surface temperature distribution: and solving the matrix equation set by adopting a Gauss-Seidel method to obtain the node temperature and the outer wall surface temperature distribution of the multilayer thermal protection structure.
The invention adopts the basic steps to solve the technical problem that the heat insulation performance of the multilayer thermal protection structure cannot be accurately calculated, and simultaneously, the invention also provides the following specific implementation modes:
example 1: in the embodiment, the multilayer thermal protection structure comprises five layers; the method comprises the following specific steps:
1) aiming at the multilayer thermal protection structure, establishing a one-dimensional transient multilayer thermal conduction differential equation set:
layer 1
Figure BDA0003473371310000151
Layer 2
Figure BDA0003473371310000152
Layer 3
Figure BDA0003473371310000153
Layer 4
Figure BDA0003473371310000154
Layer 5
Figure BDA0003473371310000155
Wherein: rho, c and lambda are respectively density, specific heat and thermal conductivity of the material, and are generally functions of temperature, namely rho (T), c (T) and lambda (T), however, when the multilayer thermal protection material is subjected to high temperature in a short time, the functions of time, namely rho (T), c (T) and lambda (T); generally, the wall of the combustion chamber quickly becomes extremely high under the action of high-temperature combustion gas, for example, the combustion gas temperature is 2200K, and the wall may become 1800K in a matter of tens of seconds, but the material is subjected to such high temperature at a moment, the internal structure and material properties of the material change, and the change is mainly time-dependent, unlike some cases, the temperature gradually becomes high, and the material properties change as the temperature becomes high.
Figure BDA0003473371310000156
The endothermic heat of the pyrolysis gas is generally calculated by the mass change rate of the pyrolysis gas and the enthalpy of pyrolysis, i.e.
Figure BDA0003473371310000157
That is, the present method can solve the problem with pyrolized materials.
The initial conditions and the two boundary conditions of the multilayer thermal protection structure are expressed as follows:
initial conditions:
t＝0：T(x，0)＝Tf (6)
the boundary conditions adopted in this embodiment are:
x＝x1：T(x1，t)＝Twg (7)
Figure BDA0003473371310000161
wherein, TwgIs the high temperature side wall surface temperature, h is the convective heat transfer coefficient, TwfIs the low temperature side wall surface temperature, TfIs the ambient temperature.
The interfacial continuity conditions for the multilayer thermal protective structure are expressed as follows:
Figure BDA0003473371310000162
Figure BDA0003473371310000163
Figure BDA0003473371310000164
Figure BDA0003473371310000165
wherein q (t) is the heat quantity of the pyrolysis gas for forming thermal blockage between layers, psi can be taken as the thermal blockage coefficient,
Figure BDA0003473371310000166
if an air gap is added between the metal shell and the gradient thermal insulation layer in consideration of the influence of thermal expansion of the metal shell, and the thickness of the air gap varies with temperature, the calculation method is as follows:
(a) adding a differential equation of an air gap layer between the formula (4) and the formula (5), and adding a section continuous condition between the formula (12) and the formula (13);
(b) in the embodiment, the air gap is simulated by adopting the dynamic grid, and the thermophysical parameters of the mixed gas in the air gap are applied, so that calculation can be carried out. Let the coordinate of the last node of the gradient thermal insulation layer be x0Meanwhile, N units are assumed to be arranged on the thickness of the metal shell layer, and the length of each unit is leThe coordinate of each node is x1，x2，...，xI. When air gaps are generated among the layers, the movable grids divide part of grids (such as N) of the metal layers into the air gaps, the thickness of the metal shell layer is kept unchanged, and the number of the grids of the metal layers is changed into N-N. The correction of the coordinates of each node is as follows:
xn＝xn-1+δ(T)(n＞1)
Figure BDA0003473371310000171
wherein N is the number of units of the metal layer, N is the number of units of the gas layer, and leδ (T) is the thickness of the air gap for the original cell length of the metal layer.
2) Dispersing the equation set by adopting a finite element method, dividing each layer of structure into discrete elements, and establishing a heat conduction matrix and a heat capacity matrix of the elements; secondly, analyzing the heat load of the structural layer where the unit is located, such as the thermal desorption heat condition, the thermal blockage condition, the temperature of the inner boundary and the convection heat exchange condition of the outer boundary, and establishing a heat load column vector of the unit; then, the information of each unit is integrated into a system overall heat conduction matrix, a heat capacity matrix and a heat load array to form a finite element equation set. The method comprises the following specific steps:
21) discretization of the equation:
as shown in fig. 3, a heat conduction calculation model with a multilayer heat protection structure having a thickness of L is divided into N discrete units and N +1 nodes, and for a one-dimensional problem, one unit is composed of 2 nodes.
Randomly selecting 1 unit, wherein the nodes are i and j respectively, the length of the unit is l, and if the temperature on the node is Ti and Tj respectively, the temperature T (x) on the point which is at the distance of x from the node i can be interpolated by a shape function as follows:
T(x)＝NiTi+NjTj (14)
wherein Ni and Nj are shape functions expressed as:
Figure BDA0003473371310000181
the one-dimensional heat conduction problem adopts linear units, namely each unit consists of 2 nodes, the heat conduction matrix and the heat capacity matrix of each unit are both 2 x 2 matrixes, and the temperature and the heat load are both vectors, which is specifically as follows:
Figure BDA0003473371310000182
for a calculation region with n units, the finite element equation of transient heat conduction based on the variational principle is as follows:
Figure BDA0003473371310000185
wherein [ C ]](N+1)×(N+1)For a matrix of heat capacities superposed by a matrix of discrete cells, [ K ]](N+1)×(N+1)A heat conduction matrix which is a superposition of matrices of discrete elements, { T }(N+1)×1A temperature array of nodes which are superimposed on each discrete element matrix, { P }(N+1)×1A superimposed array of temperature loads for each matrix of discrete elements,
Figure BDA0003473371310000186
a first derivative array of the superimposed node temperatures of each discrete element matrix with respect to time;
22) matrix and array establishment:
with reference to fig. 1 and equation (14), the coefficient matrix and vector in the cell are calculated as follows: (221) heat capacity matrix
Figure BDA0003473371310000183
Figure BDA0003473371310000184
Wherein Ni and Nj are shape functions. For the one-dimensional problem, if the shape function of equation (14) is used, there are:
Figure BDA0003473371310000191
Figure BDA0003473371310000192
thus, there are:
Figure BDA0003473371310000193
in the formula: a is the cross-sectional area of the discrete unit, and for one-dimensional calculation, the A is 1; l is the length of the discrete unit, i.e. the distance between two nodes; rho is the density of the heat insulation layer material where the unit is located, kg/m 3; c is the specific heat capacity of the heat insulation layer material where the discrete units are located, J/(kg.K).
(222) Heat conduction matrix
Figure BDA0003473371310000194
Figure BDA0003473371310000195
Namely: the heat-conducting matrix is composed of two parts,
Figure BDA0003473371310000196
contributes to the heat transfer matrix for the cell interior domains;
Figure BDA0003473371310000197
is the contribution of the third type of boundary to the thermal conduction matrix.
Using the unit shown in fig. 3 and the shape function of (equation 14), there are:
Figure BDA0003473371310000198
Figure BDA0003473371310000199
thus, there are:
Figure BDA0003473371310000201
in the formula: a is the cross-sectional area of the discrete unit, and for one-dimensional calculation, the A is 1; l is the length of the discrete unit, i.e. the distance between two nodes; and lambda is the thermal conductivity coefficient of the heat insulation layer material where the discrete units are located, W/(mK).
For
Figure BDA0003473371310000202
On the third class of boundaries:
Figure BDA0003473371310000203
the same can be obtained:
Figure BDA0003473371310000204
thus, there are:
Figure BDA0003473371310000205
in the formula, A is a surface of a third type boundary, and for a one-dimensional unit, the calculation is that A is 1; h is the heat transfer coefficient at the third type of boundary, W/(m2. K).
Contributions within the synthesis unit and on the third class boundary are:
Figure BDA0003473371310000206
(23) heat load matrix calculation
In conjunction with the deduction in the literature, Pe consists of three parts:
Figure BDA0003473371310000207
in the formula:
Figure BDA0003473371310000208
is the heat load generated by the heat source in the cell,
Figure BDA0003473371310000209
the thermal load generated for the heat flow on the second type boundary,
Figure BDA00034733713100002010
the heat load generated for the heat exchange on the third type boundary.
Wherein
Figure BDA0003473371310000211
For node i there are:
Figure BDA0003473371310000212
for node j there is:
Figure BDA0003473371310000213
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1; l is the length of the discrete unit, rho is the density of the heat insulation layer material where the discrete unit is located, Q is the heat source density of the heat insulation layer where the discrete unit is located, W/kg, heat release is positive, and heat absorption is negative.
Considering the heat load generated by the class II boundary, if there is a heat flow density qi at node i and qj at node j, then:
Figure BDA0003473371310000214
Figure BDA0003473371310000215
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1;
if the heat transfer coefficient at the i node is hi, the heat transfer coefficient hj at the j node is hj, and the ambient temperature is Tf, the heat load generated by the class III boundary has the following components:
Figure BDA0003473371310000216
Figure BDA0003473371310000221
to sum up:
Figure BDA0003473371310000222
by calculating the heat capacity matrix, the heat conductivity coefficient matrix and the heat load array of each unit, a transient heat conduction finite element equation (36) of the whole system is formed, namely a formula (16) is obtained:
Figure BDA0003473371310000223
3) and solving the matrix equation set by adopting a Gauss-Seidel method to obtain the node temperature change and the outer wall surface temperature change of the multilayer thermal protection structure. The calculation flow chart is shown in FIG. 4.
Obtaining an array with respect to temperature:
for first derivative arrays
Figure BDA0003473371310000224
The derivation can be in differential format as follows:
assuming that the node temperature vector at the mth step { T } has been calculated after m iterationsmIterating from the mth step to the m +1 th stepIn the process of (1), let { T (T + θ Δ T) } ═ 1- θ { T }m}+θ{Tm+1And then, there are:
Figure BDA0003473371310000225
while letting { P (t + θ Δ t) } ═ 1- θ { P }m}+θ{Pm+1And substituting formula (15) to obtain:
Figure BDA0003473371310000226
the temperature array of each node (the outer wall surface is the node N +1) in the (m +1) th step is expressed as follows:
Figure BDA0003473371310000231
in the formula:
Figure BDA0003473371310000232
[C]is a heat capacity matrix, [ K ]]For a heat conduction matrix, Δ t is the time step;
Figure BDA0003473371310000233
{Tmis the temperature array of the mth step, { P }mThe m-th step of the thermal load array, { P }m+1The heat load array of the (m +1) th step, Tm+1，N+1The temperature of the outer wall surface in the step (m + 1). The calculation is changed along with the iteration of the time step, for example, the total calculation time is 500s, and the m +1 th step is the result obtained after the 500 steps are finished and the last step is finished.
Wherein, the heat capacity matrix [ C ], the heat conductivity coefficient matrix [ K ] and the temperature load array { P } are shown as follows:
Figure BDA0003473371310000234
Figure BDA0003473371310000235
Figure BDA0003473371310000241
where ρ isi、ci、λiThe density, specific heat and thermal conductivity of each node are respectively; liIs the length of the ith cell; h is the equivalent heat exchange coefficient of the outer wall surface;
Figure BDA0003473371310000242
qiand
Figure BDA0003473371310000243
internal heat source, thermal obstruction and surface heat flow density influence of each node respectively, and the unit nodes without the internal heat source, the thermal obstruction and the surface heat flow density are subjected to the internal heat source, the thermal obstruction and the surface heat flow density influence
Figure BDA0003473371310000244
qiAnd
Figure BDA0003473371310000245
are all 0.
According to the calculation process of embodiment 1, the present invention provides a specific experimental example:
the internal part of the multilayer passive thermal protection structure of a combustion chamber of an engine is eroded by high-temperature fuel gas, and the external part of the multilayer passive thermal protection structure is subjected to natural convection heat exchange. The ablation type material of the middle layer can be pyrolyzed to form an internal heat source, and high-temperature pyrolysis gas generated by the ablation type material can permeate into the inner side of the metal outer shell layer to influence the overall heat insulation effect of the metal outer shell layer. The schematic diagram of the heat transfer of the multilayer thermal protection structure is shown in FIG. 5; according to the literature, the thermal plugging coefficient generally plays a certain role in the first 30s, the thermal plugging effect on the back side is almost disappeared, and the change of the thermal plugging coefficient along with the time in the example is shown in the table 2; the thickness change of the interlayer air gap caused by the metal expansion is shown in Table 3, and the thermal physical property parameter thereof is 50And taking the value of the standard smoke of 0K. The heat protection structure in this example has five layers, the thickness and thermophysical parameters of each layer are shown in Table 1, the inner boundary condition is constant wall temperature of 1800 ℃, and the heat exchange coefficient of the outer boundary is 60 (W/m)2-K) ambient temperature of 25 ℃, internal heat source QintIs-1 x 107(W/m3) Setting the calculation time to 500s, the time step length to 1s, and selecting the Galerkin method for the time integral type. As shown in fig. 6 and 7, it can be seen that according to the method provided by the present invention, the influence of the complex thermal environment inside the multi-layer thermal protection structure, such as an internal heat source, a thermal plug, and a dynamic interlayer air gap, is considered, a thermal transfer calculation model of the thermal protection structure is accurately and conveniently established, the rule and the characteristic of the temperature change of the outer wall surface of the multi-layer thermal protection structure are predicted, and the heat insulation performance of the multi-layer thermal protection structure is realized.
TABLE 1 physical Properties of multilayer thermal protection Structure
Figure BDA0003473371310000251
TABLE 2 thermal blockage coefficient variation
Time s 0 10 20 30 After that
Coefficient of thermal blockage 0 0 0.7 0.97 0.97
TABLE 3 interlayer air gap thickness variation
Temperature K 400 500 600 700 800 900 1000
Air gap thickness (10)-6m) 4.80 9.92 15.4 20.8 27.4 34.1 41.2
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.

Claims (7)

1. A calculation method for evaluating the thermal insulation performance of a complex multi-layer thermal protection structure, characterized by comprising the following steps:
(1) establishing a heat transfer model control differential equation set of the multilayer thermal protection structure:
establishing a one-dimensional transient multilayer heat conduction differential equation set according to the initial parameters of the multilayer heat protection structure:
layer 1
Figure FDA0003473371300000011
Layer 2
Figure FDA0003473371300000012
……
The n-th layer
Figure FDA0003473371300000013
Wherein: rhon、cn、λnThe density, specific heat and heat conductivity of the nth layer of thermal protection structure material are respectively shown, T is temperature, T is time, and x is a coordinate in the thickness direction;
the initial condition and the two boundary conditions of the multilayer thermal protection structure are expressed as:
the initial condition is that T is 0T (x,0) Tf；
The boundary condition is that x is equal to x1: a first type or a second type or a third type of thermal boundary condition;
Figure FDA0003473371300000014
wherein, TwgIs the high temperature side wall surface temperature, h is the convective heat transfer coefficient, TwfIs the low temperature side wall surface temperature, TfIs the ambient temperature;
the boundary surface continuous condition of the multilayer thermal protection structure is as follows:
Figure FDA0003473371300000015
Figure FDA0003473371300000021
……
Figure FDA0003473371300000022
(2) adopting a finite element method to disperse a heat transfer model control differential equation set of the multilayer thermal protection structure:
dividing a multilayer thermal protection structure with the thickness of L into N one-dimensional heat conduction discrete units and N +1 nodes, wherein one discrete unit consists of 2 nodes, one discrete unit is arbitrarily selected, the nodes are i and j respectively, the length of the discrete unit is L, if the temperature on the nodes is Ti and Tj respectively, the distance node i is the temperature T (x) on the point x, and the distance node i is expressed by a shape function interpolation formula as follows:
T(x)＝NiTi+NjTj
wherein Ni and Nj are shape functions expressed as:
Figure FDA0003473371300000023
(3) obtaining a heat conduction matrix and a heat capacity matrix of each discrete unit according to the discrete units in the step (2):
the heat capacity matrix [ C ] and the heat conduction matrix [ K ] of each discrete unit are both 2 x 2 matrixes, the temperature { T } and the heat load { P } are both vector matrixes, and the matrixes are respectively expressed as follows:
Figure FDA0003473371300000024
wherein,
Figure FDA0003473371300000025
a first derivative matrix of the node temperature with respect to time;
(4) based on a variational principle, calculating the multilayer thermal protection structure area of the N one-dimensional heat conduction discrete units to obtain a transient heat conduction finite element equation as follows:
Figure FDA0003473371300000031
wherein [ C ]](N+1)×(N+1)For a matrix of heat capacities superposed by a matrix of discrete cells, [ K ]](N+1)×(N+1)A heat conduction matrix which is a superposition of matrices of discrete elements, { K }(N+1)×1A temperature array of nodes which are superimposed on each discrete element matrix, { P }(N+1)×1For the superimposed temperature loaded array of each matrix of discrete elements,
Figure FDA0003473371300000032
a first derivative array of the superimposed node temperatures of each discrete element matrix with respect to time;
(5) and (3) solving the finite element equation of the transient heat conduction obtained in the step (4) by adopting a Gauss-Seidel method to obtain the node temperature change and the outer wall surface temperature change of the multilayer thermal protection structure, and evaluating the heat insulation performance of the multilayer thermal protection structure.
2. The calculation method for evaluating the thermal insulation performance of a complex multilayer thermal protection structure according to claim 1, characterized in that:
among the multilayer thermal protection structure, including at least one deck can produce the pyrolysis layer of pyrolysis gas under the high temperature effect, just the pyrolysis layer not set up at multilayer thermal protection structure's innermost or outmost, the one-dimensional transient state heat conduction differential equation system that the pyrolysis layer corresponds be:
Figure FDA0003473371300000033
in the formula:
Figure FDA0003473371300000034
the endothermic heat of the pyrolysis gas is generally calculated by the mass change rate of the pyrolysis gas and the enthalpy of pyrolysis, and the calculation formula is:
Figure FDA0003473371300000041
the interface continuous condition corresponding to the pyrolysis layer is represented as:
Figure FDA0003473371300000042
Figure FDA0003473371300000043
wherein q (t) is the heat quantity of thermal blockage formed between layers after pyrolysis gas is released, and psi is the thermal blockage coefficient.
3. The calculation method for evaluating the thermal insulation performance of a complex multilayer thermal protection structure according to claim 1, characterized in that:
in the heat transfer model control differential equation set of the multilayer thermal protection structure, the density rho, the specific heat c and the heat conductivity coefficient lambda of the multilayer thermal protection material are temperature functions, namely rho (T), c (T) and lambda (T), however, when the surface temperature of the multilayer thermal protection material becomes 1800K within ten seconds, the density rho, the specific heat c and the heat conductivity coefficient lambda of the multilayer thermal protection material are time functions, namely rho (T), c (T) and lambda (T).
4. The calculation method for evaluating the thermal insulation performance of a complex multi-layered thermal protection structure according to claim 1, wherein when the metal shell of the multi-layered thermal protection structure is affected by thermal expansion, an air gap layer is added between the metal shell and the multi-layered thermal protection structure, and the thickness of the air gap layer varies with temperature, the calculation method further comprises the following steps:
(a) adopting a dynamic grid to simulate an air gap, and assuming that the coordinate of the last node of the multilayer thermal protection structure is x according to the thermophysical parameters of mixed gas in the air gap0Meanwhile, the thickness of the metal shell layer is assumed to have I units, and the length of each unit is leThe coordinate of each node is x1，x2，...，xI(ii) a When air gaps are generated among the layers, the movable grids divide the grids of the I metal layers into the air gaps, the thickness of the metal shell layer is kept unchanged, and the number of the grids of the metal layers is I-I;
(b) according to the grid number of the metal layer obtained in the step a, a one-dimensional transient heat conduction differential equation of the air gap layer is obtained as follows:
Figure FDA0003473371300000051
wherein I is the number of metal layers, I is the number of gas layers, leThe original unit length of the metal layer;
(c) the differential equation of the interlayer interface continuous condition of adding the air gap layer between the outermost layer and the adjacent layer of the outermost layer is as follows:
xn＝xn-1+δ(T) (n＞1)
in the formula: δ (T) is the thickness of the air gap.
5. The calculation method for evaluating the thermal insulation performance of a complex multi-layer thermal protection structure according to claim 1, wherein the step (4) is based on the principle of variation, and the specific steps for calculating the multi-layer thermal protection structure area of the N one-dimensional heat conduction discrete units are as follows:
(41) calculating the heat capacity matrix [ C ] in the step (3) according to the shape function interpolation formula in the step (2), and concretely, the following steps are carried out:
(1) calculating heat capacity matrix
Figure FDA0003473371300000052
Figure FDA0003473371300000053
Wherein N isiAnd NjFor the shape function, there are:
Figure FDA0003473371300000061
Figure FDA0003473371300000062
thus, there are:
Figure FDA0003473371300000063
in the formula: a is the cross-sectional area of the discrete unit, and for one-dimensional calculation, the A is 1; l is the length of the unit, i.e. the distance between two nodes; rho is the density kg/m of the heat insulation layer material in which the unit is arranged3(ii) a c is the specific heat capacity J/(kg.K) of the heat insulation layer material where the unit is located;
(42) calculating the heat conduction matrix [ K ] in the step (3) according to the shape function interpolation formula in the step (2), and concretely, the following steps are carried out:
heat conduction matrix
Figure FDA0003473371300000064
Figure FDA0003473371300000065
Namely: the heat-conducting matrix is composed of two parts,
Figure FDA0003473371300000066
contributes to the heat transfer matrix for the cell interior domains;
Figure FDA0003473371300000067
contribution to the thermal conduction matrix for a third type of boundary;
then there are:
Figure FDA0003473371300000068
Figure FDA0003473371300000071
thus, there are:
Figure FDA0003473371300000072
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1; l is the length of the unit, i.e. the distance between two nodes; lambda is the thermal conductivity coefficient of the heat insulation layer material where the unit is located, W/(mK);
for the
Figure FDA0003473371300000073
On the third class of boundaries:
Figure FDA0003473371300000074
the same can be obtained:
Figure FDA0003473371300000075
thus, there are:
Figure FDA0003473371300000076
in the formula, A is a surface of a third type boundary, and for a one-dimensional unit, the calculation is that A is 1; h is the heat transfer coefficient on the third type boundary, W/(m2. K);
integrating the contributions within the discrete elements of the multilayer thermal protection structure and at the boundary of the third type yields:
Figure FDA0003473371300000077
(43) calculating the heat load { P } in the step (3) according to the shape function interpolation formula in the step (2), specifically as follows:
Pethe medicine consists of three parts:
Figure FDA0003473371300000081
in the formula:
Figure FDA0003473371300000082
is the heat load generated by the heat source in the cell,
Figure FDA0003473371300000083
the thermal load generated for the heat flow on the second type boundary,
Figure FDA0003473371300000084
for heat exchange on a boundary of the third kindExchanging the generated heat load;
wherein
Figure FDA0003473371300000085
For the node i, the following steps are carried out:
Figure FDA0003473371300000086
for the node j, the following are:
Figure FDA0003473371300000087
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1; l is the length of the unit, rho is the density of the heat insulation layer material where the unit is located, Q is the heat source density W/kg of the heat insulation layer where the unit is located, heat release is positive, and heat absorption is negative.
Considering the heat load generated by the second type of boundary, if there is a density of flow q at node iiThe density of heat flow q existing at j nodejThen, there are:
Figure FDA0003473371300000088
Figure FDA0003473371300000089
in the formula: a is the cross-sectional area of the unit, and for one-dimensional calculation, the value of A is 1;
considering the heat load generated by the third type boundary, if the heat exchange coefficient at the i node is hiHeat transfer coefficient h at j nodejAt an ambient temperature of TfThen, there are:
Figure FDA0003473371300000091
Figure FDA0003473371300000092
in conclusion, the following results are obtained:
Figure FDA0003473371300000093
6. the calculation method for evaluating the thermal insulation properties of a complex multilayer thermal protection structure according to claim 5, characterized in that: the heat capacity matrix [ C ], the heat conduction matrix [ K ] and the heat load { P } are formed by superposing the N discrete unit matrixes or vector matrixes, and are shown as follows:
Figure FDA0003473371300000094
Figure FDA0003473371300000095
Figure FDA0003473371300000101
wherein ρi、ci、λiThe density, specific heat and thermal conductivity of each node are respectively; l. theiIs the length of the ith cell; h is the equivalent heat exchange coefficient of the outer wall surface;
Figure FDA0003473371300000102
qiand
Figure FDA0003473371300000103
internal heat source, thermal chokes and surface heat flow density effects for each node respectivelyFor a unit node without internal heat source, no thermal blockage and no surface heat flux density, it
Figure FDA0003473371300000104
qiAnd
Figure FDA0003473371300000105
are all 0.
7. Calculation method for evaluating the thermal insulation properties of a complex multilayer thermal protection structure according to any one of claims 1 to 6, characterized in that: the specific steps of the step (5) are as follows:
according to the transient heat conduction finite element equation obtained in the step (4), for the first derivative array
Figure FDA0003473371300000106
The following can be derived using the differential format:
supposing that after m iterations, the node temperature vector { T ] at the mth step is calculatedmStep (1-1) { T (T + θ Δ T) } is performed during the iteration from step (m) to step (m +1)m}+θ{Tm+1And then, there are:
Figure FDA0003473371300000107
while making { P (t + θ Δ t) } (1- θ) { P }m}+θ{Pm+1And (5) simplifying to obtain:
Figure FDA0003473371300000108
then, at each node of the (m +1) th step, the outer wall surface is a node n +1, and the temperature array is expressed as follows:
Figure FDA0003473371300000111
in the formula:
Figure FDA0003473371300000112
then the following results are obtained:
Figure FDA0003473371300000113
wherein [ C ]]Is a heat capacity matrix, [ K ]]Is a heat conduction matrix, Δ t is a time step; { TmIs the temperature array of the mth step, { P }mThe m-th step of the thermal load array, { P }m+1The heat load array of the (m +1) th step, Tm+1，N+1The temperature of the outer wall surface in the step (m + 1).
CN202210049435.XA 2022-01-17 2022-01-17 Calculation method for evaluating heat insulation performance of complex multilayer heat protection structure Active CN114547790B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210049435.XA CN114547790B (en) 2022-01-17 2022-01-17 Calculation method for evaluating heat insulation performance of complex multilayer heat protection structure

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210049435.XA CN114547790B (en) 2022-01-17 2022-01-17 Calculation method for evaluating heat insulation performance of complex multilayer heat protection structure

Publications (2)

Publication Number Publication Date
CN114547790A true CN114547790A (en) 2022-05-27
CN114547790B CN114547790B (en) 2024-02-23

Family

ID=81671311

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210049435.XA Active CN114547790B (en) 2022-01-17 2022-01-17 Calculation method for evaluating heat insulation performance of complex multilayer heat protection structure

Country Status (1)

Country Link
CN (1) CN114547790B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115408758A (en) * 2022-09-19 2022-11-29 西安交通大学 Node division and heat conduction calculation method for distributed multilayer material thermal component in containment
CN115544818A (en) * 2022-12-05 2022-12-30 中国空气动力研究与发展中心低速空气动力研究所 Grid division and heat conduction calculation method for multilayer heterogeneous thin-wall solid heat conduction calculation
CN115577566A (en) * 2022-11-15 2023-01-06 中国空气动力研究与发展中心计算空气动力研究所 Processing method, device, equipment and medium for continuous ablation of multilayer heat-proof structure
CN118551483A (en) * 2024-07-22 2024-08-27 中国人民解放军国防科技大学 Annular heat-resistant layer variable thickness design method, device, equipment and medium

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101078698A (en) * 2007-06-26 2007-11-28 东南大学 Positive type detection method for protecting integral heat-insulation property of structure
KR20150124724A (en) * 2014-04-29 2015-11-06 현대자동차주식회사 Device and method for detecting thermal behavior of engine insulation coating layer
CN113071718A (en) * 2021-02-26 2021-07-06 北京空间飞行器总体设计部 Thermal protection device for lunar surface takeoff riser and heat insulation performance calculation method thereof

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101078698A (en) * 2007-06-26 2007-11-28 东南大学 Positive type detection method for protecting integral heat-insulation property of structure
KR20150124724A (en) * 2014-04-29 2015-11-06 현대자동차주식회사 Device and method for detecting thermal behavior of engine insulation coating layer
CN113071718A (en) * 2021-02-26 2021-07-06 北京空间飞行器总体设计部 Thermal protection device for lunar surface takeoff riser and heat insulation performance calculation method thereof

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张津, 吴护林: "发动机喷管隔热涂层的设计和模拟计算", 兵工学报, no. 02, 20 June 2002 (2002-06-20) *
马秀峰;李志永;郑日恒;李立翰;陈静敏;陈玉峰;石兴;: "冲压发动机用梯度隔热层隔热性能研究", 推进技术, no. 03, 2 March 2015 (2015-03-02) *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115408758A (en) * 2022-09-19 2022-11-29 西安交通大学 Node division and heat conduction calculation method for distributed multilayer material thermal component in containment
CN115408758B (en) * 2022-09-19 2023-05-09 西安交通大学 Distributed multi-layer material thermal component node division and heat conduction calculation method in containment
CN115577566A (en) * 2022-11-15 2023-01-06 中国空气动力研究与发展中心计算空气动力研究所 Processing method, device, equipment and medium for continuous ablation of multilayer heat-proof structure
CN115577566B (en) * 2022-11-15 2023-03-10 中国空气动力研究与发展中心计算空气动力研究所 Processing method, device, equipment and medium for continuous ablation of multilayer heat-proof structure
CN115544818A (en) * 2022-12-05 2022-12-30 中国空气动力研究与发展中心低速空气动力研究所 Grid division and heat conduction calculation method for multilayer heterogeneous thin-wall solid heat conduction calculation
CN115544818B (en) * 2022-12-05 2023-02-28 中国空气动力研究与发展中心低速空气动力研究所 Grid division and heat conduction calculation method for multilayer heterogeneous thin-wall solid heat conduction calculation
CN118551483A (en) * 2024-07-22 2024-08-27 中国人民解放军国防科技大学 Annular heat-resistant layer variable thickness design method, device, equipment and medium

Also Published As

Publication number Publication date
CN114547790B (en) 2024-02-23

Similar Documents

Publication Publication Date Title
CN114547790A (en) Calculation method for evaluating heat insulation performance of complex multi-layer thermal protection structure
Daryabeigi Thermal analysis and design optimization of multilayer insulation for reentry aerodynamic heating
Daryabeigi Thermal analysis and design of multi-layer insulation for re-entry aerodynamic heating
Peng et al. Numerical simulation of natural convection in a concentric annulus between a square outer cylinder and a circular inner cylinder using the Taylor-series-expansion and least-squares-based lattice Boltzmann method
CN113326564A (en) Method for obtaining transient temperature field of gradient composite heat insulation structure
CN109960878A (en) A kind of active/passive thermal protection system coupling design method towards hypersonic aircraft totality
CN111460709A (en) Method for predicting temperature distribution and deformation of part in fused deposition manufacturing process
Zhang et al. Preliminary study on the thermal insulation of a multilayer passive thermal protection system with carbon-phenolic composites in a combustion chamber
CN115809515B (en) Multi-layer heat insulation structure optimization design method for high-speed aircraft
Li et al. Design and modeling of a multiscale porous ceramic heat exchanger for high temperature applications with ultrahigh power density
Mohammadiun Time-dependent heat flux estimation in multi-layer systems by inverse method
Kang et al. Thermomechanical characterization of hot surface ignition device using phenomenological heat flux model
Naved et al. Numerical simulation of transpiration cooling for a high-speed vehicle with substructure
CN113139324B (en) Finite element method for calculating effective heat conductivity of all-ceramic micro-encapsulated fuel pellet
Cross Coupled simulations of finite-rate ablation with pyrolysis in rocket nozzles
Aljaghtham A comparative performance analysis of thermoelectric generators with a novel leg geometries
CN112632772A (en) Extreme environment-oriented multifunctional collaborative design method for thermal protection material
Zhang et al. Investigation on a CMC-aided multilayer thermal structure and the functionality of active cooling
Johnson Optimization of layer densities for multilayered insulation systems
CN109598059A (en) A kind of thermal protection system optimum design method and designing system based on agent model
CN115408758B (en) Distributed multi-layer material thermal component node division and heat conduction calculation method in containment
Yang et al. A fin analogy model for thermal analysis of regeneratively cooled rocket engines
Xiao et al. Insight into chemical reaction kinetics effects on thermal ablation of charring material
Yanxia et al. Heat transfer performance study of composite phase change materials for thermal management
Singhal et al. Convective and radiative thermal analysis of composite wall with nonlinear temperature-dependent properties

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant