CN103366065A - Size optimization design method for aircraft thermal protection system based on section reliability - Google Patents

Abstract

A size optimization design method for an aircraft thermal protection system based on section reliability comprises the following steps: 1. mathematical modeling, according to the structural form of the thermal protection system, determining a size related design variable, and according to the temperature requirement, establishing a heat proof structure weight reduction optimization model; 2. describing the material, external load and uncertainty in an initial/boundary value condition by utilizing the section in a quantifying manner; 3. carrying out reliability switch on the temperature constraint condition containing an uncertain parameter based on a section possibility degree; 4. utilizing the section limited volume method, realizing quickly solving of the structure transient temperature response range, and changing a two-layer nesting optimization problem into a conventional single-layer optimization problem; 5. carrying out programming calculation on switched single-layer optimization problem containing the reliability restraint, and determining an optimal solution to reaching an optimal weight reduction effect. The invention can systematically solve the thermal protection system structure optimization problem containing multisource section uncertainty, improves the computational efficiency, and ensures the reliability and the safety of the structure use.

Description

A kind of size optimization design method for aircraft thermal protection system based on section reliability

Technical field

The present invention relates to aircraft thermal protection system optimum structure design method field, more particularly to a kind of thermal protection system sizing method optimized based on interval heat conduction analysis and section reliability.

Background technology

With the fast development of Chinese Space technology, spacecraft is entered the orbit or returned ground, to be got through the earth's atmosphere, and serious Aerodynamic Heating effect is especially will suffer from return phase.The various pieces of spacecraft structure have different requirements to temperature, and most of instrument and equipment requirement is within the scope of certain temperature, it is therefore necessary to using thermal protection system with to spacecraft progress temperature control.It is, in general, that thermal protection system thickness is bigger, thermal protection ability is better；But on the other hand, the architecture quality increase that thickness increase is brought greatly reduces the overall performance of aircraft, therefore accurate temperature field analysis is carried out to thermal protection system, just had great importance to obtain suitable thermal protection system thickness by effective dimensionally-optimised design.

It is well-known, uncertainty is widely present in objective world, aircraft thermal protection structure system its thermal force, physical dimension, material property etc. during the manufacturing and use are inevitably influenceed by uncertain factors such as production defect, measurement errors, these all can produce influence to the calculating of thermal protection system heat transfer process, cause structure temperature field distribution to produce fluctuation, or even the possibility failed occur.Traditional thermal protection system thermodynamic analysis and optimization design are all based on deterministic models implementation, it is impossible to embody practical problem and contain probabilistic objective essence, usually these designs can cause the increase of architecture quality and certain unsafe factor.

To reduce the various uncertain influences to thermal protection system heat-proof quality as far as possible, designer should just predict the change that may occur in the design phase, and take effective computational methods and design, the insensitivity that enhancing temperature field is fluctuated to Parameters variation, so as to improve aircraft body structure and the safety in utilization of instrument and equipment, here it is carrying out the original intention of size optimization design method for aircraft thermal protection system research based on reliability theory.For practical structures many parameter uncertainty quantification the problem of, to obtain enough sample informations, often seem that extremely difficult or cost is too high to construct its accurate probability-distribution function or fuzzy membership function., only need to be by the bound of less information acquisition variable and interval model is a class relatively new uncertain quantitative method, therefore embody in terms of modeling more preferable economy.In addition, being required for becoming more meticulous for structural reliability problem research so that traditional design concept based on method of safety coefficients seems overly conservative.Therefore, various probabilistic influences are just taken into full account from the establishment stage of calculating and Optimized model, propose the size optimization design method for aircraft thermal protection system based on section reliability, deficiency for making up existing thermal protection system transient state temperature field numerical computations and optimization design, with important engineering application value.

The content of the invention

The technology of the present invention solves problem：Overcome existing design technology not enough present in aircraft thermal protection system structure optimization, a kind of size optimization design method for aircraft thermal protection system based on section reliability is provided, in the numerical computations that interval Finite Volume Method is incorporated into the thermal protection system transient state temperature field containing interval uncertain parameter, on the premise of ensureing that temperature field meets aircraft body structure use requirement, a kind of physical dimension Robust-Design scheme of reduction system oeverall quality has been obtained.

The technology of the present invention solution：A kind of size optimization design method for aircraft thermal protection system based on section reliability, comprises the following steps：

Step one：It is determined that needing the design parameter for optimizing the Basic Design variable of the aircraft thermal protection system of design and correlation, design variable x=(x₁,x₂,x₃,x₄)^TTo represent, wherein：

x₁、x₂、x₃、x₄Radiation coating, upper surface nonwoven fabric layer, heat insulation layer, the thickness of lower surface nonwoven fabric layer are represented respectively；According to actual requirement of engineering, the initial range of above design variable is determined；

Design parameter includes the physical attribute of three kinds of materials, such as density of material ρ_i, thermal conductivity factor k_i, specific heat c_iI=1,2,3；The heat load born, is represented with heat flow density q；For convenience, the relevant parameter in this calculating and Optimized model is expressed as to vectorial α form, i.e.,：

α=(ρ₁,ρ₂,ρ₃,k₁,k₂,k₃,c₁,c₂,c₃,q)^T；

Step 2：Using thermal protection system oeverall quality as the object function of optimization, the temperature in use scope of aircraft body structure sets up following Optimized model as constraints：

$\underset{x}{\min }M(\mathrm{\α},x)$

s.t.T_j(α,x)≤T^maxJ=1,2 ..., m

$\underset{\‾}{x}\≤x\≤\stackrel{\‾}{x}$

(α x) represents structure gross mass to wherein M；T^maxThe temperature upper limit that can bear for main structure；M is the number of constraints；

It is the bound of the design variable initial range defined in step one；

Step 3：The heat transfer problem of thermal protection system is a transient, the hot attribute of heat-barrier material can be changed with temperature, take into full account the fluctuation that temperature change is brought to material properties, and the uncertain factor such as the measurement error of hot-fluid, each uncertain parameter in this thermal protection system is described with interval vector, i.e.,：

$\mathrm{\α}\∈{\mathrm{\α}}^I=[\underset{\‾}{\mathrm{\α}},\stackrel{\‾}{\mathrm{\α}}]=[{\mathrm{\α}}^c-\mathrm{\Δ\α},{\mathrm{\α}}^c+\mathrm{\Δ\α}]={\mathrm{\α}}^c+[-\mathrm{\Δ\α},\mathrm{\Δ\α}]={\mathrm{\α}}^c+\mathrm{\Δ\α\δ}$

Wherein

The upper bound and the lower bound of vector, α are represented respectively^c, Δ α is respectively the nominal value and radius of vector, and is met：

${\mathrm{\α}}^c={({{\mathrm{\α}}_1}^c,...,{{\mathrm{\α}}_n}^c)}^T={(({\stackrel{\‾}{\mathrm{\α}}}_1+{\underset{\‾}{\mathrm{\α}}}_1)/2,...,({\stackrel{\‾}{\mathrm{\α}}}_n+{\underset{\‾}{\mathrm{\α}}}_n)/2)}^T$

$\mathrm{\Δ\α}={(\mathrm{\Δ}{\mathrm{\α}}_1,...,\mathrm{\Δ}{\mathrm{\α}}_n)}^T={(({\stackrel{\‾}{\mathrm{\α}}}_1-{\underset{\‾}{\mathrm{\α}}}_1)/2,...,({\stackrel{\‾}{\mathrm{\α}}}_n-{\underset{\‾}{\mathrm{\α}}}_n)/2)}^T;$

δ=[- 1,1]

Step 4：Constraints conversion based on section reliability, when the design parameter vector α in the dimensionally-optimised model of thermal protection system changes in its interval range, the temperature-responsive T in step 2 constraints_j(α, x) is no longer traditional fixed function, and can be by interval function T_j(α^I, x) substituted；, might as well be by the reliability index η of temperature restraint according to the definition of reliability, it is necessary to ensure the certain requirement of the possibility sexual satisfaction that constraints is set up to improve the safety and stability that structure is used_jTo represent, then the constraints based on section reliability is then converted to：

Poss(T_j(α^I,x)≤T^max)≥η_j

Wherein Poss represents the probability that inequality is set up.Assuming that interval variable is equally distributed in the small perturbation range of its nominal value, then this probability can be calculated by general interval possibility degree calculation formula：

$\mathrm{Poss}(T_j({\mathrm{\α}}^I,x)\≤T^{\max })=\left\{\begin{array}{cc}1& {\stackrel{\‾}{T}}_j({\mathrm{\α}}^I,x)\≤T^{\max }\\ \frac{T^{\max }-{\underset{\‾}{T}}_j({\mathrm{\α}}^I,x)}{{\stackrel{\‾}{T}}_j({\mathrm{\α}}^I,x)-{\underset{\‾}{T}}_j({\mathrm{\α}}^I,x)}& {\underset{\‾}{T}}_j({\mathrm{\α}}^I,x)\≤T^{\max }\≤{\stackrel{\‾}{T}}_j({\mathrm{\α}}^I,x)\\ 0& {\underset{\‾}{T}}_j({\mathrm{\α}}^I,x)\≥T^{\max }\end{array}\right.$

Wherein

T _j(α^I, x) it is respectively interval function T_j(α^I, upper bound x) and lower bound, i.e.,：

${\stackrel{\‾}{T}}_j({\mathrm{\α}}^I,x)=\underset{\mathrm{\α}\∈{\mathrm{\α}}^I}{\max }{T}_j(\mathrm{\α},x)$ ${\underset{\‾}{T}}_j({\mathrm{\α}}^I,x)=\underset{\mathrm{\α}\∈{\mathrm{\α}}^I}{\max }{T}_j(\mathrm{\α},x)$

Pass through the processing of this step so that the Optimized model set up in step 2 is converted to the nested Optimized model of the complexity containing reliability index, including inside and outside bilevel optimization, wherein outer layer optimizes for design vector x=(x₁,x₂,x₃,x₄)^TOptimizing, can be realized by step 6；And internal layer optimization is then used to calculate temperature field T_j(α^I, x) on block design parameter alpha^IThe bound of response range, can pass through step 5 Equivalent realization；

Step 5：The rapid solving of temperature field response range, in addition to optimization method, can also try to achieve temperature field on block design parameter alpha by interval numerical computations^IResponse range.For the solution of the algebraic equation containing interval parameter, traditional perturbation method can approximately try to achieve response range, but due to only remaining linear term during matrix inversion, therefore often bring than larger deviation.The present invention is by means of improved Neumann expansion techniques, establish the interval Finite Volume Method suitable for the equation of heat conduction, thermal protection system temperature field response range can fast and accurately be determined, optimize so as to instead of the internal layer described in step 4 in nested Optimized model with interval arithmetic, become two layers of nested optimization problem into conventional individual layer optimization problem, substantially increase optimization computational efficiency.Specific implementation method is as follows：

The limited configurations volume discrete model of thermal protection system is initially set up, using six point symmetry forms of second order accuracy, the following transient state temperature field limited bulk algebraic equation containing interval parameter can be obtained：

A(α^I)T^k+1=B (α^I)T^k+F(α^I)

Wherein T^kRepresent temperature vector at all nodes of kth time step；

The coefficient matrix and right-hand-side vector of above-mentioned equation are carried out into Taylor expansion at parameter intermediate value to obtain：

$A\left({\mathrm{\α}}^I\right)=A\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂A}{\∂{\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}({\mathrm{\α}}_i-{\mathrm{\α}}_i^c)=A\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂A}{{\∂\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}\mathrm{\Δ}{\mathrm{\α}}_i{\mathrm{\δ}}_i=A^c+\mathrm{\Δ}{A}^I$

$B\left({\mathrm{\α}}^I\right)=B\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂B}{\∂{\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}({\mathrm{\α}}_i-{\mathrm{\α}}_i^c)=B\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂B}{{\∂\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}\mathrm{\Δ}{\mathrm{\α}}_i{\mathrm{\δ}}_i=B^c+\mathrm{\Δ}{B}^I$

$F\left({\mathrm{\α}}^I\right)=F\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂F}{\∂{\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}({\mathrm{\α}}_i-{\mathrm{\α}}_i^c)=F\left({\mathrm{\α}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂F}{{\∂\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}\mathrm{\Δ}{\mathrm{\α}}_i{\mathrm{\δ}}_i=F^c+\mathrm{\Δ}{F}^I$

Further understood using Neumann series：

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\∞}{\mathrm{\Σ}}}{(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂A}{{\∂\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}\mathrm{\Δ}{\mathrm{\α}}_i^I{\left(A^c\right)}^{-1})}^r$

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\∞}{\mathrm{\Σ}}}{(-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_i^I{A}_i)}^r$

Wherein

When the symbolic variable r in sum formula takes different value, in above formula

Item can specifically deploy, and can be obtained after merging similar terms：

Ensureing norm | | Δ α_iA_i| | on the premise of ＜ 1 is set up, then level is several

It is convergent；Therefore, cast out the cross term in above formula, can obtain：

${(A^c+{\mathrm{\ΔA}}^I)}^{-1}\≈{\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\underset{r=1}{\overset{\∞}{\mathrm{\Σ}}}{\left({-\mathrm{\Δ\α}}_i^I{A}_i\right)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^I$

Wherein $E_i^I=\frac{{-\mathrm{\Δ\α}}_i^I{A}_i}{I+{\mathrm{\Δ\α}}_i^I{A}_i}$

It is updated to formula：

(A^c+ΔA^I)((T^k+1)^c+Δ(T^k+1)^I)=(B^c+ΔB^I)((T^k)^c+Δ(T^k)^I)+(F^c+ΔF^I)

In, using the basic operation rule of intervl mathematics, it can obtain：

${\left(T^{k+1}\right)}^c={\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^c)[B^c{\left(T^k\right)}^c+F^c]$

$\mathrm{\Δ}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{E}_i^I(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^c)(B^c\mathrm{\Δ}{\left(T^k\right)}^I+{\mathrm{\ΔB}}^I{\left(T^k\right)}^c+{\mathrm{\ΔF}}^I)]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{E}_i\·\mathrm{\δ}\·(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^c)(B^c\mathrm{\Δ}\left(T^k\right)\·\mathrm{\δ}+\mathrm{\ΔB}\·\mathrm{\δ}{\·\left(T^k\right)}^c+\mathrm{\ΔF}\·\mathrm{\δ})]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{\Δ}{E}_i(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^c)(B^c\mathrm{\Δ}\left(T^k\right)+\mathrm{\ΔB}{\left(T^k\right)}^c+\mathrm{\ΔF})]\·\mathrm{\δ}$

$=\mathrm{\Δ}\left(T^{k+1}\right)\·\mathrm{\δ}=\mathrm{\Δ}\left(T^{k+1}\right)\·[-\mathrm{1,1}]$

Wherein：

$\mathrm{\Δ}\left(T^{k+1}\right)=|\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left(A^c\right)}^{-1}{\mathrm{\ΔE}}_i(B^c{\left(T^k\right)}^c+F^c)+{\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{E}_i^c)(B^c{\mathrm{\ΔT}}^k+\mathrm{\ΔB}{\left(T^k\right)}^c+\mathrm{\ΔF})|$

Then the interval bound of thermal protection structure transient temperature response is：

${\stackrel{\‾}{T}}^{k+1}={\left(T^{k+1}\right)}^c+\mathrm{\Δ}\left(T^{k+1}\right)$ T ^k+1=(T^k+1)^c-Δ(T^k+1)

In this step, conventional finite volume method is combined with interval mathematical theory, establishes a kind of interval Finite Volume Method calculated suitable for transient heat conduction containing interval parameter；By retaining the part higher order term in Neumann series, computational accuracy is greatly improved in the case where amount of calculation allows, this just provides a strong instrument for the approximate processing of internal layer optimization problem in nesting Optimized model；

Step 6：The solution of deterministic optimization problem, interval limited bulk side in step 5 is converted into individual layer deterministic optimization problem to the rapid solving of structure transient state temperature field, the former nested optimization problem containing interval parameter；Using simulated annealing, calculation procedure is write, maximum cycle Iter is defined_maxWith convergence factor ε, when any one in following three condition is met, calculates and terminate：

（1）Loop iteration frequency n ＞ Iter_max；

（2）In double iterative process, object function nominal value relative variation is met：

$\left|\frac{M^{(i+1)}\left({\mathrm{\&alpha;}}^c\right)-M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}{M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}\right|<\mathrm{\&epsiv;};$

（3）   ||x⁽ⁱ⁺¹⁾-x⁽ⁱ⁾||₂＜ ε

Wherein | | | |₂Represent 2 norms of vector；

When reaching condition（1）When, the new initial value of design variable is given, and be brought into algorithm and recalculate；When algorithm is because of condition（2）Or（3）During termination, the result of calculation x of ith iteration process is taken⁽ⁱ⁾As the optimal value of design variable, the dimensionally-optimised design process of thermal protection system based on section reliability is completed, with the weight loss effect being optimal.

The advantage of the present invention compared with prior art is：

（1）Compared with traditional structural optimization problems, the Optimized model set up takes into full account that the multi-source caused by the complicated thermal force environment of aircraft thermal protection system is uncertain, so as to improve thermal protection system safety in utilization and stability, result of calculation has prior directive significance to its structure design.

（2）Uncertain factor in calculating and optimization problem is characterized with interval model, the rigors of traditional stochastic model and fuzzy model to amounts of specimen information are greatly reduced.

（3）By retaining the higher order term in Neumann series, the interval Finite Volume Method suitable for heat conduction analysis is established, the response range of each moment structure temperature can be quickly determined, so as to instead of the internal layer optimization in nested Optimized model with interval arithmetic.And compared with traditional Novel Interval Methods, computational accuracy is significantly improved.

Brief description of the drawings

Fig. 1 thermal protection systems containing interval parameter reliability Optimum Design flow；

Fig. 2 aircraft thermal protection system structural representations；

Fig. 3 constraints reliable realization principle schematics.

Embodiment

The present invention will be further described with reference to the accompanying drawings and examples.

In order to which the present invention is discussed in detail, definition and its algorithm of the interval analysis operation used in the present invention are introduced first.Real number field is let R be, for two given real numbersAnd

Then：

$x^I=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]=\{x\&Element;R|\underset{\&OverBar;}{x}\&le;x\&le;\stackrel{\&OverBar;}{x}\}$

Referred to as bounded closed interval, is also interval number, referred to as interval.Wherein claimxFor interval lower bound or lower extreme point, claim

For the interval upper bound or upper extreme point.If two intervals

WithThe corresponding difference of end points up and down is equal, then claims this two intervals equal, evenx=yAnd

Then x^I=y^I.Claim in addition

With

Respectively interval x^INominal value and radius.

For arbitrary two intervals in real number field

Its interval arithmetic is defined as：

$x^I+y^I=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]+[\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}]=[\underset{\&OverBar;}{x}+\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{x}+\stackrel{\&OverBar;}{y}]$

$x^I-y^I=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]-[\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}]=[\underset{\&OverBar;}{x}-\stackrel{\&OverBar;}{y},\stackrel{\&OverBar;}{x}-\underset{\&OverBar;}{y}]$

$x^I\&CenterDot;y^I=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]\&CenterDot;[\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}]=\left[\min \right\{\underset{\&OverBar;}{x}\underset{\&OverBar;}{y},\underset{\&OverBar;}{x}\stackrel{\&OverBar;}{y},\stackrel{\&OverBar;}{x}\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{\mathrm{xy}}\},\max \{\underset{\&OverBar;}{x}\underset{\&OverBar;}{y},\underset{\&OverBar;}{x}\stackrel{\&OverBar;}{y},\stackrel{\&OverBar;}{x}\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{\mathrm{xy}}\left\}\right]$

$x^I/y^I=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]/[\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}]=[\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}]\&CenterDot;[1/\stackrel{\&OverBar;}{y},1/\underset{\&OverBar;}{y}]0\&NotElement;y^I$

Size optimization design method for aircraft thermal protection system described in detail below based on section reliability：

The present invention is applied to the dimensionally-optimised problem of thermal protection system containing interval uncertain parameter.Present embodiment illustrates proposed section reliability optimization method by taking the ceramic blanket thermal protection system scheme " AFRSI " that NASA AMES research centers are proposed as an example.In addition, the section reliability optimization method of this aircraft thermal protection system can be generalized in the uncertain optimization design of other labyrinths containing interval parameter.

The structural model of this thermal protection system is as shown in Figure 2, it is considered to which size is 0.3m × 0.3m daughter board.The structure is followed successively by from outer surface to main structure：C9 protects radiation coating, AB312 upper surfaces nonwoven fabric layer, Q-fiber Felt insulation material layers, AB312 lower surface nonwoven fabric layers, using hot-fluid as outside thermal force.Structure carries out discretization using hexahedral element, extracts lower surface nonwoven fabric layer and actual work temperature of the temperature maximum at all common nodes of main structure as main structure.

The section reliability optimization process of this thermal protection system is as shown in Figure 1, on the basis of traditional Optimized model, take into full account system in itself and external applied load uncertainty, quantitative description is carried out to uncertain parameter using interval, it is introduced into reliability index and the transformation model of constraints in optimization problem is set up based on interval possibility degree, while quickly tries to achieve the bound of structure transient state temperature field response using the interval Finite Volume Method of proposition.Using simulated annealing, calculation procedure is write, can be required to choose optimal thermal protection system design size according to designer.Following several steps can be divided into carry out：

Step one：It is determined that needing the design parameter for optimizing the Basic Design variable of the thermal protection system of design and correlation, design variable x=(x₁,x₂,x₃,x₄)^TTo represent, wherein：

x₁、x₂、x₃、x₄：Radiation coating, upper surface nonwoven fabric layer, heat insulation layer, the thickness of lower surface nonwoven fabric layer are represented respectively；

During initial designs, the thickness of design variable is set as x₁=8mm, x₂=5mm, x₃=25mm, x₄=5mm.In order to meet actual requirement of engineering, above design variable has the size requirement of itself, i.e.,：

10mm≤x₃≤40mm

2mm≤x_i≤ 15mm i=1,2,4

In this thermal protection system structural model, the density p of C9 coating materials₁=2.0 × 10³kg/m³, thermal conductivity factor k₁=0.64W/ (m DEG C), specific heat c₁=628J/ (kg DEG C)；AB312 weaves cotton cloth the density p of layer material₂=985.15kg/m³, thermal conductivity factor k₂=7.1 × 10^-2W/ (m DEG C), specific heat c₂=558J/ (kg DEG C)；The density p of Q-fiber Felt heat-insulating materials₃=56.1kg/m³, thermal conductivity factor k₃=8.2 × 10^-3W/ (m DEG C), specific heat c₃=605J/ (kg DEG C)；The function that heat flow density is changed over time is q (t)=[12000- (t-800)²/60]W/m²；

For convenience, the relevant parameter in this calculating and Optimized model is expressed as to vectorial α form, i.e.,：

α=(ρ₁,ρ₂,ρ₃,k₁,k₂,k₃,c₁,c₂,c₃,q)^T；

Step 2：During transient heat conduction, the characteristic time is extracted（Respectively 200s, 500s, 800s, 1200s, 1600s）Lower surface nonwoven fabric layer, as the actual work temperature of main structure, uses T respectively with the temperature maximum at all common nodes of main structure_j(α, x) j=1,2 ..., 5 represent.The maximum temperature that constraints is taken as main structure in the working time is no more than 150 DEG C, i.e. T^max=150 DEG C.Under this temperature constraint, with the gross mass M of thermal protection system structure, (α is x) design object, can set up such as next Optimized model：

$\underset{x}{\min }M(\mathrm{\&alpha;},x)$

s.t.T_j(α,x)≤T^maxJ=1,2 ..., 5；

10mm≤x₃≤40mm

2mm≤x_i≤ 15mm i=1,2,4

Step 3：Because during Transient Heat Transfer, the hot attribute of structural material can be changed with temperature, and temperature change brings certain fluctuation to material properties, while there is certain error in the measurement of hot-fluid.Consider less on probabilistic information content in Practical Project problem, describe each uncertain parameter of this thermal protection system in the present invention using interval.It there is no harm in the perturbation that each design parameter listed in setting procedure one has 5% near its nominal value, i.e.,：

α∈α^I=α^c* [0.95,1.05]=α^c+Δα·[-1,1]

The α of wherein Δ α=0.05^c；

Step 4：Constraints conversion based on section reliability, on the premise of can tolerating constraints destruction to a certain degree in view of designer, the reliability index provided for designer, utilize interval possibility degree calculation formula, the transformation model of constraints is set up, under the conditions of meter and various parameters fluctuating change so that design point is still in feasible zone, the requirement of reliability is met, as shown in Figure 3.Transverse and longitudinal coordinate x₁And x₂Two design variables in two-dimensional space are represented respectively, and the critical condition of constraints uniformly uses g₁=0 and g (x)₂=0 (x) represent.A and B represent traditional optimal solution and the optimal solution obtained based on reliability optimization respectively in figure, and solid line and dotted line are to represent tradition optimization and the feasible zone border corresponding to reliability optimization respectively.It can be seen that, traditional optimal solution A be often positioned in feasible zone border or its near, but it is due to the influence of uncertain factor, constraints can change, one of which situation is exactly that feasible zone border changes to dotted line from solid line, so traditional optimal solution A is located at outside new feasible zone, does not meet design requirement；And the optimal solution B obtained with reliability optimization then still meets the requirement of new constraints.In the present embodiment, in order to reach more preferable weight loss effect, it is allowed to which reliability index is less than 1, might as well be set to η_j=0.95 j=1,2 ..., 5, therefore the constraints based on reliability is converted to：

Poss(T_j(α^I, x)≤150 DEG C) and >=0.95 j=1,2 ..., 5

Wherein Poss represents the probability that inequality is set up, and can specifically be solved by following interval possibility degree calculation formula：

$\mathrm{Poss}(T_j({\mathrm{\&alpha;}}^I,x)\&le;150)=\left\{\begin{array}{cc}1& {\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&le;150\\ \frac{150-{\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)}{{\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)-{\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)}& {\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&le;150\&le;{\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\\ 0& {\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&GreaterEqual;150\end{array}\right.$

Wherein

T _j(α^I, it is respectively x) upper bound of main structure temperature-responsive and lower bound under each characteristic time, i.e.,：

${\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)=\underset{\mathrm{\&alpha;}\&Element;{\mathrm{\&alpha;}}^I}{\max }{T}_j(\mathrm{\&alpha;},x)$ ${\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)=\underset{\mathrm{\&alpha;}\&Element;{\mathrm{\&alpha;}}^I}{\max }{T}_j(\mathrm{\&alpha;},x)j=\mathrm{1,2},...,5$

Pass through the processing of this step so that the Optimized model set up in step 2 is converted to the nested Optimized model of the complexity containing reliability index, including inside and outside bilevel optimization, wherein outer layer optimizes for design vector x=(x₁,x₂,x₃,x₄)^TOptimizing, can be realized by step 6；And internal layer optimization is then used to calculate temperature field T_j(α^I, x) on block design parameter alpha^IThe bound of response range, can pass through step 5 Equivalent realization；

Step 5：The rapid solving of temperature field response range, in addition to optimization method, can also try to achieve temperature field on block design parameter alpha by interval numerical computations^IResponse range.For the solution of the algebraic equation containing interval parameter, traditional perturbation method can approximately try to achieve response range, but due to only remaining linear term during matrix inversion, therefore often bring than larger deviation.The present invention is by means of improved Neumann expansion techniques, establish the interval Finite Volume Method suitable for the equation of heat conduction, thermal protection system temperature field response range can fast and accurately be determined, optimize so as to instead of the internal layer described in step 4 in nested Optimized model with interval arithmetic, become two layers of nested optimization problem into conventional individual layer optimization problem, substantially increase computational efficiency.Specific implementation method is as follows：

First by the structural model discretization shown in Fig. 2, it is Δ x=2mm to make spatial mesh size, using six point symmetry forms of second order accuracy, to ensure its stability, time step is taken as Δ t=0.2s, can then obtain the following transient state temperature field limited bulk algebraic equation containing interval parameter：

A(α^I)T^k+1=B (α^I)T^k+F(α^I)

Wherein T^kRepresent temperature vector at all nodes of kth time step.

The coefficient matrix and right-hand-side vector of above-mentioned equation are carried out into Taylor expansion at parameter intermediate value to obtain：

$A\left({\mathrm{\&alpha;}}^I\right)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=A^c+\mathrm{\&Delta;}{A}^I$

$B\left({\mathrm{\&alpha;}}^I\right)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;B}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;B}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=B^c+\mathrm{\&Delta;}{B}^I$

$F\left({\mathrm{\&alpha;}}^I\right)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;F}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=F^c+\mathrm{\&Delta;}{F}^I$

Further understood using Neumann series：

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i^I{\left(A^c\right)}^{-1})}^r$

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i)}^r$

Wherein

When the symbolic variable r in sum formula takes different value, in above formula

Item can specifically deploy, and can be obtained after merging similar terms：

Ensureing norm | | Δ α_iA_i| | on the premise of ＜ 1 is set up, then level is several

It is convergent.Therefore, cast out the cross term in above formula, can obtain：

${(A^c+{\mathrm{\&Delta;A}}^I)}^{-1}\&ap;{\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\left({-\mathrm{\&Delta;\&alpha;}}_i^I{A}_i\right)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^I$

Wherein $E_i^I=\frac{{-\mathrm{\&Delta;\&alpha;}}_i^I{A}_i}{I+{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i}$

It is updated to formula：

(A^c+ΔA^I)((T^k+1)^c+Δ(T^k+1)^I)=(B^c+ΔB^I)((T^k)^c+Δ(T^k)^I)+(F^c+ΔF^I)

In, using the basic operation rule of intervl mathematics, it can obtain：

${\left(T^{k+1}\right)}^c={\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^c)[B^c{\left(T^k\right)}^c+F^c]$

$\mathrm{\&Delta;}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i^I(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}{\left(T^k\right)}^I+{\mathrm{\&Delta;B}}^I{\left(T^k\right)}^c+{\mathrm{\&Delta;F}}^I)]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i\&CenterDot;\mathrm{\&delta;}\&CenterDot;(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}\left(T^k\right)\&CenterDot;\mathrm{\&delta;}+\mathrm{\&Delta;B}\&CenterDot;\mathrm{\&delta;}{\&CenterDot;\left(T^k\right)}^c+\mathrm{\&Delta;F}\&CenterDot;\mathrm{\&delta;})]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}\left(T^k\right)+\mathrm{\&Delta;B}{\left(T^k\right)}^c+\mathrm{\&Delta;F})]\&CenterDot;\mathrm{\&delta;}$

$=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\mathrm{\&delta;}=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;[-\mathrm{1,1}]$

Wherein：

$\mathrm{\&Delta;}\left(T^{k+1}\right)=|\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{\left(A^c\right)}^{-1}{\mathrm{\&Delta;E}}_i(B^c{\left(T^k\right)}^c+F^c)+{\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{10}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c{\mathrm{\&Delta;T}}^k+\mathrm{\&Delta;B}{\left(T^k\right)}^c+\mathrm{\&Delta;F})|$

Then the interval bound of whole thermal protection system transient state temperature field response is：

${\stackrel{\&OverBar;}{T}}^{k+1}={\left(T^{k+!}\right)}^c+\mathrm{\&Delta;}\left(T^{k+1}\right)$ T ^k+1=(T^k+1)^c-Δ(T^k+1)

Further, extract lower surface nonwoven fabric layer and the temperature maximum at all common nodes of main structure, the upper bound of main structure temperature-responsive and lower bound under each characteristic time can be quickly obtained, the internal layer that be instead of with interval arithmetic described in step 4 in nested Optimized model optimizes.Due to remaining part higher order term during matrix inversion, therefore compared with traditional Novel Interval Methods, computational accuracy is significantly improved.

Step 6：The solution of deterministic optimization problem, by being handled in step 5 simplifying for reliability constraint, object function takes the nominal value of construction weight simultaneously, former double-deck uncertain optimization problem can be converted into the conventional single layer deterministic optimization problem containing reliability index, i.e.,：

$\underset{x}{\min }M({\mathrm{\&alpha;}}^c,x)$

s.t.Poss(T_j(α^I, x)≤150 DEG C) and >=0.95 j=1,2 ..., 5

10mm≤x₃≤40mm

2mm≤x_i≤ 15mm i=1,2,4

Using simulated annealing, calculation procedure is write.Relation between being expended in view of computational accuracy and calculating, defines maximum cycle Iter_max=3000 and convergence factor ε=10^-4, when any one in following 3 conditions is met, calculates and terminate：

（1）Loop iteration frequency n ＞ Iter_max；

（2）In double iterative process, object function relative variation is met $\left|\frac{M^{(i+1)}\left({\mathrm{\&alpha;}}^c\right)-M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}{M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}\right|<\mathrm{\&epsiv;};$

（3）   ||x⁽ⁱ⁺¹⁾-x⁽ⁱ⁾||₂＜ ε

Wherein | | | |₂Represent 2 norms of vector.

When reaching condition（1）When, the new initial value of design variable is given, and be brought into algorithm and recalculate；When algorithm is because of condition（2）Or（3）During termination, the result of calculation x of ith iteration process is taken⁽ⁱ⁾It is used as the optimal value of design variable.

Among the present embodiment, by 847 iterative calculation, the end condition shown in above-mentioned 2nd article is reached, the dimensionally-optimised design process of thermal protection system based on section reliability has been completed, optimal solution is respectively x₁=5.1mm, x₂=3.2mm, x₃=18.3mm, x₄=2.7mm, now the nominal value of construction weight is 1.534kg, has reached optimal weight loss effect.

Above-described is only presently preferred embodiments of the present invention, and the present invention is not limited solely to above-described embodiment, and local change, equivalent substitution, improvement for being made within the spirit and principles of the invention etc. should be included in the scope of the protection.

Claims (2)

1. a kind of size optimization design method for aircraft thermal protection system based on section reliability, it is characterised in that comprise the following steps：

Step one：It is determined that the design parameter for optimizing the Basic Design variable of the aircraft thermal protection system of design and correlation is needed, wherein the Basic Design variable x=(x₁,x₂,x₃,x₄)^T, x₁、x₂、x₃、x₄Radiation coating, upper surface nonwoven fabric layer, heat insulation layer, the thickness of lower surface nonwoven fabric layer are represented respectively, according to actual requirement of engineering, determine the initial range of above design variable；The design parameter includes the hot property parameters of material, outside heat and carried；Convenient for statement, all design parameters are uniformly written as vectorial α=(α₁,...,α_n)^TForm, wherein n represents the quantity of parameter；

Step 2：Using thermal protection system oeverall quality as the object function of optimization, the temperature in use scope of aircraft body structure sets up such as next Non-linear Optimal Model as constraints：

$\underset{x}{\min }M(\mathrm{\&alpha;},x)$

s.t.T_j(α,x)≤T^maxJ=1,2 ..., m

$\underset{\&OverBar;}{x}\&le;x\&le;\stackrel{\&OverBar;}{x}$

(α x) represents structure gross mass to wherein M；T_j(α x) is temperature-responsive function；For T^maxThe temperature upper limit that can bear for main structure；M is the number of constraints；

It is the bound of the design variable initial range defined in step one；

Step 3：Take into full account the fluctuation that temperature change is brought to material properties, and hot-fluid the uncertain factor such as measurement error, the uncertain parameter in this thermal protection system is described with interval vector, i.e.,：

$\mathrm{\&alpha;}\&Element;{\mathrm{\&alpha;}}^I=[\underset{\&OverBar;}{\mathrm{\&alpha;}},\stackrel{\&OverBar;}{\mathrm{\&alpha;}}]=[{\mathrm{\&alpha;}}^c-\mathrm{\&Delta;\&alpha;},{\mathrm{\&alpha;}}^c+\mathrm{\&Delta;\&alpha;}]={\mathrm{\&alpha;}}^c+[-\mathrm{\&Delta;\&alpha;},\mathrm{\&Delta;\&alpha;}]={\mathrm{\&alpha;}}^c+\mathrm{\&Delta;\&alpha;\&delta;}$

Wherein

The upper bound and the lower bound of vector, α are represented respectively^c, Δ α is respectively the nominal value and radius of vector, and is met：

${\mathrm{\&alpha;}}^c={({{\mathrm{\&alpha;}}_1}^c,...,{{\mathrm{\&alpha;}}_n}^c)}^T={(({\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_1+{\underset{\&OverBar;}{\mathrm{\&alpha;}}}_1)/2,...,({\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_n+{\underset{\&OverBar;}{\mathrm{\&alpha;}}}_n)/2)}^T$

$\mathrm{\&Delta;\&alpha;}={(\mathrm{\&Delta;}{\mathrm{\&alpha;}}_1,...,\mathrm{\&Delta;}{\mathrm{\&alpha;}}_n)}^T={(({\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_1-{\underset{\&OverBar;}{\mathrm{\&alpha;}}}_1)/2,...,({\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_n-{\underset{\&OverBar;}{\mathrm{\&alpha;}}}_n)/2)}^T;$

δ=[- 1,1]

Step 4：Constraints conversion based on section reliability, when the uncertainty of system uses interval vector α^IDuring portraying, the temperature-responsive T of aircraft body structure in step 2_j(α, x) will be by its interval function T_j(α^I, x) substituted, to improve the safety and stability that structure is used, according to the definition of reliability, the certain requirement of the possibility sexual satisfaction that temperature constraint is set up required in the design phase, i.e.,：

Poss(T_j(α^I,x)≤T^max)≥η_j

Wherein η_jFor reliability index, value is between 0 to 1；Poss represents the probability that condition is set up；Assuming that interval variable is equally distributed in its given range, then this probability can be tried to achieve by following interval possibility degree calculation formula：

$\mathrm{Poss}(T_j({\mathrm{\&alpha;}}^I,x)\&le;T^{\max })=\left\{\begin{array}{cc}1& {\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&le;T^{\max }\\ \frac{T^{\max }-{\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)}{{\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)-{\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)}& {\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&le;T^{\max }\&le;{\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\\ 0& {\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)\&GreaterEqual;T^{\max }\end{array}\right.$

Wherein

T _j(α^I, x) it is respectively interval function T_j(α^I, upper bound x) and lower bound, i.e.,：

${\stackrel{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)=\underset{\mathrm{\&alpha;}\&Element;{\mathrm{\&alpha;}}^I}{\max }{T}_j(\mathrm{\&alpha;},x)$ ${\underset{\&OverBar;}{T}}_j({\mathrm{\&alpha;}}^I,x)=\underset{\mathrm{\&alpha;}\&Element;{\mathrm{\&alpha;}}^I}{\max }{T}_j(\mathrm{\&alpha;},x)$

Pass through the processing of this step so that the Optimized model set up in step 2 is converted to the nested Optimized model of the complexity containing reliability index, including inside and outside bilevel optimization, wherein outer layer optimizes for design vector x=(x₁,x₂,x₃,x₄)^TOptimizing, can be realized by step 6；And internal layer optimization is then used to calculate temperature field T_j(α^I, x) on block design parameter alpha^IThe bound of response range, can pass through step 5 Equivalent realization；

Step 5：The rapid solving of temperature field response range, in addition to optimization method, temperature field is on block design parameter alpha^IResponse range can be tried to achieve by interval numerical computations, initially set up the limited configurations volume-based model of thermal protection system, using six point symmetry discrete schemes, the following transient state temperature field limited bulk equation containing interval parameter can be obtained：

A(α^I)T^k+1=B (α^I)T^k+F(α^I)

Wherein T^kRepresent temperature vector at all nodes of kth time step；

The coefficient matrix and right-hand-side vector of above-mentioned equation are carried out into Taylor expansion at parameter nominal value to obtain：

$A\left({\mathrm{\&alpha;}}^I\right)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=A^c+\mathrm{\&Delta;}{A}^I$

$B\left({\mathrm{\&alpha;}}^I\right)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;B}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;B}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=B^c+\mathrm{\&Delta;}{B}^I$

$F\left({\mathrm{\&alpha;}}^I\right)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;F}{\&PartialD;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;F}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i{\mathrm{\&delta;}}_i=F^c+\mathrm{\&Delta;}{F}^I$

Further understood using Neumann series：

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{\&PartialD;A}{{\&PartialD;\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_i^I{\left(A^c\right)}^{-1})}^r$

$={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i)}^r$

Wherein

When the symbolic variable r in sum formula takes different value, in above formula

Item can specifically deploy, and can be obtained after merging similar terms：

If norm condition | | Δ α_iA_i| | ＜ 1 is set up, then level is several

It is convergent, therefore, if casting out the cross term in above formula, can obtains：

${(A^c+{\mathrm{\&Delta;A}}^I)}^{-1}\&ap;{\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{r=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{\left({-\mathrm{\&Delta;\&alpha;}}_i^I{A}_i\right)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^I$

Wherein $E_i^I=\frac{{-\mathrm{\&Delta;\&alpha;}}_i^I{A}_i}{I+{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i};$

It is updated to formula：

(A^c+ΔA^I)((T^k+1)^c+Δ(T^k+1)^I)=(B^c+ΔB^I)((T^k)^c+Δ(T^k)^I)+(F^c+ΔF^I)

In, using the basic operation rule of intervl mathematics, it can obtain：

${\left(T^{k+1}\right)}^c={\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^c)[B^c{\left(T^k\right)}^c+F^c]$

$\mathrm{\&Delta;}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i^I(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}{\left(T^k\right)}^I+{\mathrm{\&Delta;B}}^I{\left(T^k\right)}^c+{\mathrm{\&Delta;F}}^I)]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i\&CenterDot;\mathrm{\&delta;}\&CenterDot;(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}\left(T^k\right)\&CenterDot;\mathrm{\&delta;}+\mathrm{\&Delta;B}\&CenterDot;\mathrm{\&delta;}{\&CenterDot;\left(T^k\right)}^c+\mathrm{\&Delta;F}\&CenterDot;\mathrm{\&delta;})]$

$={\left(A^c\right)}^{-1}[\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{E}_i(B^c{\left(T^k\right)}^c+F^c)+(I+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c\mathrm{\&Delta;}\left(T^k\right)+\mathrm{\&Delta;B}{\left(T^k\right)}^c+\mathrm{\&Delta;F})]\&CenterDot;\mathrm{\&delta;}$

$=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\mathrm{\&delta;}=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;[-\mathrm{1,1}]$

Wherein：

$\mathrm{\&Delta;}\left(T^{k+1}\right)=|\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(A^c\right)}^{-1}{\mathrm{\&Delta;E}}_i(B^c{\left(T^k\right)}^c+F^c)+{\left(A^c\right)}^{-1}(I+\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{E}_i^c)(B^c{\mathrm{\&Delta;T}}^k+\mathrm{\&Delta;B}{\left(T^k\right)}^c+\mathrm{\&Delta;F})|$

Then the interval bound of thermal protection structure transient temperature response is：

${\stackrel{\&OverBar;}{T}}^{k+1}={\left(T^{k+!}\right)}^c+\mathrm{\&Delta;}\left(T^{k+1}\right)$ T ^k+1=(T^k+1)^c-Δ(T^k+1)

Using the interval Finite Volume Method proposed in this step, thermal protection structure temperature field response range can be quickly determined, so that the internal layer that be instead of with interval arithmetic described in step 4 in nested Optimized model optimizes, computational efficiency is improved；

Step 6：The solution of deterministic optimization problem, interval limited bulk side in step 5 is converted into individual layer deterministic optimization problem to the rapid solving of structure transient state temperature field, the former nested optimization problem containing interval parameter；Using simulated annealing, calculation procedure is write, maximum cycle Iter is defined_maxWith convergence factor ε, when any one in following three condition is met, calculates and terminate：

（1）Loop iteration frequency n ＞ Iter_max；

（2）In double iterative process, object function nominal value relative variation is met：

$\left|\frac{M^{(i+1)}\left({\mathrm{\&alpha;}}^c\right)-M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}{M^{\left(i\right)}\left({\mathrm{\&alpha;}}^c\right)}\right|<\mathrm{\&epsiv;};$

（3）   ||x⁽ⁱ⁺¹⁾-x⁽ⁱ⁾||₂＜ ε

Wherein | | | |₂Represent 2 norms of vector；

When reaching condition（1）When, the new initial value of design variable is given, and be brought into algorithm and recalculate；When algorithm is because of condition（2）Or（3）During termination, the result of calculation x of ith iteration process is taken⁽ⁱ⁾As the optimal value of design variable, the dimensionally-optimised design process of thermal protection system based on section reliability is completed, with the weight loss effect being optimal.

2. a kind of size optimization design method for aircraft thermal protection system based on section reliability according to claim 1, it is characterised in that：Required main structure temperature in use scope is no more than 150 DEG C, i.e. T in the step 2^max=150 DEG C；Engineering requirements design variable meets 2mm≤x₁,x₂,x₄≤15mm 10mm≤x₃≤40mm。

Images