- Email: [email protected]
Improvement of mixed time implicit-explicit algorithms for thermal analysis of structures
COMPUTER METHODS IN APPLIED IvIECI-IANICS AND ENGINEERING NORTH-HOLLAND PUBLISHING COMPANY
37 (1983) 207-223
IMPROVEMENT OF MIXED TIME IMPLICIT-EXPLICIT ALGORITHMS FOR THERMAL ANALYSIS OF STRUCTURES Wing Kam LIU Northwestern
University,
Department of Mechanical and Nuclear Engineering, Evanston, IL 60201, U.S.A.
and Yi Fei ZHANG Tianjin Institute of Light Industry,
Tianjin, China
Received 23 July 1982 Revised manuscript received 25 October 1982
Computer implementation aspects and numerical evaluation of the recently introduced mixed time implicit-explicit algorithms in thermal analysis of structures are presented. A computationally useful method of estimating the critical time step for a linear quadrilateral element is given herein for the methods introduced by Liu and co-workers. Numerical tests confirm the stability criterion and accuracy characteristics of the methods. The superiority of these mixed time methods to the fully implicit method or the fully explicit method is also demonstrated.
1. Introduction The demand for an optimal design of high temperature structure requires effective coupled thermal and stress analysis techniques. Finite element methods facilitate the modeling of such complicated problems. However, the resulting semi-discrete equations may involve many thousand degrees of freedom. When the problem to be solved is transient and nonlinear, due to the different time scales of the thermal response, the selection of a time integration ;nethod can be a difficult yet critical factor in the efficient solution of such problems. For this reason Adelman and Hafka [l] have conducted a survey study on the performance of explicit and implicit algorithms for transient thermal analysis of structures. Based upon their studies they concluded that, generally, implicit algorithms are preferable to explicit algorithms for ‘stiff’ problems, though non-convergence wide-banding and/or demanding computer storage of the resulting matrix equations may decrease the advantage of the implicit methods. To circumvent these difficulties methods have been developed in which it is attempted to simultaneously achieve the attributes of both classes of algorithms [2,3]. Recently, Liu and co-workers [4-6] developed a family of mixed time integration schemes in which different time integration methods (implicit/explicit) with different time steps can be used in each element group. The stability characteristic of these mixed time partition algorithms has been studied in [5]. The energy balance technique is employed to carry out the stability analysis. It is found
00457825/83/0000-0000/$03.00
@ 1983 North-Holland
20X
W. K. Liu,
Y.F. Zhang,
Mixed
implicit-explicit
ulgorithms
that the critical time-step restrictions are governed by each individual explicit element groul) only. The one-dimensiona computer impIement~tion aspects of these proposed techni~lues arc investigated in 143. A number of simple one-dimensional numerical examples are presented in ]4-51 to demonstrate the feasibility (accuracy and stability behavior) of the algorithms. In the present paper, we extend these mixed time implicit-explicit concepts to transient analysis of large scale multi-dimensional structures. Emphasis is placed on the effective computer implementation aspects, stability, and improvement of these recently introduced partition procedures. The ‘active column equation solver’ (see [3-h] for a discussion) is the key to success of these techniques. The equations system for each element group is constructed in such a way that the groups are uncoupled and hence each group can be integrated at its own group time step with different integration methods. The active column equation solver used in 13-61 is modified to enhance the efficiency of the mixed time implicit-explicit algorithms. Of more significance is that the proposed methods provided a natural framework for the further development of efficient, clean, and modularized computer codes. In Section 2 we outline the computer implementation aspects of these mixed time methods. In Section 3 the critical time-step estimates for these finite element solutions to twodimensional heat-conduction transients are given. A simple computationally useful formula is also given for estimating the critical time-step calculation. In Section 3 an illustrative example problem is described to demonstrate the practicability and usefulness of our proposed approach. In addition, the selection of a time integration method and the selection of an element group time step for each group are illustrated. Three numerical example solutions are presented in Section 5 to evaiuate the performance (i.e., accuracy and stability behavior, computer storage, and solution time, etc.) of these mixed time finite element al~~~rithnls. This represents the first comprehensive study of the effectivenesss of the proposed methods. Related discussion and conclusions are presented in Section A.
2. Computer
impIementatio~
of mixed time implicit~xpIicit
algorithms
The purpose of this section is to introduce a more versatile mixed time integration method for transient analysis of large scale multi-dimensional structures based on the previous developments by Liu [4,5]. The one-dimensional computer implementation of these mixed time methods has been presented in [4]. Much of the present formulation in this section follows the development of [4]. Also, the notations used in this paper are the same as those in [4]. As a result a preliminary reading of [4] is helpful for a better understanding. To facilitate the presentation, we simply write down the following definitions [4]: TM: master time which is incremented by At; TNEG: element group time which is incremented by AtNecj: TNoDE: nodal time which is incremented by AtNonF.; TNEO: equation time which is incremented by hfNEO. With these definitions the mixed time implicit-explicit algorithms are to proceed over the time interval [O, T,,,]. The procedures are as follows. SE/, 1. Initialization. VO.
Set TM, 7k3.
7NOr)E and TNFLo= 0, and define
the initial data f& and
W. K. Liu, Y.F. Zhang, Mixed implicit-explicit
Step 2. Form and factorize K*
=
K* where NEMEL
NUMEG c NEG=
209
algorithms
K*NEG,
K’NEG
=
c
K*‘,
e=l
1
M’ + ffhtNE,$
if implicit ,
M”
if explicit
K*e = i
(lumped
capitance
matrix
is assumed).
Step 3. Time step loop. Step 4. TMtlTM+ At, if TM 2 TM- stop. Step 5. Set F* = a AtNEoFDISCRETE. Step 6. Loop on element groups NEG = 1, . . . , NUMEG. go to Step 6b. If TNEG+ A tNEG- TM>Torr’ Step 7. Loop on elements e = 1, . . . , NUMEL; Step 7a. Define predictor values &. If TeNoDE+ At&,DE(1 - ff )A t&ODE &ODEOtherwise 0” = &&DE + (1 - a)( TM - T&DE) V&D,. Step 7b. Form element effective force f*“.
Me&
TM
then
& = 8&,DES
if implicit ,
f*e=
I Me&
- a
where
WKe&
if explicit ,
W = diagonal matrix with At&,DE along the diagonals. 7c. Sum up effective force from elemental contributions F* + F* + f*e. 7d. End of element loop. 6a. Update TNEG. TNEG = TNEG+ AtNEG. 6b. End of element group loop. 8. Solve for 8. (In order to enhance the efficiency of these mixed time methods the active column equation solver used in [3-61 is modified in the forward reduction and backsubstitution procedures as follows.) Step 8a. Loop on equation number N = 1, . . . , NEQ. If T:EEQ+ At&, - TM > Torr go to Step 8b. Forward reduction and backsubstitution for equation N. Step 8b. End of equation number loop. Step 9. Update V and 8. Step 9a. Loop on N = 1,. . . , NUMNP. If T&E+ AtKoDE- TM > Torr go to Step 9b. eNcsolution from Step 8. V” t (8” - h”)/aAt&,,. Step 9b. End of nodal number loop. Step 10. Update TNODE. step 10a. Loop on N = 1, . . . , NUMNP. If T&DE f AtF;ODE - TM > Torr go to Step lob.
Step Step Step Step Step
1 + means ‘is replaced by’. ’ We use Torr = 1.0 E - 9.
W.K. Liu, Y.F. Zhung, Mixed implicit-explicit
210
algorithms
+Kam + Wkm~. Step lob. End of nodal number loop. Step 11. Update TNEU. Step lla. Loop on N = I,. . . , NEQ. If T:EEO+ A&, - TM > Torr go to Step 1lb. T&o +- 7’KEIEQ + A&o. Step llb. End of equation number loop. Step 3a. End of time step loop.
T&mE
3. Critical time-step estimates for bilinear quadrilateral
element
A detailed
study of the eigenvalues of the heat conduction bilinear quadrifateraf element based on one-point quadrature is given in 171. A stabilization procedure is also developed in [7] for this one-point quadrature element because of the additional singular mode (often called hourglass mode) other than the singular mode associated with the constant temperature field. For the one-point quadrature element there are four eigenvalues of which only two of them are nonzero. One of these zero eigenvalues can be removed if the stabilization stiffness is added to the one-point quadrature stiffness. This third eigenvalue can be shown to be bounded between the other two nonzero eigenvalues. In particular, for the case of rectangular element these three eigenvalues are identical to those three obtained using a two-by-two quadrature element. Hence a computationally useful method of estimating the critical time step for linear quadrilateral element can easily be derived for our mixed time methods. and w,,, is the maximum Following f4] AtGritG 2/wmax is the critical time-step restriction eigenvafue of the eth element subjected to insulated boundary conditions. For an arbitrary four-node quadrilateral as shown in Fig. l(a), let Xt
=
($1,
X27
-x3,
and
x4)
are the x and y-coordinates
yf
=
(YE,
YZ, ~3, ~4)
of the four nodes.
Fig. i(a). Arbitrary
four-node
Define
quadrilateral
element:
(b) rectangular
element.
WK. Liu, Y.F. Zhang, Mixed implicit-exp~~cjt algorithms XIJ = XI -
XJ
and
YIJ
=
YI
-
YJ
for&J=1
,...,
211
4.
Let Y = y& + $3, and A is the area of the quadrilateral. shown that Atcrit d mm /L
2 = -(yz4y31+ X&3),
Combining
the results given in [4] and [7] it can be
2A2 a((X + Y) r V/(x - Y)’ + 422}
where a is the thermai diffusivity and Al, is between 0.95 and 1.0. We have performed a numerical study of the eigenvalues of the various shapes (for all practical purposes) of two-by-two quadrature elements. It is found empirically that p can be picked to be 0.95 if the element is really skewed, and g can be picked to be 1.0 if the element is rectangular. In the special case of a rectangular element (i.e. with or.= 1.0)
where 1, and 1, are defined in Fig. l(b).
4. Illustration of the selections of a time integration method and an element group time step We have formulated herein a two-dimensional finite element model for the thermal protection system (TPS) of the space shuttle [f] to illustrate the selection procedures for a time integration method and an element group time step. The practicability and usefulness of the algorithms can also be demonstrated with this model as will be described in Section 5. Since the emphasis in this paper is on the evaluation of the methodologies described in previous sections, nonlinear effects such as internal and external radiation are ignored. Further only the assumed mean temperature material properties of the various components of the TPS are considered. The finite element mesh is depicted in Fig. 7. Due to symmetry only half of the TPS is modeled. Its material properties are tabulated in Table 1. The number of elements, the minimum characteristic length Imingthe estimated explicit critical time step Atmi,, the proposed integration method, and the proposed element group time step for each group are included in Table 2. As can be seen from Table 2 due to the various thermal time scales (e.g. Atmi, = 0.041 for Al, and At,, = 16.71 for RSI), a single integration method is definitely not effective. For example, if an explicit method is employed, a time step of 0.041 has to be used; while if an implicit method is emptoyed, there is no stability-imposed limitation on At; however, widebanding and/or demanding computer storage of the resulting matrix equations may decrease its advantage. The family of mixed time integration schemes developed is best suited for this type of problem. The attributes of the various time integration methods are fully achieved using the proposed approach as can readily be seen from Table 2. It should also be observed that Atmincan be set as high as 16.71 in this mixed model. Subsequently At = 16.71 is employed for element groups 1 and 2 (042 coating and RSI), At = 66.84 for element group 3 (RSI), and At = 66.84 for element groups 4 and 5 (RTV, felt, Al and air). However, judging from the heat
212
W. K. Liu, Y.F. Zhang,
Table 1 Mean temperature
material
Mixed
implicit-explicit
algorithms
properties Material 042 coating
Properties
RSI
RTV
Felt
AI
Air ____
A w/m-k conductivity
1.4338
0.1354
0.3113
0.0363
135373
0.0271
density C J/kg-k
1665.92
144.17
1409.62
96.11
285 1.28
1.126
specific heat a = h/pCm’/sec
1317.96
1255.20
1213.36
1213.26
962.32
1009.0
6.52 x lo-’
7.48 x lo-’
1.82x lo- ’
3.11 x lo-’
4.78 x 10-S
2.39 x IO ’
P
kg/m3
thermal
diffusivity
The mean temperature 333 k and 300 k.
Table 2 Characteristic
length
material
4 5 I = implicit
are computed
from the average
of those at 1200 k. 950 k. 500 k, 477 k.
and time scales
Element group number
properties
Material type
Number of elements (total: 120)
Estimated I,i,
(cm)
ALi”
Integration method’
Element group time step
042 coating RSI RSI
12
0.04
0.1225
I
16.7 I (=At)
30 18
0.5
1.o
16.71 66.84
E E
16.71 66.84 (=4AT)
RTV Felt Al Air
18 6 16 20
0.02 0.40 0.20 2.0
0.106 25.70 0.041 8.36
I
66.84
I
66.84
integration;
E = explicit
integration.
loading versus time curve (see Fig. 3) it is advisable to pick At = 25.7/4 = 6.425 during the heat loading period (i.e. time = 0 to 835.5) and At = 16.71 in later time. Of course, an automatic selection of the step size based on truncation error and accuracy considerations is desirable. However, these latter two are not being considered here. The importance of this proposed mixed time implicit-explicit finite element concepts can further be visualized if a nonlinear analysis of the above finite element model is assumed. The thermal responses of the various components of the TPS may be divided into regions of slowly and rapidly varying temperatures as follows. (1) 042 coating. The thermal diffusivity is almost independent of temperature and the estimate At is so small (0.1225) that explicit calculation is not cost effective at all. Implicit
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms
213
calculation is best suited since only limited partial reformation and refactorization of this element group’s effective stiffness will be needed. This is one of the major advantages of the mixed time implicit-explicit concept. (2) RSI. The thermal diffusivity changes rapidly with temperature, therefore reformation and refactorization of this element group’s effective stiffness are frequently required if the implicit method is employed. Nonconvergence and/or wide-banding of the resulting element group effective stiffness may decrease the advantage of applying the implicit method. The estimated At (based on mean temperature value) is 16.71; therefore the explicit method is best suited (no matrices operations). (3) RTV. Its thermal diffusivity is independent of temperature and the estimated At is so small (0.106) that implicit method is recommended. (4) AZ. Its thermal diffusivity is fairly independent of temperature (at least in the operational temperature range of the space shuttle), and its estimated At is so stringent (0.041) that implicit method is again recommended. (5) Felt. Its thermal diffusivity is independent of temperature and its estimated At is large (25.7) therefore the implicit method or explicit method can be used. For this example we propose to employ the implicit method. However, for 3D calculation in order to reduce the bandwidth and the unnecessary nonlinear calculations (due to the fact that the adjacent groups might be implicit groups too), the explicit method will be recommended. (6) Air. The variations of the thermal diffusivity with temperature are small in our range of interest. Also the estimated At is small, therefore the implicit method is again well suited. REMARKS.
(1)
In the above At calculations
lmin= min{l,, I,} where 1, and 1, are defined in
Section 3. (2) The estimated At is defined to be 1$iJ2a. It is a safe critical time-step calculation. See [4,5] for a discussion. (3) By virtue of the fact that each element group’s effective stiffness is uncoupled to the global assembled matrix equations system any element group can be reformed and refactorized at any instant if required without affecting the global equations system. This partial factorization procedure can further be enhanced if it is to be combined with an iterative update procedure. (4) In the case of short time duration transients mixed time explicit-explicit is most effective. It is even more attractive if the mixed time explicit-explicit procedures are unconditionally stable [8].
5. Numerical results
A two-dimensional finite element pilot computer code incorporating the methodologies described in previous sections has been written to evaluate the performance of these mixed time finite element algorithms. Three numerical examples are presented to demonstrate the accuracy, stability, and efficiency of these proposed methods. All computations are performed on a CDC Cyber 170/730 computer in single precision (60 bits per floating point word), and lumped capacitance matrices are used throughout.
214
W.K. Liu, Y. F. Zhang, Mixed implicit-explicit
algorithms
EXAMPLE 5.1.Comparison with finite diflerence solution. The finite element mesh as shown in Fig. 2 consists of 100 uniform square elements (each with 1, = I, = 10). The mesh is deliberately divided into four element groups. (1) Group 1 is integrated implicitly with a group time step of 2At. (2) Group 2 is integrated explicitly with a group time step of 2At. (3) Group 3 is integrated implicitly with a group time step of 4At. (4) Group 4 is integrated explicitly with a group time step of At. We called the above partition a 21-2E-41-1E partition. Observe that group 3 (41) is chosen to be at the center of the square. Its purpose is to confirm the generality of our modified column equation solver. As described in [3-61 the variable column heights (which affect the solution time) depend on the nodal numbering. Furthermore, no forward reduction and backsubstitution are performed until every 4At time interval. As can be pictured from Fig. 2 the column heights of the ‘4At’ equations are coupled to the ‘At’ and ‘2At‘ equations, therefore unlike the one-dimensional case (see [4] for a discussion) a special forward reduction and backsubstitution procedure as described in Section 2 is required for efficiency. The thermal diffusivity a is chosen to be 6.25 for all the elements. Node 1 is kept at a constant temperature of 0.0 while all the other nodes are subjected to a constant initial temperature of 0.1. All four sides are insulated. A time step of At equal to 1 is employed for the transient analysis. The temperature-time histories of nodes 3, 46, 57, and 67 are depicted in Fig. 3. In order to confirm the accuracy of these solutions the mixed time finite element solutions are compared to those obtained using a finite difference 21 X 21 (400 zones) grid with a single time step of At equal to 1. The finite difference solutions of these four nodes are also depicted in Fig. 3. As can be seen the comparison is quite good despite the crudeness of the mesh (as compared to the finite difference method) and the larger time steps used in the finite The to 1, 2, and 4 for groups 4, 1 and 2, and 3 respectively). element mesh (1 as compared maximum relatively difference between the two solutions is no more than 2%.
Fig. 2. Problem
EXAMPLE
statement
and finite element
mesh of numerical
Example
5.1.
5.2. Demonstration of accuracy, stability, and ejJiciency. The problem statement is depicted in Fig. 4. The finite element mesh consists of 400 uniform size rectangular elements (each has 1, = 5 and 1, = 10) and 486 nodes (81 in x-direction and 6 in y-direction) where the x-axis is defined by joining node 3 to node 17. The material properties of this finite element
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit
0
100
800
TEE
7.ooo 4
0
Fig. 3. Comparison
215
algorithms
100
lEl3
between finite element and finite difference solutions.
TIME
Fig. 4. Problem statement
and finite element mesh of numerical Example 5.2.
800
216
W.K. Liu,
Y.F. Zhang.
Mixed
~~~fj~it-~~pficit
u~gorit~~.s
model are composed of four different thermal diffusivities of O.l2Sa, 0.25~. a and 40~ respectively as shown in Fig. 4. a is chosen to be 12.5. In order to simulate a realistic three-dimensional bandwidth using this two-dimensi~~nal finite element model we number the nodes along the x-axis. The heat load which is also shown in Fig. 4 is applied at node I. The initial temperature for all the nodes is 0.1. All four sides are insulated. AS mentioned earlier, due to the different thermal di~usivities, a single time integration method is inefficient for this problem. For example, if the explicit method is employed, a time step of 0.025 is required so that 16000 steps (as compared to 400 steps) are required for this problem. Therefore the purely explicit integration method is not being considered here. If the implicit method is employed, we can employ a much larger time step though the wide bandwidth and demanding computer storage (e.g. 3288 vs. 33771 single precision words) may decrease the advantage of the implicit method. This can readily be seen in Table 4. Hence the mixed time implicit-explicit method is best suited for this exampfe. A totaf of five computer runs are made to demonstrate the accuracy, stability, and efficiency of the proposed mixed time methods and they are: (1) lE-4E-81-8E with At = I; (2) 11 with At = 1; (3) 11-lE-21-2E with At = 4; (4) 11 with At = 4; (5) 11 with At = 8. In run 1, groups 1, 2 and 4 are integrated expiicitly with lE, 4E and 8E respectively; and group 3 is integrated implicitly. In run 3, groups 2 and 4 are integrated explicitly with 1E and 2E respectively; while groups 1 and 3 are integrated implicitly with 11 and 21 respectively. The time step At and the element group time steps are chosen according to the critical time-step formula given in Section 3 with Al;equal to 1.0. The temperature-time histories at six different locations for run number 1, i.e. lE-4E-SI-8E, are presented in Fig. 5. These results are virtually identical to those obtained from run number 2, i.e. 11. It should be pointed out that the curves shown in Fig. 5 are plotted for every time interval of 4 except for run number 5 in which the time interval is 8. Since the plot of node 327 is integrated with a time step of 8, therefore the temperature-time history shown is ‘step-wise’ (from time -4X to 240). The solution times, and storage requirements for the above five runs are tabulated in Tables 3 and Table 3 Solution times (in seconds) Method
Item Formulation of stiffness and load Factorization Forward reduction and backsubstitution Total solution time
II (At = 1)
1I- IE-2E-21
78.96 (1.34
21 I.13 11.21
58.13 0.80
52.98 11.55
28.81
10.95 90.25
293.58 515.92
10.90 69.83
74.24 138.77
37.50 77.84
1E-4E-XE-81 (At = 1)
(At = 4)
(AZ
4)
I t.s?l
217
a+
TIME
TIME
TIME
.Tlm
.ooil
TIME Fig. 5. Temperature-time
histories of iE-4E-8I-8E
mixed time method with At = 1.
218
W. K. Liu, Y.F. Zhang,
Mixed implicit-explicit
16
-lI-lE-2I-2E AND 4I(DT=4)
12--
+
r
62
T-
167
_______ ..
lI(DT=8)
4
66
ulgorithms
1
3
2
1
0
3.000
327
1.300
l.ooo
.500
4
KuI
4 0.000t 0
Fig. 6. Comparison of solutions of five integration methods.
!
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit Table 4 Storage (in single precision
219
algorithms
words) Method
Item Number of equations (NEQ) Number of terms in stiffness (NA) Average half bandwidth
lE-4E-8E-81 (At = 1)
(At’: 1)
lI-lE-2E-21 (At = 4)
(At’! 4)
(At’: 8)
486
486
486
486
486
3288
33771
6089
33771
33771
6
69
12
69
69
MB zNA,NEQ
4 respectively. From these tables the advantages of these mixed time integration methods are apparent. The relative accuracy of these five runs are compared and they are shown in Fig. 6. These curves are plotted at every time interval of 8. In summary, from accuracy consideration, runs 1 and 3 are compatible to runs 2 and 4 respectively while run 5 is not as accurate as the others. However, judging from the solution times and storage requirements (e.g. 90.25 vs. 515.92 and 3288 vs. 33771 respectively), the superiority of these mixed time methods to the implicit method is fully illustrated. The potential savings in both computer time and computer storage are apparent in three-dimensional and/or nonlinear calculations. Finally, the critical time-step estimate is also confirmed by this numerical example. EXAMPLE
5.3. Applications to TPS. The space shuttle orbiter is designed to be used for at least 100 missions. Currently, various metallic TPS concepts are being investigated at the NASA Langley Research Center (see e.g. [9, lo] for a discussion). One of the proposed candidates is the titanium multiwall tile. A titanium multiwall tile consists of alternating layers of superplastically formed dimpled sheets and flat septum sheets of titanium foil. As described in [lo] this multiwall concept impedes all three modes of heat transfer--conduction, radiation, and convection. The superplastically formed dimpled sheets and the long thin conduction path tend to minimize heat conduction. The flat septum sheets of titanium foil impede radiation. The small individual volumes created by the dimpled layers virtually eliminate air convection. We employed two integration methods to compute the time histories of the top (node 84) and bottom (node 64) surface temperatures of the 042 coatings and the lower (nodes 11 and 137) and upper (nodes 15 and 57) aluminum skin temperatures for the two-dimensional finite element model described in Section 4. Its dimensions (not to scale and in centimeters) are depicted in Fig. 7. The heat loads which are also shown in Fig. 7 are applied at the top and bottom of the TPS model. The initial temperature for all the nodes is -20. Sides AD and BC are insulated.
220
Fig. 7. Problem
W.K. Liu, Y.F. Zhang,
statement
and finite element
Mixed implicit-explicit
mesh of a two-dimensional
ulgorithrns
thermal
protection
systems
(TPS).
As described in Section 4 407317 steps (as compared to 1000 steps) are required if the purely explicit method is employed. Therefore, we do not consider the purely explicit time integration method. We proposed to use lI-lE-4E-41-41 time integration with a time step of 16.71 for the 042 coating and part of RSI, and a time step of 66.84 for the other part of RSI, RTV, felt, Al and air. The computed results are then compared to those using purely implicit method with a single time step of 16.71. The temperature-time histories (plotted at every time interval of 66.84) at six different locations (042 coatings and aluminum skin) are presented in Fig. 8. These results are virtually identical to those obtained from the purely implicit case, i.e. 11. The solution times and computer storage requirements for these two runs are tabulated in Tables 5 and 6 respectively. Even though we use only a total of 147 nodes and 120 elements (a very small mesh as compared to a typical 3D finite element model of 3000 nodes and 4000 elements) we gain a factor of about 1.7 in solution time and a factor of about 1.5 in computer storage. Again, the advantages as well as the stability of these mixed time methods are fully ilfustrated.
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms loo1X-lE-4E-41-41
85--
30-m
-S--
_I
-40-r 0.000
LOWER
-5-w
AL
NODE 137 4.000
8.000
TIME
LOWER
AL
J
12.000
16.000
20.000
I,
-40 0.000
4.000
(Xl 0”)
8.000
TIME
12.000
16.000
20.000
(Xl 03)
00
40
20
0
UPPER
AL
-20
UPPER
AL
4.000
8.000
If 0
4.000
8.000
TIME
12.000
16.000
21
J,
10
3
(Xl 03)
TIME
12.000
16.000
20 00
(Xl 03)
3ooo-
I
I 170--
8840~-
lOo--
1480~-
LOWER
042
UPPER so-
720--
64
-40-c 0.000
04’2
4.000
8.000
TIME
12.000
16.000
2o.oOa
,
o.ooo 4.000
(Xl 03)
Fig. 8. Comparison
84
-40 8.000
TIME of solutions
of two integration
methods.
12000
(Xl 0’)
l&m0
20
222
WK. Liu, Y.F. Zhang, Mixed ~rnp~ic~t-explicit ~lgo~~thrns
Table 5 Solution times (in seconds) Method II-lE-4E-41-41 (At = 16.71)
Item Formulation
(At =1:6.71)
of
stiffness and load Factorization Forward reduction and backsubstitution Total solution time
Table 6 Storage (in single precision
105.12 0.36
162.54 0.57
45.05 160.42
91.00 263.35
words) Method
Item Number of equations (NEQ) Number of terms in stiffness (NA) Average half bandwidth MB “=”NA/NEQ
II-lE-4E-41-41 (At = 16.71)
(Al =1:6.71)
147
147
2121
2933
14
19
6. Conclusions The computer implementation aspects and the evaluation of the performance of mixed time implicit-explicit finite element as applied to multi-dimensional transient thermal analysis are presented. The paper focuses on applications where the thermal model would normally have different time scales of thermal response. Three two-dimensional examples are presented to ihustrate the approach: (1) comparison with finite difference method, (2) demonstration of accuracy, stability, and efficiency, and (3) application to thermal protection systems. The examples show: (1) the superiority of mixed time methods to a single integration method (either implicit or explicit), (2) potential savings in computer time and computer storage, (3) the accuracy and stability behavior of mixed time finite element, and (4) the importance of these mixed time methods as applied to three-dimensional and/or nonlinear thermal analysis.
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms
223
Acknowledgment
We wish to thank Professor T. Belytschko for his many helpful discussions during the course of the present study. We also wish to express our gratitude to the computer center at Northwestern University for providing computer services. The support of NASA under Grant No. NAG-l-210 to this research is gratefully acknowledged.
References [I] H. Adelman and R. Haftka, On the performance of explicit and implicit algorithms for transient thermal analysis of structures, NASA TM 81880, 1980. [2] T. Belytschko and R. Mullen, Stability of explicit-implicit mesh partitions in time integration, Internat. J. Numer. Meths. Engrg. 12 (1978) 1575-1586. [3] T. Hughes and W. Liu, Implicit-explicit finite elements in transient analysis, J. Appl. Mech. 45 (1978) 371-378. [4] W. Liu, Development of mixed time partition procedures for thermal analysis of structures, Internat. J. Numer. Meths. Engrg. (1983) to appear. [5] W. Liu and J. Lin, Stability of mixed time integration schemes for transient thermal analysis, Numer. Heat Transfer J. 5 (1982) 445450. [6] W. Liu and T. Belytschko, Mixed time implicit-explicit finite elements for transient analysis, Comput. & Structures 15 (1982) to appear. [7] W. Liu andT. Belyschko, Efficient linear and nonlinear heat conduction with a quadrilateral element, Internat. J. Numer. Meths. Engrg. to appear. [8] W. Liu and Y. Zhang, Unconditionally stable implicit-explicit algorithms for coupled thermal stress waves, Comput. & Structures (1983) to appear. [9] H. Kelly, D. Rummler and R. Jackson, Research in structures and material for future space transportation systems-an overview, AIAA Paper 79-0859, AIAA Conf. Advanced Technology for Future Space Systems, 1979. [lo] R. Jackson and S. Dixon, A design assessment of multiwall, metallic stand-off and RSI reusable thermal protection system including space shuttle applications, NASA TM 81780, 1980.
37 (1983) 207-223
IMPROVEMENT OF MIXED TIME IMPLICIT-EXPLICIT ALGORITHMS FOR THERMAL ANALYSIS OF STRUCTURES Wing Kam LIU Northwestern
University,
Department of Mechanical and Nuclear Engineering, Evanston, IL 60201, U.S.A.
and Yi Fei ZHANG Tianjin Institute of Light Industry,
Tianjin, China
Received 23 July 1982 Revised manuscript received 25 October 1982
Computer implementation aspects and numerical evaluation of the recently introduced mixed time implicit-explicit algorithms in thermal analysis of structures are presented. A computationally useful method of estimating the critical time step for a linear quadrilateral element is given herein for the methods introduced by Liu and co-workers. Numerical tests confirm the stability criterion and accuracy characteristics of the methods. The superiority of these mixed time methods to the fully implicit method or the fully explicit method is also demonstrated.
1. Introduction The demand for an optimal design of high temperature structure requires effective coupled thermal and stress analysis techniques. Finite element methods facilitate the modeling of such complicated problems. However, the resulting semi-discrete equations may involve many thousand degrees of freedom. When the problem to be solved is transient and nonlinear, due to the different time scales of the thermal response, the selection of a time integration ;nethod can be a difficult yet critical factor in the efficient solution of such problems. For this reason Adelman and Hafka [l] have conducted a survey study on the performance of explicit and implicit algorithms for transient thermal analysis of structures. Based upon their studies they concluded that, generally, implicit algorithms are preferable to explicit algorithms for ‘stiff’ problems, though non-convergence wide-banding and/or demanding computer storage of the resulting matrix equations may decrease the advantage of the implicit methods. To circumvent these difficulties methods have been developed in which it is attempted to simultaneously achieve the attributes of both classes of algorithms [2,3]. Recently, Liu and co-workers [4-6] developed a family of mixed time integration schemes in which different time integration methods (implicit/explicit) with different time steps can be used in each element group. The stability characteristic of these mixed time partition algorithms has been studied in [5]. The energy balance technique is employed to carry out the stability analysis. It is found
00457825/83/0000-0000/$03.00
@ 1983 North-Holland
20X
W. K. Liu,
Y.F. Zhang,
Mixed
implicit-explicit
ulgorithms
that the critical time-step restrictions are governed by each individual explicit element groul) only. The one-dimensiona computer impIement~tion aspects of these proposed techni~lues arc investigated in 143. A number of simple one-dimensional numerical examples are presented in ]4-51 to demonstrate the feasibility (accuracy and stability behavior) of the algorithms. In the present paper, we extend these mixed time implicit-explicit concepts to transient analysis of large scale multi-dimensional structures. Emphasis is placed on the effective computer implementation aspects, stability, and improvement of these recently introduced partition procedures. The ‘active column equation solver’ (see [3-h] for a discussion) is the key to success of these techniques. The equations system for each element group is constructed in such a way that the groups are uncoupled and hence each group can be integrated at its own group time step with different integration methods. The active column equation solver used in 13-61 is modified to enhance the efficiency of the mixed time implicit-explicit algorithms. Of more significance is that the proposed methods provided a natural framework for the further development of efficient, clean, and modularized computer codes. In Section 2 we outline the computer implementation aspects of these mixed time methods. In Section 3 the critical time-step estimates for these finite element solutions to twodimensional heat-conduction transients are given. A simple computationally useful formula is also given for estimating the critical time-step calculation. In Section 3 an illustrative example problem is described to demonstrate the practicability and usefulness of our proposed approach. In addition, the selection of a time integration method and the selection of an element group time step for each group are illustrated. Three numerical example solutions are presented in Section 5 to evaiuate the performance (i.e., accuracy and stability behavior, computer storage, and solution time, etc.) of these mixed time finite element al~~~rithnls. This represents the first comprehensive study of the effectivenesss of the proposed methods. Related discussion and conclusions are presented in Section A.
2. Computer
impIementatio~
of mixed time implicit~xpIicit
algorithms
The purpose of this section is to introduce a more versatile mixed time integration method for transient analysis of large scale multi-dimensional structures based on the previous developments by Liu [4,5]. The one-dimensional computer implementation of these mixed time methods has been presented in [4]. Much of the present formulation in this section follows the development of [4]. Also, the notations used in this paper are the same as those in [4]. As a result a preliminary reading of [4] is helpful for a better understanding. To facilitate the presentation, we simply write down the following definitions [4]: TM: master time which is incremented by At; TNEG: element group time which is incremented by AtNecj: TNoDE: nodal time which is incremented by AtNonF.; TNEO: equation time which is incremented by hfNEO. With these definitions the mixed time implicit-explicit algorithms are to proceed over the time interval [O, T,,,]. The procedures are as follows. SE/, 1. Initialization. VO.
Set TM, 7k3.
7NOr)E and TNFLo= 0, and define
the initial data f& and
W. K. Liu, Y.F. Zhang, Mixed implicit-explicit
Step 2. Form and factorize K*
=
K* where NEMEL
NUMEG c NEG=
209
algorithms
K*NEG,
K’NEG
=
c
K*‘,
e=l
1
M’ + ffhtNE,$
if implicit ,
M”
if explicit
K*e = i
(lumped
capitance
matrix
is assumed).
Step 3. Time step loop. Step 4. TMtlTM+ At, if TM 2 TM- stop. Step 5. Set F* = a AtNEoFDISCRETE. Step 6. Loop on element groups NEG = 1, . . . , NUMEG. go to Step 6b. If TNEG+ A tNEG- TM>Torr’ Step 7. Loop on elements e = 1, . . . , NUMEL; Step 7a. Define predictor values &. If TeNoDE+ At&,DE(1 - ff )A t&ODE &ODEOtherwise 0” = &&DE + (1 - a)( TM - T&DE) V&D,. Step 7b. Form element effective force f*“.
Me&
TM
then
& = 8&,DES
if implicit ,
f*e=
I Me&
- a
where
WKe&
if explicit ,
W = diagonal matrix with At&,DE along the diagonals. 7c. Sum up effective force from elemental contributions F* + F* + f*e. 7d. End of element loop. 6a. Update TNEG. TNEG = TNEG+ AtNEG. 6b. End of element group loop. 8. Solve for 8. (In order to enhance the efficiency of these mixed time methods the active column equation solver used in [3-61 is modified in the forward reduction and backsubstitution procedures as follows.) Step 8a. Loop on equation number N = 1, . . . , NEQ. If T:EEQ+ At&, - TM > Torr go to Step 8b. Forward reduction and backsubstitution for equation N. Step 8b. End of equation number loop. Step 9. Update V and 8. Step 9a. Loop on N = 1,. . . , NUMNP. If T&E+ AtKoDE- TM > Torr go to Step 9b. eNcsolution from Step 8. V” t (8” - h”)/aAt&,,. Step 9b. End of nodal number loop. Step 10. Update TNODE. step 10a. Loop on N = 1, . . . , NUMNP. If T&DE f AtF;ODE - TM > Torr go to Step lob.
Step Step Step Step Step
1 + means ‘is replaced by’. ’ We use Torr = 1.0 E - 9.
W.K. Liu, Y.F. Zhung, Mixed implicit-explicit
210
algorithms
+Kam + Wkm~. Step lob. End of nodal number loop. Step 11. Update TNEU. Step lla. Loop on N = I,. . . , NEQ. If T:EEO+ A&, - TM > Torr go to Step 1lb. T&o +- 7’KEIEQ + A&o. Step llb. End of equation number loop. Step 3a. End of time step loop.
T&mE
3. Critical time-step estimates for bilinear quadrilateral
element
A detailed
study of the eigenvalues of the heat conduction bilinear quadrifateraf element based on one-point quadrature is given in 171. A stabilization procedure is also developed in [7] for this one-point quadrature element because of the additional singular mode (often called hourglass mode) other than the singular mode associated with the constant temperature field. For the one-point quadrature element there are four eigenvalues of which only two of them are nonzero. One of these zero eigenvalues can be removed if the stabilization stiffness is added to the one-point quadrature stiffness. This third eigenvalue can be shown to be bounded between the other two nonzero eigenvalues. In particular, for the case of rectangular element these three eigenvalues are identical to those three obtained using a two-by-two quadrature element. Hence a computationally useful method of estimating the critical time step for linear quadrilateral element can easily be derived for our mixed time methods. and w,,, is the maximum Following f4] AtGritG 2/wmax is the critical time-step restriction eigenvafue of the eth element subjected to insulated boundary conditions. For an arbitrary four-node quadrilateral as shown in Fig. l(a), let Xt
=
($1,
X27
-x3,
and
x4)
are the x and y-coordinates
yf
=
(YE,
YZ, ~3, ~4)
of the four nodes.
Fig. i(a). Arbitrary
four-node
Define
quadrilateral
element:
(b) rectangular
element.
WK. Liu, Y.F. Zhang, Mixed implicit-exp~~cjt algorithms XIJ = XI -
XJ
and
YIJ
=
YI
-
YJ
for&J=1
,...,
211
4.
Let Y = y& + $3, and A is the area of the quadrilateral. shown that Atcrit d mm /L
2 = -(yz4y31+ X&3),
Combining
the results given in [4] and [7] it can be
2A2 a((X + Y) r V/(x - Y)’ + 422}
where a is the thermai diffusivity and Al, is between 0.95 and 1.0. We have performed a numerical study of the eigenvalues of the various shapes (for all practical purposes) of two-by-two quadrature elements. It is found empirically that p can be picked to be 0.95 if the element is really skewed, and g can be picked to be 1.0 if the element is rectangular. In the special case of a rectangular element (i.e. with or.= 1.0)
where 1, and 1, are defined in Fig. l(b).
4. Illustration of the selections of a time integration method and an element group time step We have formulated herein a two-dimensional finite element model for the thermal protection system (TPS) of the space shuttle [f] to illustrate the selection procedures for a time integration method and an element group time step. The practicability and usefulness of the algorithms can also be demonstrated with this model as will be described in Section 5. Since the emphasis in this paper is on the evaluation of the methodologies described in previous sections, nonlinear effects such as internal and external radiation are ignored. Further only the assumed mean temperature material properties of the various components of the TPS are considered. The finite element mesh is depicted in Fig. 7. Due to symmetry only half of the TPS is modeled. Its material properties are tabulated in Table 1. The number of elements, the minimum characteristic length Imingthe estimated explicit critical time step Atmi,, the proposed integration method, and the proposed element group time step for each group are included in Table 2. As can be seen from Table 2 due to the various thermal time scales (e.g. Atmi, = 0.041 for Al, and At,, = 16.71 for RSI), a single integration method is definitely not effective. For example, if an explicit method is employed, a time step of 0.041 has to be used; while if an implicit method is emptoyed, there is no stability-imposed limitation on At; however, widebanding and/or demanding computer storage of the resulting matrix equations may decrease its advantage. The family of mixed time integration schemes developed is best suited for this type of problem. The attributes of the various time integration methods are fully achieved using the proposed approach as can readily be seen from Table 2. It should also be observed that Atmincan be set as high as 16.71 in this mixed model. Subsequently At = 16.71 is employed for element groups 1 and 2 (042 coating and RSI), At = 66.84 for element group 3 (RSI), and At = 66.84 for element groups 4 and 5 (RTV, felt, Al and air). However, judging from the heat
212
W. K. Liu, Y.F. Zhang,
Table 1 Mean temperature
material
Mixed
implicit-explicit
algorithms
properties Material 042 coating
Properties
RSI
RTV
Felt
AI
Air ____
A w/m-k conductivity
1.4338
0.1354
0.3113
0.0363
135373
0.0271
density C J/kg-k
1665.92
144.17
1409.62
96.11
285 1.28
1.126
specific heat a = h/pCm’/sec
1317.96
1255.20
1213.36
1213.26
962.32
1009.0
6.52 x lo-’
7.48 x lo-’
1.82x lo- ’
3.11 x lo-’
4.78 x 10-S
2.39 x IO ’
P
kg/m3
thermal
diffusivity
The mean temperature 333 k and 300 k.
Table 2 Characteristic
length
material
4 5 I = implicit
are computed
from the average
of those at 1200 k. 950 k. 500 k, 477 k.
and time scales
Element group number
properties
Material type
Number of elements (total: 120)
Estimated I,i,
(cm)
ALi”
Integration method’
Element group time step
042 coating RSI RSI
12
0.04
0.1225
I
16.7 I (=At)
30 18
0.5
1.o
16.71 66.84
E E
16.71 66.84 (=4AT)
RTV Felt Al Air
18 6 16 20
0.02 0.40 0.20 2.0
0.106 25.70 0.041 8.36
I
66.84
I
66.84
integration;
E = explicit
integration.
loading versus time curve (see Fig. 3) it is advisable to pick At = 25.7/4 = 6.425 during the heat loading period (i.e. time = 0 to 835.5) and At = 16.71 in later time. Of course, an automatic selection of the step size based on truncation error and accuracy considerations is desirable. However, these latter two are not being considered here. The importance of this proposed mixed time implicit-explicit finite element concepts can further be visualized if a nonlinear analysis of the above finite element model is assumed. The thermal responses of the various components of the TPS may be divided into regions of slowly and rapidly varying temperatures as follows. (1) 042 coating. The thermal diffusivity is almost independent of temperature and the estimate At is so small (0.1225) that explicit calculation is not cost effective at all. Implicit
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms
213
calculation is best suited since only limited partial reformation and refactorization of this element group’s effective stiffness will be needed. This is one of the major advantages of the mixed time implicit-explicit concept. (2) RSI. The thermal diffusivity changes rapidly with temperature, therefore reformation and refactorization of this element group’s effective stiffness are frequently required if the implicit method is employed. Nonconvergence and/or wide-banding of the resulting element group effective stiffness may decrease the advantage of applying the implicit method. The estimated At (based on mean temperature value) is 16.71; therefore the explicit method is best suited (no matrices operations). (3) RTV. Its thermal diffusivity is independent of temperature and the estimated At is so small (0.106) that implicit method is recommended. (4) AZ. Its thermal diffusivity is fairly independent of temperature (at least in the operational temperature range of the space shuttle), and its estimated At is so stringent (0.041) that implicit method is again recommended. (5) Felt. Its thermal diffusivity is independent of temperature and its estimated At is large (25.7) therefore the implicit method or explicit method can be used. For this example we propose to employ the implicit method. However, for 3D calculation in order to reduce the bandwidth and the unnecessary nonlinear calculations (due to the fact that the adjacent groups might be implicit groups too), the explicit method will be recommended. (6) Air. The variations of the thermal diffusivity with temperature are small in our range of interest. Also the estimated At is small, therefore the implicit method is again well suited. REMARKS.
(1)
In the above At calculations
lmin= min{l,, I,} where 1, and 1, are defined in
Section 3. (2) The estimated At is defined to be 1$iJ2a. It is a safe critical time-step calculation. See [4,5] for a discussion. (3) By virtue of the fact that each element group’s effective stiffness is uncoupled to the global assembled matrix equations system any element group can be reformed and refactorized at any instant if required without affecting the global equations system. This partial factorization procedure can further be enhanced if it is to be combined with an iterative update procedure. (4) In the case of short time duration transients mixed time explicit-explicit is most effective. It is even more attractive if the mixed time explicit-explicit procedures are unconditionally stable [8].
5. Numerical results
A two-dimensional finite element pilot computer code incorporating the methodologies described in previous sections has been written to evaluate the performance of these mixed time finite element algorithms. Three numerical examples are presented to demonstrate the accuracy, stability, and efficiency of these proposed methods. All computations are performed on a CDC Cyber 170/730 computer in single precision (60 bits per floating point word), and lumped capacitance matrices are used throughout.
214
W.K. Liu, Y. F. Zhang, Mixed implicit-explicit
algorithms
EXAMPLE 5.1.Comparison with finite diflerence solution. The finite element mesh as shown in Fig. 2 consists of 100 uniform square elements (each with 1, = I, = 10). The mesh is deliberately divided into four element groups. (1) Group 1 is integrated implicitly with a group time step of 2At. (2) Group 2 is integrated explicitly with a group time step of 2At. (3) Group 3 is integrated implicitly with a group time step of 4At. (4) Group 4 is integrated explicitly with a group time step of At. We called the above partition a 21-2E-41-1E partition. Observe that group 3 (41) is chosen to be at the center of the square. Its purpose is to confirm the generality of our modified column equation solver. As described in [3-61 the variable column heights (which affect the solution time) depend on the nodal numbering. Furthermore, no forward reduction and backsubstitution are performed until every 4At time interval. As can be pictured from Fig. 2 the column heights of the ‘4At’ equations are coupled to the ‘At’ and ‘2At‘ equations, therefore unlike the one-dimensional case (see [4] for a discussion) a special forward reduction and backsubstitution procedure as described in Section 2 is required for efficiency. The thermal diffusivity a is chosen to be 6.25 for all the elements. Node 1 is kept at a constant temperature of 0.0 while all the other nodes are subjected to a constant initial temperature of 0.1. All four sides are insulated. A time step of At equal to 1 is employed for the transient analysis. The temperature-time histories of nodes 3, 46, 57, and 67 are depicted in Fig. 3. In order to confirm the accuracy of these solutions the mixed time finite element solutions are compared to those obtained using a finite difference 21 X 21 (400 zones) grid with a single time step of At equal to 1. The finite difference solutions of these four nodes are also depicted in Fig. 3. As can be seen the comparison is quite good despite the crudeness of the mesh (as compared to the finite difference method) and the larger time steps used in the finite The to 1, 2, and 4 for groups 4, 1 and 2, and 3 respectively). element mesh (1 as compared maximum relatively difference between the two solutions is no more than 2%.
Fig. 2. Problem
EXAMPLE
statement
and finite element
mesh of numerical
Example
5.1.
5.2. Demonstration of accuracy, stability, and ejJiciency. The problem statement is depicted in Fig. 4. The finite element mesh consists of 400 uniform size rectangular elements (each has 1, = 5 and 1, = 10) and 486 nodes (81 in x-direction and 6 in y-direction) where the x-axis is defined by joining node 3 to node 17. The material properties of this finite element
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit
0
100
800
TEE
7.ooo 4
0
Fig. 3. Comparison
215
algorithms
100
lEl3
between finite element and finite difference solutions.
TIME
Fig. 4. Problem statement
and finite element mesh of numerical Example 5.2.
800
216
W.K. Liu,
Y.F. Zhang.
Mixed
~~~fj~it-~~pficit
u~gorit~~.s
model are composed of four different thermal diffusivities of O.l2Sa, 0.25~. a and 40~ respectively as shown in Fig. 4. a is chosen to be 12.5. In order to simulate a realistic three-dimensional bandwidth using this two-dimensi~~nal finite element model we number the nodes along the x-axis. The heat load which is also shown in Fig. 4 is applied at node I. The initial temperature for all the nodes is 0.1. All four sides are insulated. AS mentioned earlier, due to the different thermal di~usivities, a single time integration method is inefficient for this problem. For example, if the explicit method is employed, a time step of 0.025 is required so that 16000 steps (as compared to 400 steps) are required for this problem. Therefore the purely explicit integration method is not being considered here. If the implicit method is employed, we can employ a much larger time step though the wide bandwidth and demanding computer storage (e.g. 3288 vs. 33771 single precision words) may decrease the advantage of the implicit method. This can readily be seen in Table 4. Hence the mixed time implicit-explicit method is best suited for this exampfe. A totaf of five computer runs are made to demonstrate the accuracy, stability, and efficiency of the proposed mixed time methods and they are: (1) lE-4E-81-8E with At = I; (2) 11 with At = 1; (3) 11-lE-21-2E with At = 4; (4) 11 with At = 4; (5) 11 with At = 8. In run 1, groups 1, 2 and 4 are integrated expiicitly with lE, 4E and 8E respectively; and group 3 is integrated implicitly. In run 3, groups 2 and 4 are integrated explicitly with 1E and 2E respectively; while groups 1 and 3 are integrated implicitly with 11 and 21 respectively. The time step At and the element group time steps are chosen according to the critical time-step formula given in Section 3 with Al;equal to 1.0. The temperature-time histories at six different locations for run number 1, i.e. lE-4E-SI-8E, are presented in Fig. 5. These results are virtually identical to those obtained from run number 2, i.e. 11. It should be pointed out that the curves shown in Fig. 5 are plotted for every time interval of 4 except for run number 5 in which the time interval is 8. Since the plot of node 327 is integrated with a time step of 8, therefore the temperature-time history shown is ‘step-wise’ (from time -4X to 240). The solution times, and storage requirements for the above five runs are tabulated in Tables 3 and Table 3 Solution times (in seconds) Method
Item Formulation of stiffness and load Factorization Forward reduction and backsubstitution Total solution time
II (At = 1)
1I- IE-2E-21
78.96 (1.34
21 I.13 11.21
58.13 0.80
52.98 11.55
28.81
10.95 90.25
293.58 515.92
10.90 69.83
74.24 138.77
37.50 77.84
1E-4E-XE-81 (At = 1)
(At = 4)
(AZ
4)
I t.s?l
217
a+
TIME
TIME
TIME
.Tlm
.ooil
TIME Fig. 5. Temperature-time
histories of iE-4E-8I-8E
mixed time method with At = 1.
218
W. K. Liu, Y.F. Zhang,
Mixed implicit-explicit
16
-lI-lE-2I-2E AND 4I(DT=4)
12--
+
r
62
T-
167
_______ ..
lI(DT=8)
4
66
ulgorithms
1
3
2
1
0
3.000
327
1.300
l.ooo
.500
4
KuI
4 0.000t 0
Fig. 6. Comparison of solutions of five integration methods.
!
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit Table 4 Storage (in single precision
219
algorithms
words) Method
Item Number of equations (NEQ) Number of terms in stiffness (NA) Average half bandwidth
lE-4E-8E-81 (At = 1)
(At’: 1)
lI-lE-2E-21 (At = 4)
(At’! 4)
(At’: 8)
486
486
486
486
486
3288
33771
6089
33771
33771
6
69
12
69
69
MB zNA,NEQ
4 respectively. From these tables the advantages of these mixed time integration methods are apparent. The relative accuracy of these five runs are compared and they are shown in Fig. 6. These curves are plotted at every time interval of 8. In summary, from accuracy consideration, runs 1 and 3 are compatible to runs 2 and 4 respectively while run 5 is not as accurate as the others. However, judging from the solution times and storage requirements (e.g. 90.25 vs. 515.92 and 3288 vs. 33771 respectively), the superiority of these mixed time methods to the implicit method is fully illustrated. The potential savings in both computer time and computer storage are apparent in three-dimensional and/or nonlinear calculations. Finally, the critical time-step estimate is also confirmed by this numerical example. EXAMPLE
5.3. Applications to TPS. The space shuttle orbiter is designed to be used for at least 100 missions. Currently, various metallic TPS concepts are being investigated at the NASA Langley Research Center (see e.g. [9, lo] for a discussion). One of the proposed candidates is the titanium multiwall tile. A titanium multiwall tile consists of alternating layers of superplastically formed dimpled sheets and flat septum sheets of titanium foil. As described in [lo] this multiwall concept impedes all three modes of heat transfer--conduction, radiation, and convection. The superplastically formed dimpled sheets and the long thin conduction path tend to minimize heat conduction. The flat septum sheets of titanium foil impede radiation. The small individual volumes created by the dimpled layers virtually eliminate air convection. We employed two integration methods to compute the time histories of the top (node 84) and bottom (node 64) surface temperatures of the 042 coatings and the lower (nodes 11 and 137) and upper (nodes 15 and 57) aluminum skin temperatures for the two-dimensional finite element model described in Section 4. Its dimensions (not to scale and in centimeters) are depicted in Fig. 7. The heat loads which are also shown in Fig. 7 are applied at the top and bottom of the TPS model. The initial temperature for all the nodes is -20. Sides AD and BC are insulated.
220
Fig. 7. Problem
W.K. Liu, Y.F. Zhang,
statement
and finite element
Mixed implicit-explicit
mesh of a two-dimensional
ulgorithrns
thermal
protection
systems
(TPS).
As described in Section 4 407317 steps (as compared to 1000 steps) are required if the purely explicit method is employed. Therefore, we do not consider the purely explicit time integration method. We proposed to use lI-lE-4E-41-41 time integration with a time step of 16.71 for the 042 coating and part of RSI, and a time step of 66.84 for the other part of RSI, RTV, felt, Al and air. The computed results are then compared to those using purely implicit method with a single time step of 16.71. The temperature-time histories (plotted at every time interval of 66.84) at six different locations (042 coatings and aluminum skin) are presented in Fig. 8. These results are virtually identical to those obtained from the purely implicit case, i.e. 11. The solution times and computer storage requirements for these two runs are tabulated in Tables 5 and 6 respectively. Even though we use only a total of 147 nodes and 120 elements (a very small mesh as compared to a typical 3D finite element model of 3000 nodes and 4000 elements) we gain a factor of about 1.7 in solution time and a factor of about 1.5 in computer storage. Again, the advantages as well as the stability of these mixed time methods are fully ilfustrated.
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms loo1X-lE-4E-41-41
85--
30-m
-S--
_I
-40-r 0.000
LOWER
-5-w
AL
NODE 137 4.000
8.000
TIME
LOWER
AL
J
12.000
16.000
20.000
I,
-40 0.000
4.000
(Xl 0”)
8.000
TIME
12.000
16.000
20.000
(Xl 03)
00
40
20
0
UPPER
AL
-20
UPPER
AL
4.000
8.000
If 0
4.000
8.000
TIME
12.000
16.000
21
J,
10
3
(Xl 03)
TIME
12.000
16.000
20 00
(Xl 03)
3ooo-
I
I 170--
8840~-
lOo--
1480~-
LOWER
042
UPPER so-
720--
64
-40-c 0.000
04’2
4.000
8.000
TIME
12.000
16.000
2o.oOa
,
o.ooo 4.000
(Xl 03)
Fig. 8. Comparison
84
-40 8.000
TIME of solutions
of two integration
methods.
12000
(Xl 0’)
l&m0
20
222
WK. Liu, Y.F. Zhang, Mixed ~rnp~ic~t-explicit ~lgo~~thrns
Table 5 Solution times (in seconds) Method II-lE-4E-41-41 (At = 16.71)
Item Formulation
(At =1:6.71)
of
stiffness and load Factorization Forward reduction and backsubstitution Total solution time
Table 6 Storage (in single precision
105.12 0.36
162.54 0.57
45.05 160.42
91.00 263.35
words) Method
Item Number of equations (NEQ) Number of terms in stiffness (NA) Average half bandwidth MB “=”NA/NEQ
II-lE-4E-41-41 (At = 16.71)
(Al =1:6.71)
147
147
2121
2933
14
19
6. Conclusions The computer implementation aspects and the evaluation of the performance of mixed time implicit-explicit finite element as applied to multi-dimensional transient thermal analysis are presented. The paper focuses on applications where the thermal model would normally have different time scales of thermal response. Three two-dimensional examples are presented to ihustrate the approach: (1) comparison with finite difference method, (2) demonstration of accuracy, stability, and efficiency, and (3) application to thermal protection systems. The examples show: (1) the superiority of mixed time methods to a single integration method (either implicit or explicit), (2) potential savings in computer time and computer storage, (3) the accuracy and stability behavior of mixed time finite element, and (4) the importance of these mixed time methods as applied to three-dimensional and/or nonlinear thermal analysis.
W.K. Liu, Y.F. Zhang, Mixed implicit-explicit algorithms
223
Acknowledgment
We wish to thank Professor T. Belytschko for his many helpful discussions during the course of the present study. We also wish to express our gratitude to the computer center at Northwestern University for providing computer services. The support of NASA under Grant No. NAG-l-210 to this research is gratefully acknowledged.
References [I] H. Adelman and R. Haftka, On the performance of explicit and implicit algorithms for transient thermal analysis of structures, NASA TM 81880, 1980. [2] T. Belytschko and R. Mullen, Stability of explicit-implicit mesh partitions in time integration, Internat. J. Numer. Meths. Engrg. 12 (1978) 1575-1586. [3] T. Hughes and W. Liu, Implicit-explicit finite elements in transient analysis, J. Appl. Mech. 45 (1978) 371-378. [4] W. Liu, Development of mixed time partition procedures for thermal analysis of structures, Internat. J. Numer. Meths. Engrg. (1983) to appear. [5] W. Liu and J. Lin, Stability of mixed time integration schemes for transient thermal analysis, Numer. Heat Transfer J. 5 (1982) 445450. [6] W. Liu and T. Belytschko, Mixed time implicit-explicit finite elements for transient analysis, Comput. & Structures 15 (1982) to appear. [7] W. Liu andT. Belyschko, Efficient linear and nonlinear heat conduction with a quadrilateral element, Internat. J. Numer. Meths. Engrg. to appear. [8] W. Liu and Y. Zhang, Unconditionally stable implicit-explicit algorithms for coupled thermal stress waves, Comput. & Structures (1983) to appear. [9] H. Kelly, D. Rummler and R. Jackson, Research in structures and material for future space transportation systems-an overview, AIAA Paper 79-0859, AIAA Conf. Advanced Technology for Future Space Systems, 1979. [lo] R. Jackson and S. Dixon, A design assessment of multiwall, metallic stand-off and RSI reusable thermal protection system including space shuttle applications, NASA TM 81780, 1980.