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

Abstract

Based on a size optimization design method for aircraft thermal protection system for section reliability, step is as follows: 1, mathematical modeling, according to thermal protection system version, determines the design variable that size is relevant, sets up thermal protection struc ture loss of weight Optimized model according to temperature requirement; 2, the uncertainty in quantitative description material, external applied load and just/boundary value condition is carried out in utilization interval; 3, based on interval possibility degree, reliability conversion is carried out to the temperature constraint containing uncertain parameter; 4, utilize interval Finite Volume Method, realize the rapid solving of structure transient temperature responding range, become two-layer nested optimization problem into conventional individual layer optimization problem; 5, program calculation is carried out to the individual layer optimization problem containing reliability constraint after conversion, determine optimum solution, to reach best weight loss effect.The present invention systematization can solve the thermal protection system structural optimization problems containing multi-source bounded-but-unknown uncertainty, improves counting yield, ensure that the reliability that structure uses and security.

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, particularly a kind of thermal protection system sizing method optimized based on interval heat conduction analysis and section reliability.

Background technology

Along with the fast development of Chinese Space technology, spacecraft is entered the orbit or is returned to ground, get through the earth's atmosphere, and especially will suffer serious Aerodynamic Heating effect at return phase.The various piece of spacecraft structure has different requirements to temperature, and most of instrument and equipment requires to be in certain temperature range, and thermal protection system therefore must be used to carry out temperature control to spacecraft.In general, thermal protection system thickness is larger, and thermal protection ability is better; But on the other hand, the architecture quality increase that thickness increase brings 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.

As everyone knows, uncertainty is extensively present in objective world, aircraft thermal protection structure system is in the impact of manufacturing and in use procedure, its thermal force, physical dimension, material behavior etc. are inevitably subject to the uncertain factors such as production defect, measuring error, these all can have an impact to the calculating of thermal protection system heat transfer process, cause structure temperature field distribution to produce fluctuation, even occur the possibility lost efficacy.Traditional thermal protection system thermodynamic analysis and optimal design are all implemented based on deterministic models, can not embody practical problems and contain probabilistic objective essence, usually these design proposals can cause the increase of architecture quality and certain unsafe factor.

For reducing the impact of various uncertainty on thermal protection system heat-proof quality as far as possible, deviser should just predict contingent change in the design phase, and take effective computing method and design proposal, strengthen the insensitivity fluctuated to Parameters variation in temperature field, thus improve the safety in utilization of primary aircraft structure and instrument and equipment, Here it is based on reliability theory carries out the original intention of size optimization design method for aircraft thermal protection system research.For the problem of many parameter uncertainty quantification of practical structures, enough sample informations be obtained, construct its accurate probability distribution function or fuzzy membership function and often to seem very difficult or high cost.And the interval model uncertain quantitative method that to be a class relatively new, only by the bound of less information acquisition variable, therefore need embody better economy in modeling.In addition, becoming more meticulous of structural reliability problem research is required to make tradition seem too conservative based on the design concept of method of safety coefficients.Therefore, various probabilistic impact is just taken into full account from the establishment stage of calculating and Optimized model, size optimization design method for aircraft thermal protection system based on section reliability is proposed, for the deficiency making up existing thermal protection system transient state temperature field numerical evaluation and optimal design, there is important engineer applied and be worth.

Summary of the invention

Technology of the present invention is dealt with problems: overcome the deficiency that existing designing technique exists in the structure optimization of aircraft thermal protection system, a kind of size optimization design method for aircraft thermal protection system based on section reliability is provided, interval Finite Volume Method is incorporated in the numerical evaluation of the thermal protection system transient state temperature field containing interval uncertain parameter, under ensureing that temperature field meets the prerequisite of primary aircraft structure request for utilization, obtain a kind of physical dimension Robust-Design scheme reducing overall system quality.

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: determine the Basic Design variable of the aircraft thermal protection system needing to be optimized design and relevant design parameter, design variable x=(x ₁, x ₂, x ₃, x ₄) ^trepresent, wherein:

X ₁, x ₂, x ₃, x ₄represent the thickness of radiation coating, upper surface nonwoven fabric layer, heat insulation layer, lower surface nonwoven fabric layer respectively; According to actual requirement of engineering, determine the initial range of above design variable;

Design parameter comprises the physical attribute of three kinds of materials, as density of material ρ _i, coefficient of heat conductivity k _i, specific heat c _ii=1,2,3; The heat born is carried, and represents with heat flow density q; For simplicity, the correlation parameter in this calculating and Optimized model is expressed as the form of vectorial α, that is:

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

Step 2: using thermal protection system oeverall quality as the objective function optimized, the serviceability temperature scope of primary aircraft structure, as constraint condition, sets up following Optimized model:

$\begin{array}{cc}\underset{x}{min}& M(\mathrm{\α},x)\end{array}$

s.t.T _j(α,x)≤T ^maxj＝1,2,...,m

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

Wherein M (α, x) represents structure gross mass; T ^maxfor the upper temperature limit that main structure can bear; M is the number of constraint condition; 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 change with temperature, take into full account the fluctuation that temperature variation is brought to material properties, and the uncertain factor such as the measuring error of hot-fluid, each uncertain parameter in this thermal protection system is described, that is: with interval vector

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

Wherein represent the upper bound and the lower bound of vector respectively, α ^c, Δ α is respectively nominal value and the radius of vector, and meets:

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

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

δ＝[-1,1]

Step 4: the constraint condition based on section reliability is changed, 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 constraint condition _j(α, x) is no longer traditional fixed function, and can by interval function T _j(α ⁱ, x) substitute; For improving the safety and stability that structure uses, according to the definition of reliability, need to ensure that the possibility that constraint condition is set up meets certain requirement, might as well by the reliability index η to temperature restraint _jrepresent, the constraint condition so 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.Suppose that interval variable is equally distributed in the little perturbation range of its nominal value, then this probability calculates by general interval possibility degree computing formula:

$Poss\left(T_j\right({\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 be respectively interval function T _j(α ⁱ, upper bound x) and lower bound, that is:

$\begin{array}{cc}{\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}{min}{T}_j(\mathrm{\α},x)\end{array}$

By the process of this step, make the Optimized model set up in step 2 be converted to the nested Optimized model of complexity containing reliability index, comprise inside and outside bilevel optimization, its ectomesoderm optimization is used for design vector x=(x ₁, x ₂, x ₃, x ₄) ^toptimizing, realize by step 6; Internal layer optimization is then for accounting temperature field T _j(α ⁱ, x) about block design parameter alpha ⁱthe bound of responding range, by step 5 Equivalent realization;

Step 5: the rapid solving of temperature field responding range, except optimization method, can also try to achieve temperature field about block design parameter alpha by interval numerical evaluation ⁱresponding range.For solving containing interval parameter algebraic equation, traditional perturbation method can be similar to and try to achieve responding range, but owing to only remaining linear term in matrix inversion process, therefore often brings larger deviation.The present invention is by means of the Neumann expansion technique improved, establish the interval Finite Volume Method being applicable to heat-conduction equation, thermal protection system temperature field responding range can be determined fast and accurately, thus the internal layer optimization that instead of with interval arithmetic in nested Optimized model described in step 4, become two-layer nested optimization problem into conventional individual layer optimization problem, substantially increase optimization counting yield.Specific implementation method is as follows:

First set up the limited configurations volume discrete model of thermal protection system, adopt six point symmetry forms of second order accuracy, the transient state temperature field limited bulk algebraic equation as follows containing interval parameter can be obtained:

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

Wherein T ^krepresent all Nodes temperature vectors of kth time step;

The matrix of coefficients of above-mentioned equation and right-hand-side vector are carried out Taylor expansion at parameter intermediate value place can obtain:

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

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

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

Utilize Neumann progression known further:

$\begin{array}{c}{(A^c+{\mathrm{\ΔA}}^I)}^{-1}={\left(A^c\right)}^{-1}+\underset{r=1}{\overset{\mathrm{\∞}}{\Σ}}{\left(A^c\right)}^{-1}{(-{\mathrm{\ΔA}}^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}{\Σ}}\frac{\∂A}{\∂{\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}{\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{\Δ\α}}_i^I{A}_i)}^r\end{array}$

Wherein $A_i=\frac{\∂A}{\∂{\mathrm{\α}}_i}{|}_{{\mathrm{\α}}^c}{\left(A^c\right)}^{-1}$

When the symbolic variable r in sum formula gets different value, in above formula item can specifically launch, and can obtain after merging like terms:

In guarantee norm || Δ α _ia _i|| under the prerequisite that <1 sets up, then level is several it is convergence; 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{n}{\&Sigma;}}\underset{r=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^I$

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

Be 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, utilize the fundamental operation rule of intervl mathematics, can obtain:

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

$\begin{array}{c}\mathrm{\&Delta;}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}\&lsqb;\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)\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\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)\left(B^c\mathrm{\&Delta;}\right(T^k)\&CenterDot;\mathrm{\&delta;}+\mathrm{\&Delta;}B\&CenterDot;\mathrm{\&delta;}\&CenterDot;{\left(T^k\right)}^c+\mathrm{\&Delta;}F\&CenterDot;\mathrm{\&delta;})\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\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)\left(B^c\mathrm{\&Delta;}\right(T^k)+\mathrm{\&Delta;}B{\left(T^k\right)}^c+\mathrm{\&Delta;}F)\&rsqb;\&CenterDot;\mathrm{\&delta;}\\ =\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\mathrm{\&delta;}=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\&lsqb;-1,1\&rsqb;\end{array}$

Wherein:

$\mathrm{\&Delta;}\left(T^{k+1}\right)=|\underset{i=1}{\overset{n}{\&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}{\&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:

$\begin{array}{cc}{\stackrel{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c+\mathrm{\&Delta;}\left(T^{k+1}\right)& {\underset{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c-\mathrm{\&Delta;}\left(T^{k+1}\right)\end{array}$

In this step, conventional finite volume method is combined with interval mathematical theory, establish a kind of interval Finite Volume Method being applicable to calculate containing interval parameter transient heat conduction; By retaining the part higher order term in Neumann progression, improve computational accuracy greatly when calculated amount allows, this just provides a strong instrument for the approximate processing of internal layer optimization problem in nested Optimized model;

Step 6: solving of deterministic optimization problem, according to the interval Finite Volume Method in step 5 to the rapid solving of structure transient state temperature field, the former nested optimization problem containing interval parameter is converted into individual layer deterministic optimization problem; Adopt simulated annealing, write calculation procedure, definition maximum cycle Iter _maxwith converging factor ε, when any one in following three conditions is met, calculates and stop:

(1) loop iteration frequency n >Iter _max;

(2) in double iterative process, objective function nominal value relative variation meets:

$\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), the initial value that given design variable is new, and be brought in algorithm and recalculate; When algorithm stops because of condition (2) or (3), get the result of calculation x of i-th iterative process ⁽ⁱ⁾as the optimal value of design variable, complete the dimensionally-optimised design process of thermal protection system based on section reliability, to reach optimum weight loss effect.

The present invention's advantage is compared with prior art:

(1) compared with traditional structural optimization problems, the Optimized model set up takes into full account that the multi-source that the complicated thermal force environment of aircraft thermal protection system causes is uncertain, thus improving thermal protection system safety in utilization and stability, result of calculation has prior directive significance to its structural design.

(2) characterize uncertain factor in calculating and optimization problem by interval model, greatly reduce traditional probabilistic model and fuzzy model to the rigors of amounts of specimen information.

(3) by retaining the higher order term in Neumann progression, establish the interval Finite Volume Method being applicable to heat conduction analysis, the responding range of each moment structure temperature field can be determined fast, thus 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.

Accompanying drawing explanation

Fig. 1 is containing interval parameter thermal protection system reliability Optimum Design flow process;

Fig. 2 aircraft thermal protection system structural representation;

Fig. 3 constraint condition reliable realization principle schematic.

Embodiment

Below in conjunction with drawings and Examples, the present invention will be further described.

In order to introduce the present invention in detail, first introduce definition and the algorithm thereof of the interval analysis operation used in the present invention.If R is real number field, for given two real numbers and then:

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

Be called bounded closed interval, be also interval number, be called for short interval.Wherein claim for lower bound or the lower extreme point in interval, claim for the upper bound or the upper extreme point in interval.If two intervals with corresponding end points is up and down equal respectively, then claim these two intervals equal, even and then x ⁱ=y ⁱ.Claim in addition with be respectively interval x ⁱnominal value and radius.

For two intervals arbitrary in real number field its interval arithmetic is defined as:

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

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

$x^I\&CenterDot;y^I=\&lsqb;\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}\&rsqb;\&CenterDot;\&lsqb;\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}\&rsqb;=\&lsqb;\min \{\underset{\&OverBar;}{x}\underset{\&OverBar;}{y},\underset{\&OverBar;}{x}\stackrel{\&OverBar;}{y},\stackrel{\&OverBar;}{x}\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{xy}\},\max \{\underset{\&OverBar;}{x}\underset{\&OverBar;}{y},\underset{\&OverBar;}{x}\stackrel{\&OverBar;}{y},\stackrel{\&OverBar;}{x}\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{xy}\}\&rsqb;$

$\begin{array}{cc}{x}^I/y^I=\&lsqb;\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}\&rsqb;/\&lsqb;\underset{\&OverBar;}{y},\stackrel{\&OverBar;}{y}\&rsqb;=\&lsqb;\underset{\&OverBar;}{x},\stackrel{\&OverBar;}{x}\&rsqb;\&CenterDot;\&lsqb;1/\stackrel{\&OverBar;}{y},1/\underset{\&OverBar;}{y}\&rsqb;& 0\&NotElement;y^I\end{array}$

Introduce the size optimization design method for aircraft thermal protection system based on section reliability below in detail:

The present invention is applicable to the dimensionally-optimised problem of thermal protection system containing interval uncertain parameter.The ceramic blanket thermal protection system scheme " AFRSI " that present embodiment proposes for NASAAMES research centre, illustrates proposed section reliability optimization method.In addition, the section reliability optimization method of this aircraft thermal protection system can be generalized in other uncertain optimization design containing the labyrinth of interval parameter.

The structural model of this thermal protection system as shown in Figure 2, considers that size is the daughter board of 0.3m × 0.3m.This structure is followed successively by from outside surface to main structure: C9 protects radiation coating, AB312 upper surface nonwoven fabric layer, Q-fiberFelt insulation material layer, AB312 lower surface nonwoven fabric layer, take hot-fluid as external heat load.Structure adopts hexahedral element to carry out discretize, extracts the actual work temperature of temperature maximum as main structure at lower surface nonwoven fabric layer and all common node places of main structure.

The section reliability optimizing process of this thermal protection system as shown in Figure 1, on the basis of traditional Optimized model, take into full account the uncertainty of system itself and external applied load, interval is utilized to carry out quantitative description to uncertain parameter, introduce reliability index and set up the transformation model of constraint condition in optimization problem based on interval possibility degree, utilizing the interval Finite Volume Method of proposition to try to achieve the bound of structure transient state temperature field response fast simultaneously.Adopt simulated annealing, write calculation procedure, can require to choose optimum thermal protection system design size according to deviser.Following several step can be divided into carry out:

Step one: determine the Basic Design variable of the thermal protection system needing to be optimized design and relevant design parameter, design variable x=(x ₁, x ₂, x ₃, x ₄) ^trepresent, wherein:

X ₁, x ₂, x ₃, x ₄: the thickness representing radiation coating, upper surface nonwoven fabric layer, heat insulation layer, lower surface nonwoven fabric layer 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 self dimensional requirement, that is:

10mm≤x ₃≤40mm

2mm≤x _i≤15mmi＝1,2,4

In this thermal protection system structural model, the density p of C9 coating material ₁=2.0 × 10 ³kg/m ³, coefficient of heat conductivity k ₁=0.64W/ (m DEG C), specific heat c ₁=628J/ (kg DEG C); The density p of AB312 nonwoven fabric layer material ₂=985.15kg/m ³, coefficient of heat conductivity k ₂=7.1 × 10 ^-2w/ (m DEG C), specific heat c ₂=558J/ (kg DEG C); The density p of Q-fiberFelt thermal insulation material ₃=56.1kg/m ³, coefficient of heat conductivity k ₃=8.2 × 10 ^-3w/ (m DEG C), specific heat c ₃=605J/ (kg DEG C); The time dependent function of heat flow density is q (t)=[12000-(t-800) ²/ 60] W/m ²;

For simplicity, the correlation parameter in this calculating and Optimized model is expressed as the form of vectorial α, that is:

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

Step 2: in transient heat conduction process, the extraction characteristic time (is respectively 200s, 500s, 800s, 1200s, 1600s) temperature maximum at lower surface nonwoven fabric layer and all common node places of main structure as the actual work temperature of main structure, use T respectively _j(α, x) j=1,2 ..., 5 represent.Constraint condition is taken as the maximum temperature of main structure in the working time and is no more than 150 DEG C, i.e. T ^max=150 DEG C.Under this temperature constraint, with the gross mass M (α, x) of thermal protection system structure for design object, can set up as next Optimized model:

$\begin{array}{cc}\underset{x}{min}& M(\mathrm{\&alpha;},x)\end{array}$

s.t.T _j(α,x)≤T ^maxj＝1,2,...,5；

10mm≤x ₃≤40mm

2mm≤x _i≤15mmi＝1,2,4

Step 3: due in Transient Heat Transfer process, the hot attribute of structured material can change with temperature, and temperature variation brings certain fluctuation to material properties, and the measurement of hot-fluid simultaneously exists certain error.That considers in Practical Project problem is less about probabilistic quantity of information, utilizes the interval each uncertain parameter describing this thermal protection system in the present invention.There is no harm in the perturbation that each design parameter listed in setting procedure one exists 5% near its nominal value, that is:

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

Wherein Δ α=0.05 α ^c;

Step 4: the constraint condition based on section reliability is changed, under considering that deviser can tolerate the prerequisite destroyed to a certain degree to constraint condition, for the reliability index that deviser provides, utilize interval possibility degree computing formula, setting up the transformation model of constraint condition, taking into account under various parameter fluctuation change condition, making design point still in feasible zone, meet the requirement of reliability, as shown in Figure 3.Transverse and longitudinal coordinate x ₁and x ₂represent two design variables in two-dimensional space respectively, the critical conditions unification g of constraint condition ₁(x)=0 and g ₂x ()=0 represents.The optimum solution that in figure, A and B represents traditional optimum solution respectively and obtain based on reliability optimization, solid line and dotted line represent tradition optimization and the feasible zone border corresponding to reliability optimization respectively.Can find out, tradition optimum solution A be often positioned at feasible zone border or its near, but due to the impact of uncertain factor, constraint condition can change, wherein a kind of situation is exactly that feasible zone border changes to dotted line from solid line, so traditional optimum solution A is positioned at outside new feasible zone, does not meet designing requirement; The optimum solution B obtained with reliability optimization then still meets the requirement of new constraint condition.In the present embodiment, in order to reach better weight loss effect, allowing reliability index to be less than 1, might as well η be set to _j=0.95j=1,2 ..., 5, the constraint condition therefore based on reliability is just converted to:

Poss(T _j(α ^I,x)≤150℃)≥0.95j＝1,2,...,5

Wherein Poss represents the probability that inequality is set up, and concrete can be solved by following interval possibility degree computing formula:

$Poss\left(T_j\right({\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 be respectively the upper bound and the lower bound of main structure temperature-responsive under each characteristic time, that is:

$\begin{array}{cc}{\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}{\min }{T}_j(\mathrm{\&alpha;},x)\end{array},j=1,2,\mathrm{...},5$

By the process of this step, make the Optimized model set up in step 2 be converted to the nested Optimized model of complexity containing reliability index, comprise inside and outside bilevel optimization, its ectomesoderm optimization is used for design vector x=(x ₁, x ₂, x ₃, x ₄) ^toptimizing, realize by step 6; Internal layer optimization is then for accounting temperature field T _j(α ⁱ, x) about block design parameter alpha ⁱthe bound of responding range, by step 5 Equivalent realization;

Step 5: the rapid solving of temperature field responding range, except optimization method, can also try to achieve temperature field about block design parameter alpha by interval numerical evaluation ⁱresponding range.For solving containing interval parameter algebraic equation, traditional perturbation method can be similar to and try to achieve responding range, but owing to only remaining linear term in matrix inversion process, therefore often brings larger deviation.The present invention is by means of the Neumann expansion technique improved, establish the interval Finite Volume Method being applicable to heat-conduction equation, thermal protection system temperature field responding range can be determined fast and accurately, thus the internal layer optimization that instead of with interval arithmetic in nested Optimized model described in step 4, become two-layer nested optimization problem into conventional individual layer optimization problem, substantially increase counting yield.Specific implementation method is as follows:

First by the structural model discretize shown in Fig. 2, making spatial mesh size be Δ x=2mm, adopting six point symmetry forms of second order accuracy, for ensureing its stability, time step is taken as Δ t=0.2s, so can obtain the transient state temperature field limited bulk algebraic equation as follows containing interval parameter:

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

Wherein T ^krepresent all Nodes temperature vectors of kth time step,

The matrix of coefficients of above-mentioned equation and right-hand-side vector are carried out Taylor expansion at parameter intermediate value place can obtain:

$A\left({\mathrm{\&alpha;}}^I\right)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\&Sigma;}}\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\&Sigma;}}\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&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}{\&Sigma;}}\frac{\&part;B}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\&Sigma;}}\frac{\&part;B}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&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}{\&Sigma;}}\frac{\&part;F}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{10}{\&Sigma;}}\frac{\&part;F}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&alpha;}}_i{\mathrm{\&delta;}}_i=F^c+{\mathrm{\&Delta;F}}^I$

Utilize Neumann progression known further:

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

Wherein $A_i=\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\left(A^c\right)}^{-1}$

When the symbolic variable r in sum formula gets different value, in above formula item can specifically launch, and can obtain after merging like terms:

In guarantee norm || Δ α _ia _i|| under the prerequisite that <1 sets up, then level is several it is convergence.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}{\&Sigma;}}\underset{r=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^I$

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

Be 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, utilize the fundamental operation rule of intervl mathematics, can obtain:

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

$\begin{array}{c}\mathrm{\&Delta;}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}\&lsqb;\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)\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\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)\left(B^c\mathrm{\&Delta;}\right(T^k)\&CenterDot;\mathrm{\&delta;}+\mathrm{\&Delta;}B\&CenterDot;\mathrm{\&delta;}\&CenterDot;{\left(T^k\right)}^c+\mathrm{\&Delta;}F\&CenterDot;\mathrm{\&delta;})\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\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)\left(B^c\mathrm{\&Delta;}\right(T^k)+\mathrm{\&Delta;}B{\left(T^k\right)}^c+\mathrm{\&Delta;}F)\&rsqb;\&CenterDot;\mathrm{\&delta;}\\ =\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\mathrm{\&delta;}=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\&lsqb;-1,1\&rsqb;\end{array}$

Wherein:

$\mathrm{\&Delta;}\left(T^{k+1}\right)=|\underset{i=1}{\overset{n}{\&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}{\&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:

$\begin{array}{cc}{\stackrel{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c+\mathrm{\&Delta;}\left(T^{k+1}\right)& {\underset{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c-\mathrm{\&Delta;}\left(T^{k+1}\right)\end{array}$

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

Step 6: solving of deterministic optimization problem, by in step 5 to the simplify processes of reliability constraint, objective function gets the nominal value of construction weight simultaneously, can be the conventional single layer deterministic optimization problem containing reliability index by former double-deck uncertain optimization question variation, that is:

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

s.t.Poss(T _j(α ^I,x)≤150℃)≥0.95j＝1,2,...,5

10mm≤x ₃≤40mm

2mm≤x _i≤15mmi＝1,2,4

Adopt simulated annealing, write calculation procedure.Consider computational accuracy and calculating expend between relation, definition maximum cycle Iter _max=3000 and converging factor ε=10 ^-4, when any one in following 3 conditions is met, calculates and stop:

(1) loop iteration frequency n >Iter _max;

(2) in double iterative process, objective function relative variation meets

(3)||x ⁽ⁱ⁺¹⁾-x ⁽ⁱ⁾|| ₂<ε

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

When reaching condition (1), the initial value that given design variable is new, and be brought in algorithm and recalculate; When algorithm stops because of condition (2) or (3), get the result of calculation x of i-th iterative process ⁽ⁱ⁾as the optimal value of design variable.

In the middle of the present embodiment, through 847 iterative computation, reach the end condition shown in above-mentioned 2nd article, complete the dimensionally-optimised design process of thermal protection system based on section reliability, optimum 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, reaches optimum weight loss effect.

Above-describedly be only preferred embodiment of the present invention, the present invention is not only confined to above-described embodiment, and all local done within the spirit and principles in the present invention are changed, equivalent replacement, improvement etc. all should be included within protection scope of the present invention.

Claims (2)

1., based on a size optimization design method for aircraft thermal protection system for section reliability, it is characterized in that comprising the following steps:

Step one: determine the Basic Design variable of the aircraft thermal protection system needing to be optimized design and relevant design parameter, wherein said Basic Design variable x=(x ₁, x ₂, x ₃, x ₄) ^t, x ₁, x ₂, x ₃, x ₄represent the thickness of radiation coating, upper surface nonwoven fabric layer, heat insulation layer, lower surface nonwoven fabric layer respectively, according to actual requirement of engineering, determine the initial range of above design variable; Described design parameter comprises the hot property parameters of material, external heat is carried; For stating conveniently, all design parameters are unified is written as vectorial α=(α ₁..., α _n) ^tform, wherein n represents the quantity of parameter;

Step 2: using thermal protection system oeverall quality as the objective function optimized, the serviceability temperature scope of primary aircraft structure, as constraint condition, is set up as next Non-linear Optimal Model:

$\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}$

Wherein M (α, x) represents structure gross mass; T _j(α, x) is temperature-responsive function; For T ^maxfor the upper temperature limit that main structure can bear; M is the number of constraint condition; x, it is the bound of the design variable initial range defined in step one;

Step 3: take into full account the fluctuation that temperature variation is brought to material properties, and these two uncertain factors of the measuring error of hot-fluid, describe the uncertain parameter in this thermal protection system with interval vector, that is:

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

Wherein αrepresent the upper bound and the lower bound of vector respectively, α ^c, △ α is respectively nominal value and the radius of vector, and meets:

${\mathrm{\&alpha;}}^c={({{\mathrm{\&alpha;}}_1}^c,...,{{\mathrm{\&alpha;}}_n}^c)}^T={\left(\right({\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;}\mathrm{\&alpha;}={({\mathrm{\&Delta;\&alpha;}}_1,...,{\mathrm{\&Delta;\&alpha;}}_n)}^T={\left(\right({\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: the constraint condition based on section reliability is changed, as the uncertainty interval vector α of system ⁱwhen portraying, the temperature-responsive T of primary aircraft structure in step 2 _j(α, x) will by its interval function T _j(α ⁱ, x) substitute, be improve structure use safety and stability, according to the definition of reliability, the design phase require temperature constraint set up possibility meet certain requirement, that is:

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; Suppose that interval variable is equally distributed in its given range, then this probability is tried to achieve by following interval possibility degree computing formula:

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

Wherein t _j(α ⁱ, x) be respectively interval function T _j(α ⁱ, upper bound x) and lower bound, that is:

$\begin{array}{cc}{\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}{\min }{T}_j(\mathrm{\&alpha;},x)\end{array}$

By the process of this step, make the Optimized model set up in step 2 be converted to the nested Optimized model of complexity containing reliability index, comprise inside and outside bilevel optimization, its ectomesoderm optimization is used for design vector x=(x ₁, x ₂, x ₃, x ₄) ^toptimizing, realize by step 6; Internal layer optimization is then for accounting temperature field T _j(α ⁱ, x) about block design parameter alpha ⁱthe bound of responding range, by step 5 Equivalent realization;

Step 5: the rapid solving of temperature field responding range, except optimization method, temperature field is about block design parameter alpha ⁱresponding range can be tried to achieve by interval numerical evaluation, first set up the limited configurations volume-based model of thermal protection system, adopt six point symmetry discrete scheme, can obtain as follows containing the transient state temperature field limited bulk equation of interval parameter:

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

Wherein T ^krepresent all Nodes temperature vectors of kth time step;

The matrix of coefficients of above-mentioned equation and right-hand-side vector are carried out Taylor expansion at parameter nominal value place can obtain:

$A\left({\mathrm{\&alpha;}}^I\right)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\left({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c\right)=A\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&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}{\&Sigma;}}\frac{\&part;B}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}\left({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c\right)=B\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\&part;B}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&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}{\&Sigma;}}\frac{\&part;F}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^c)=F\left({\mathrm{\&alpha;}}^c\right)+\underset{i=1}{\overset{n}{\&Sigma;}}\frac{\&part;F}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\mathrm{\&Delta;\&alpha;}}_i{\mathrm{\&delta;}}_i=F^c+{\mathrm{\&Delta;F}}^I$

Utilize Neumann progression known further:

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

Wherein $A_i=\frac{\&part;A}{\&part;{\mathrm{\&alpha;}}_i}{|}_{{\mathrm{\&alpha;}}^c}{\left(A^c\right)}^{-1};$

When the symbolic variable r in sum formula gets different value, in above formula item can specifically launch, and can obtain after merging like terms:

If norm condition || △ α _ia _i|| <1 sets up, then level is several be convergence, therefore, if 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{n}{\&Sigma;}}\underset{r=1}{\overset{\mathrm{\&infin;}}{\&Sigma;}}{(-{\mathrm{\&Delta;\&alpha;}}_i^I{A}_i)}^r={\left(A^c\right)}^{-1}+{\left(A^c\right)}^{-1}\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^I$

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

Be 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, utilize the fundamental operation rule of intervl mathematics, can obtain:

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

$\begin{array}{c}\mathrm{\&Delta;}{\left(T^{k+1}\right)}^I={\left(A^c\right)}^{-1}\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&Delta;E}}_i^I\left(B^c{\left(T^k\right)}^c+F^c\right)+\left(I+\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^c\right)\left(B^c\mathrm{\&Delta;}{\left(T^k\right)}^I+{\mathrm{\&Delta;B}}^I{\left(T^k\right)}^c+{\mathrm{\&Delta;F}}^I\right)\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&Delta;E}}_i\&CenterDot;\mathrm{\&delta;}\&CenterDot;\left(B^c{\left(T^k\right)}^c+F^c\right)+\left(I+\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^c\right)\left(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;}\right)\&rsqb;\\ ={\left(A^c\right)}^{-1}\&lsqb;\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&Delta;E}}_i\left(B^c{\left(T^k\right)}^c+F^c\right)+\left(I+\underset{i=1}{\overset{n}{\&Sigma;}}{E}_i^c\right)\left(B^c\mathrm{\&Delta;}\left(T^k\right)+\mathrm{\&Delta;}B{\left(T^k\right)}^c+\mathrm{\&Delta;}F\right)\&rsqb;\&CenterDot;\mathrm{\&delta;}\\ =\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\mathrm{\&delta;}=\mathrm{\&Delta;}\left(T^{k+1}\right)\&CenterDot;\&lsqb;-1,1\&rsqb;\end{array}$

Wherein:

$\mathrm{\&Delta;}\left(T^{k+1}\right)=|\underset{i=1}{\overset{n}{\&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}{\&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:

$\begin{array}{cc}{\stackrel{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c+\mathrm{\&Delta;}\left(T^{k+1}\right)& {\underset{\&OverBar;}{T}}^{k+1}={\left(T^{k+1}\right)}^c-\mathrm{\&Delta;}\left(T^{k+1}\right)\end{array}$

Utilize the interval Finite Volume Method proposed in this step, thermal protection structure temperature field responding range can be determined fast, thus instead of the internal layer optimization in nested Optimized model described in step 4 with interval arithmetic, improve counting yield;

Step 6: solving of deterministic optimization problem, according to the interval Finite Volume Method in step 5 to the rapid solving of structure transient state temperature field, the former nested optimization problem containing interval parameter is converted into individual layer deterministic optimization problem; Adopt simulated annealing, write calculation procedure, definition maximum cycle Iter _maxwith converging factor ε, when any one in following three conditions is met, calculates and stop:

(1) loop iteration frequency n >Iter _max;

(2) in double iterative process, objective function nominal value relative variation meets:

$\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), the initial value that given design variable is new, and be brought in algorithm and recalculate; When algorithm stops because of condition (2) or (3), get the result of calculation x of i-th iterative process ⁽ⁱ⁾as the optimal value of design variable, complete the dimensionally-optimised design process of thermal protection system based on section reliability, to reach optimum weight loss effect.

2. a kind of size optimization design method for aircraft thermal protection system based on section reliability according to claim 1, is characterized in that: main structure serviceability temperature scope required in described step 2 is no more than 150 DEG C, i.e. T ^max=150 DEG C; Engineering requirements design variable meets 2mm≤x ₁, x ₂, x ₄≤ 15mm10mm≤x ₃≤ 40mm.