- Email: [email protected]
Hierarchial implicit dynamic least-square solution algorithm
Compurrrs d Srrucrures Vol. 22. No. 3. pp. 179489. Printed in Great Britain.
HIERARCHIAL
Departments
1986
0045-79491.96 53.00 + .w S 1986 Pcrgamon Press Ltd.
IMPLICIT DYNAMIC LEAST-SQUARE SOLUTION ALGORITHM*
JOSEPH PADOVANand JOSEPH LACKNEY of Mechanical and Polymer Engineering, University of Akron, Akron, Ohio 44325. U.S.A. (Beceived
10 Ocrober 1984)
Abstract-This paper develops an implicit type transient solution strategy which possesses hierarchial levels of application. In particular, due to the manner of formulation, stiffness updating, assembly inversion, solution constraint, as well as iteration are all performed at a localized level. The level of iterative calculations depends on the type of hierarchial partitioning employed, namely degree of freedom, nodal, elemental, material/nonlinear group, substructural, and so on. Since the iterative solution process and application of constraints are applied at a local level, the resulting so-called hierarchial implicit solution algorithm possesses very stable and efficient numerical properties and is highly storage efficient. To demonstrate the scheme, the results of several benchmark examples are presented. These enable comparisons with the Newton-Raphson solved implicit transient solution method. Overall the comparisons illustrate the superior stability and efticiency of the hierarchial scheme. 1. INTRODUCTION
In recent years, the finite element (FE) solution of nonlinear transient structural problems has been receiving increased attention[ 1, 2). Overall, the solution procedures employed can be categorized into several basic overall methodologies namely: (i) explicit type formulations[2]; (ii) implicit type formulations[2]; (iii) implicit-explicit (mixed) schemes[3] and; (iv) element-to-element schemes involving hybrid integration operators[4]. Each of the foregoing approaches has various advantages and disadvantages. For instance, the explicit approach is algorithmically quite simple and relatively efficient in the context of storage requirements. However, with regard to stability, very small time step increments are required to yield converged solutions. In contrast, the implicit type approach is generally quite stable with regard to step size. Unfortunately, the scheme is relatively less efficient algorithmically. This follows from the fact that the Newton-Raphson scheme used to solve the nonlinear formulation requires stiffness matrix updating assembly and global level inversion. Hence, it is very much more storage intensive. For wide-ranging fluid solid interaction problems, hybrid implicit-explicit schemes are employed[3]. They were derived to deal with the different time characteristics associated with the governing fluid and structural field equations. In this context, such schemes have no special advantages for purely structural problems. Lastly. in recent investigations Hughes et a/.[41 have been developing element-to-element-type transient formulations. The thrust of this work is to * This work supported by NASA Lewis under grant NAG3-54.
reduce the storage requirements for large-scale simulations. To date, applications to heat conduction problems have shown improved architectural-type efficiencies. As yet no large deformation applications have been reported. In more recent studies, Padovan and Lackney[S] have developed a so-called hierarchial solution strategy for static nonlinear finite element formulations. Somewhat reminiscent of the element-toelement scheme, the hierarchial approach enables termination, updating, assembly and matrix inversion as well as solution termination to be pursued at a local level, namely degree-of-freedom, nodal, elemental, material group, nonlinear group or in fact any level of partitioning. In benchmarking comparisons with the Newton-Raphson (NR) scheme, the hierarchial scheme was found to have improved stability and storage utility characteristics. In the context of the foregoing, the current paper will develop a hierarchially solved version of the implicit formulation of dynamic nonlinear structural problems. The overall thrust of the work will be to develop an implicit-type solution strategy which possesses improved numerical stability efficiency as well as reduced storage requirements. This will involve an alternative structuring of the dynamic problem. For the present purposes, a dynamic least-square representation of the dynamic field equations will be used to establish the necessary implicit formulation. The resulting nonlinear formulation will then be used to derive the necessary hierarchial solution algorithm. The main emphasis will be to develop an implicit transient solution methodology which enables localized/hierarchial: (i) stiffness updating; (ii) assembly; (iii) inversion, and (iv) iteration. Additionally, the application of individualized local/hierarchial constraints will be possible so as to control/stabilize the solution de-
479
J. PADOVANand
480
J. LACKNEV
velopment. Due to the generality of the scheme a where [Al. [B], [Cl, and [D] are coefficient mawide variety of nonlinear dynamic FE formulations trices defined by the particular scheme employed can be solved. and In the sections which follow, detailed discussion will be given on: (i) an overview of the NR-type AY = Y(t + II) - Ylr). (2.3) solution of the implicit formulation: (ii) establishment of a dynamic least-square formulation: (iii) development of hierarchial solution strategy: and (iv) Based on eqn (2.3). eqn (2.2) reduces to overview of algorithmic structure. Additionally, a series of benchm~ks will be overviewed. These iI- IMlfAlAY + v, [S*]?3 dz*, /tiAr = F(t + If) I iustrate the operational characteristics of the hierarchial implicit solution methodology. Particular - [fMl{[BlY(f) + [C]Qt, + [D]*(t)I. emphasis will be given to dealing with highly non(2.5) linear response problems. t.
O~RV~W
OF CURREM
IMPLICIT APPROACHES
Without going into details, it is worthwhile overviewing the major algorithmic features of the currently widely used implicit dynamic formulation of the governing field equations of structural mechanics. Typically, such a formulation employs a Newton-Raphson scheme to solve the resulting nonlinear equations. 3ased on the usual displacement type formulation[ I], the typical nonlinear FE equations associated with large deformation problems take the form[l, 21 [M]Y +
Iv,fB*lTS dzr, =
F,
(2.2)
To solve eqn (2.2), generally ii is approximated by a wide variety of numerical integration schemes such as in the Newmark[7] and Wilson[8] methods. Such an approximation reduces eqn (2.2) from a second-order ordinary di~erential equation to a nonlinear equation in the variable Y(i + At). In particular, for numerical integration schemes such as the Newmark or Wilson, ii is approximated by the following linear type operator. namely Yfr f Ar) = [AlilY
+ [I?]Y(t) + [C]%‘(r) + [D]ii(t,,
[[M] + [K]]A Y = F(r f At) - [M]{[B]Y(t)
where
(2.i)
where [Ml is the mass matrix, Y the nodal deflection, ( ) the acceleration, S is the second PiolaKirchhoff stress vector, V, is the volume, F the external load vector, [ ]r defines matrix transposition, and [B* ] is a matrix which defines weighting factors redistributing the internal stress state as concentrated forces on the various element nodes. For implicit-type transient solution formulations, the time history is known up to t. To obtain the dependent fields at I f 31, the next time increment, eqn (2.1) is evaluated at this point in time, namely,
= F(t -i- A?).
Lastly, the integral expression appearing in eqn (2.5) is expanded in Taylor series to yield the following implicit Newton-Raphson-type solution algorithm, that is[2]
(2.3)
such that [D] is the material stiffness, [S] is the prestress matrix, and [G] are coefftcient matrices arising from the FE methodology. Noting the form of eqn (2.6), it follows that up dating, assembly inversion, and iteration must be undertaken on a globai level. Such features obviously lead to numerical and storage inefficiencies. Recently, to enhance the convergence and stability of the NR scheme, Wempner,[9] Riks,[lO], Crisfield,[ll]. and Padovan, Tovichaikchakul, and Arechagal6, 12, 131 introduced constraint concepts. Such schemes have greatly enhanced the range of the NR scheme. This is achieved by controlling the allowable load step size through the use of various constraint functions. Regardless of the enhancement to static situations, the application to implicit transient formulation is awkward. This fohows from the fact that the apptication of constraints would cause a continuous readjustment of time step size involving an increased number of stiffness UP date, assembly, and inversion steps. While the solution would be stabilized, in general the numerical and storage efficiencies would be substantially reduced. The foregoing shortcomings of the overall family of NR-type implicit algorithms points to the need for developing alternative schemes which posses both numerical stabilityieffciency and which are not storage intensive. This then motivates the formulation of the hierarchial implicit methodoIogy outlined in detail in the next series of sections.
Hierarchialleast-squaresolution algorithm 3. DYNAMIC
VIRTUAL
LEAST-SQUARE
FE FORML~ATION
Typically, the D’Aiembert principiefl4] is used to formulate the requisite implicit dynamical soiution aigo~thm for transient mechanics problems. To enable the development of the hierarchial scheme noted earlier, the Gaussian[ IS] least-square version of the governing equations of motion is employed. For the present purposes this is achieved through the use of the virtual work principle. In particular, assuming large deformation problems, the requisite virtual work expression takes the form]151
=
av,
Ltj = 1 (tttj + =
’
Q
IMI = S(Y)‘F =
I,
J
avo
s33.
s,:.
s23.
PJ~LVI
(3.8)
S,3)
(3.9)
dv,
GltiSja(Gik+ Ui.k)Ojds
(3. IO)
[B*l = WI + f&lfGf
(3.11)
L = (@I + 4 rB,llcl)Y
(3.12)
such that
+
tgl.il~i~~,
Se = {e(~))=SYS(Y)~e~~).
(3.14)
Equating (3.7) with eqn (3.14). we see that since SY is arbitrary, the dynamic error function E(t) is
T
(3.3) E(t) = e’e =
Xj.sXi.rUsrr
vu[B*]%dv, + [M]Y - F (J
“* (6LtjStj + PoSttttti - SttfSft) dv, 2 ~LiiSj~(~i~
+
ui,t)nj
ds
>
.
(3.4)
Assuming the use of a displacement type finite eiement formulation we have that[l] u = [N]Y,
(3.5)
where [Nl, Y are, respectively, the shape function matrix and the nodal displacement matrix with @I.
To obtain the dynamic least-square error function, we note that SE takes the following quadratic form, namely
(3.2)
{J
UT =
(3.13)
given by the expression Uj.i
where p, xi, and o,, are respectively the current density, the Euler coordinates, and Cauchy stress. Based on eqn (3.1), the dynamical virtual least square expression is given by
-
ST = (Sii, s22,
(3.1)
SlfjS,~(Si~+ Iti,,+)njdS,
where S( ), S;j, Ltj, Sijv po, it, nj, V,, aV, are, respectively, the variational operator, Kronecker delta, Lagrangian strain tensor, second Piola-Kirchhoff stress tensor, initial density, displacement vector, surface normal, initial volume, and associated surface area. Note that for large deformation probiems[l61
SE =
where
LT = (Lg,, L22, L33.LIZ, L23, L,,i.
J
Sij
481
k.
u3)
X
[B*]%
U
V”
du, + [M]Y - F
)
>
.
(3.f.5)
Note the instantaneous solution which minimizes eqn (3.15) represents the true physical solution. Based on eqn (3.15), either explicit or implicit type formulations can be derived. For the present purposes, the implicit type approach will be considered. In this context, given that Y and ail the other dependent variables are known for time t or less, the error function E is evaluated at t + At so as to estimate Y(t + At). The magnitude of the time increment At is dependent on the operator used to approximate the various accelerations and ve!ocities entering E(t + Ar). Evaluating E at t + &r yields the expression E(t + At) = ([x(r + Ar)lY(r + AI) f [M]%!(t+ Ar) - F(t + A# x ([i’?(r + Ar)]Y(r + Ar)
(3.6)
such that ( IT denotes matrix transposition. Employing eqn (3.5), eqn (3.4) reduces to the form
+ [M]P(t
+
At) - F(r + Ar)), (3.16)
where. SE=
[B*J%dv, + [M]Y - F ‘SYS(Y)7 > X
(J c*e
2213-R
v,, IB*17S dv, -t- [Mfii - F
>
[i?(r + At)] = j-p*(r
(3.7)
+
Ar#.%r + At)di;. (3.17)
482
J. PADOVAN and J. LACKNEY
To estimate Y, a variety of numerical operators can be employed. For present purposes we shall use the Newmark scheme, in particular,
ii(t + At) = --& (Y(t + At) - Y(r))
k(f f At) = Y(r) f Ar(1 - i;)ii(r) + CA&t -t At),
(3.19)
where a and t are parameters which can be so adjusted as to improve the stability/accuracy of the approximation. Based on eqn (3.18), eqn (3.16) reduces to the form
E(r + At) =
+
[Ml
[x(t + At)JY(t i- At) AY -&)
-
(a
- 1) y}
- F(t f At) T (x(t + At)]Y(r f At) + [M]
>(
- F(t -I- Ar)
>
The overall algorithm development involves three major phases, namely, (i) quadratic approximation of E(t + At): (ii) Taylor series approximation of E( f + A 1); and (iii) hierarchial structuring of the algorithm. Before starting the development it should be noted that E > 0 for all real Y. In this context for small neighborhoods of the minimizing solution, E( t + At) can be approximated by the following positive definite quadratic matrix polynomial expression, that is B,(Y) = ai + (Y - Yi)‘{Hi]-‘(S - Yi),
(4.1)
where idenotes the iteration count, aithe minimum of E(r + At) after i iterations, and 8i is the ith approximation of E(t + At). The matrix [H] appearing in eqn (4.1) represents the hyperslope of the dynamic error function space. The function E(t + At) can also be approximated through the use of Taylor-type expansions. To simplify such a development, we note from eqn (3.15) that E(r) is in essentially inner product form. In this context it can be recast as follows .v E(Y) = s f%Y)‘, I= I
(4.2)
(3.20)
where N denotes the total number of degrees of freedom. Expanding each of the various E,(Y) terms in truncated Taylor series yields the following approximation of E(t + At) namely
(3.21)
ei+l(Y) = 2 (E/(Yi) I
where AY = Y(r + Ar) - Y(r).
Through the use of the Newmark operator, eqn (3.16) was converted to a nonlinear equation in the variable Y(t + Ar). The next section wili explore the development of an efficient and stable solution strategy for eqn (3.20). Once Y(r + Ar) is obtained, the process is repeated at t + 2At, and so on until the entire range r,, h t s moatis evaluated. 4. HlERARCHlAL
SOLUTION STRATEGY
Earlier we derived an implicit type transient solution scheme involving a dynamic least-square type error function, namely eqn (3.16). By employing Newmark’s operator[7], eqn (3.16) yielded a nonlinear equation in Y(t t Aht), eqn (3.20). To avoid the difftculties associated with NR-type schemes, we shall develop a modified version of the hierarchial scheme of Padovan and Lackney[Sl to calculate Y(t + Ar). As noted earlier, the main thrust of the following development will be to formulate a new solution algorithm which provides for hierarchial partitioning of: fi) stiffness updating, (ii) stiffness inversion, (iii) iteration, and (iv) solution constraint/control.
+ (Y - YifTVEr(Yji))’-I- H.O.
(4.3)
such that V( ) is the multidimensional gradient op erator and as before &+ i(Y) is the ith approximation of E(t + At) with Yi defining the value of Y after i iterations. Equations (4.1) and (4.3) represent alternative functional estimates of E(t + At). They can be used to deveiop a wide variety of different solution ap proximations. For the present purposes we will consider the development a solution estimate with degree-of-freedom level hierarchy. Namely, the approximation process wiil be undertaken on a degree-of-freedom basis. In terms of (4.1) and (4.3). the (i c lfth degree-of-freedom level estimate is structured to include the effects of the ith element of the dynamic error function as well as a scaled amount of the ith approximation of E( t -I Af ). Specifically, e;*l(Y(t + Ar)t = (&(Y,(t + At)) + (Y(t f At) - Yi(t + At))rVEi(Yi(t + At)))’ + A; O;(Y(r + At)),
(4.4)
Hierarchial least-square solution algorithm
where i is cycled through ail the degrees of freedom until convergence is obtained. The parameter Xiappearing in eqn (4.4) is used to control the level of ith approximation history admitted to the solution. Employing eqn (4.1), eqn (4.4) takes the more convenient form
433
(i) Displacement Norm Check I!Yi+l - Yi I// II Yi 11<
tOl.
(4.11)
(ii) Error Function Check (aii+f
-
Cti)iiXi C tOl.
(4.121
Bi+‘(Y(t + At)) = (Ei(Yi(t + At) + (Yft + At) and (iii) Out-of-balance load check. Note the convergence criteria are monitored continuously. Hence, if convergence is ascertained + hjoi(f + At) + (Y(f + Af) during a cycle through the various degree-of-free- Yift + Ar))‘[H;(f + At)]-’ dom level algorithms noted by (4.7)-(4.10), the x (Y(r + Ar) - Y;(r i- At)). (4.5) overall iteration process is stopped and the next time increment is then considered. - Yi(t + Ar))‘rEi(Yi(f + Jr))’
To obtain an estimate of l’i+I from the foregoing the (i + 1)th version of eqn (4.1) is equated with (4.3, specifically @ic’(Yfr + At)) =I o;c’(t + At) + (Y(t + AZ) - Yi+‘(t + At))‘[Hi+r(t + At)]-‘(Y(t + Al) - Yi+,(t -t At)).
(4.6)
Note two polynomiais are said to be equivalent if they possess equal coefficients. In this context in order to ensure identical levels of accuracy between the solution approximations defined by eqns (4.5) and (4.6), equating coefficients yields the following estimates of oi+‘(t + At), [Hi,‘(t + Al)] and Yi+‘(r + At), namely ai+t(f + At) = Xi ai(t + At) 1 + (&(r + Ar))’ Yi+tft + At) = Yi(t f At) -
’
(4.7)
Yi(t + At) > 1 yi(r + AtI
X [Hi(t + At)]Ei(t + At)VEi(r + At)
(4.8)
X V&(t + Atf([Hg(t + At)JV&(t + At)),
(4.9)
where yi(r + At) is defined by the expression Yi(f
+
At)
=
Ai + (VEi(t + At))” X fHi(t
f
Af)JVEi(r + At)
(4.10)
The foregoing relations define a complete algorithmic set which can be used to estimate the governing field variables at r + At. Note the various equations defined by (4.7)-(4.10) must be cycled for i = I through i = N, the number of degrees of freedom. If convergence is not obtained, the cycle is repeated. For the present purposes, solution convergence is defined by three separate criteria, namely
5. ALGORITHMIC
STRUCTURE
The preceding section developed the degree-offreedom level hierar~hial solution strategy required to handle the dynamic virtual ieast square FE formulation. Before going on to discuss the overall architecture of the scheme, it is perhaps worthwhile to overview the solution flow. Noting the form of eqns (4,7)-(4.10), it follows that the basic flow of control in the hierar~hial implicit scheme proceed as fohows: (i) Commence r + Ar time step. (ii) Cycle through individual degrees of freedom. . Update Y,, [Hi], oi, . . . , etc. l Monitor convergence checks. (iii) Repeat cycle through individual degrees of freedom until appropriate convergence criteria are satisfied. l If satisfied proceed to next time step. l If not satisfied after designated number ot total cycles cancel iteration process. Outwardly, the foregoing solution process apRears similar to the NR approach to implicit integration. Upon closer scrutiny we see from eqns (4.7)-(4.10) that the iteration process is occurring at an element level. Hence, rather than updating an entire stiffness matrix as in the NR scheme, only one row is dealt with at a time, namely VEi(r + Ar). In addition to improving solution efficiency, such an updating process enables a greater improvement in solution accuracy for a given amount of computational effort. Furthermore, by the judicious resealing of hi, the local level constraint, greater control is possible during the solution process. Note the actual choice of hi will be outlined in the benchmarking section of the paper. In addition to the improved control and updating characteristics, the algorithm defined by eqns (4.7)-(4.10) possesses significantly improved storage requirements. Specifically, noting the manner of use and updating of I&], only small partitions need to be maintained in memory. For the degreeof-freedom level implicit algorithm developed herein, Fig. 1 illustrates the changing storage requirements of [Hi] as i is cycled in the interval [I,
J. PADOVANand J. LACKYEY
Fig. 1. Storage requirements of [H] at various iterative stages.
NI. For current purposes, in the manner of Padovan and LackneyfS] [NJ is reset to the identity matrix at the beginning of each fuil cycle through the N degrees of freedom. Such a procedure both maintains the storage efftciency for multicycle iteration situations as well as stabilizes the algorithm. Due to the st~cture of the overalf implicit hierarchial algorithm, the foregoing architectural advantages apply both to linear and nonlinear tran-
AE = 60000X
1 O3
Fig. 2. Material properties, geometry, and ditions of truss roof.
boundarycon-
sient problems. For linear problems only one cycle through all the degrees of freedom is necessary to obtain Y(t + 1t) as we11as the various other dependent variabies. 6. BENCHSIARKING To illustrate the numerical stability and efticiency of the implicit hierarchial solution strategy developed in the foregoing sections. several benchmark examples will be treated. The main emphasis of this work will be to illustrate the inherent stability of the scheme to time step size as well as the range/ magnitude of load excursions. Such behavior will be considered for structures which exhibit significant nonlinearity. In the context of the foregoing remarks, Fig. 2 illustrates the geometry. material properties, and boundary conditions of the truss roof support which serves as the benchmark example. Before considering the numerical efftciency and stability characteristics of the implicit hierarchial algorithm, it is worthwhile overviewing the nonlinear behavior of the truss. This is best achieved by defining the static force-deflection characteristics. Figures 3 and 4 illustrate the trusses static response behavior for sev-
Hierarchial least-square
/ 80--
60--
Fig. 3. Effects of positive rise run ratios on truss roofs vertical deflection.
era1 rise height ratios. As can be seen, for trusses with negative ratios, the response is highly nonlinear of the hardening type. For positive rise ratios, softening behavior is noted. Such properties lead to frequency spectra which are deformation dependent. In this context, as the magnitude of load excursions grow, the dynamic response behavior is either hardened, softened, or combines such effects. Because of this, the amount of nonlinearity excited in a given time step of an implicit type formulation is highly dependent on kinematic, kinetic, and material properties. Of particular importance is the fact that while time increment size may be ad-
Fig. 4. Effects of negative rise run ratios on truss roofs vertical deflection.
solution algorithm
485
equate to represent inertial effects. the level of nonlinearity generated may be sufficient to cause convergence problems for the equations solver. With this in mind, the following benchmarking will indicate the level of nonlinearity excited as well as the relative stability of the Newton-Raphson and hierarchial-type equation solvers. To demonstrate the potential effects of structural nonlinearity on the dynamic response, Figs. 5 and 6 illustrate the transient behavior of the truss to step loading at the apex point. Specifically, as the load amplitude is raised, increasingly higher frequencies are excited for trusses with negative rise ratios. This is clearly seen in the series of responses given in Figs. 5 and 6. Such behavior is a direct outgrowth of the hardening characteristic of the truss. The foregoing results were obtained both by the Newton-Raphson and hierarchially solved implicit dynamic formulation. Excellent correlation was achieved over the entire range of truss behavior of engineering interest. To illustrate the inherent stability of the hierarchial implicit scheme, comparisons were made with the numerical characteristics of the NewtonRaphson approach. Specifically, for a range of given time increment sizes, the magnitude of the step loading was progressively increased. This caused increasingly higher levels of structural nonlinearity to be excited. As can be seen from Tables I and 2, while the hierarchial scheme was able to converge for the full range of load excursions depicted, the Newton-Raphson version of the implicit formulation failed at very low load excursion levels. In fact as depicted in the tables, the hierarchial implicit formulation was able to handle load excursions several orders of magnitude greater with little increase in the total number of iterations. The foregoing benchmarks illustrate the inherent numerical stability of the hierarchial implicit formulation. Such enhanced operational characteristics are an outgrowth of the following features of the scheme, that is: (i) The continuous localized updating of the various field parameters namely El, VEi. Y, . . . , etc. (ii) The imposition of localized constraints through the scaling parameters hi. (iii) The capability of independently varying the constraint level for each localized degree of freedom. Before overviewing the choice of hi, it is worthwhile clarifying the use of the term “history.” For the present purposes two uses are implied, namely: (i) time history as in the succession of events which occur in time; and (ii) iteration history as in the succession of numerically iterated calculations at a given point in time. Based on this definition of usage, it should be noted that for hie[O, l] the iterated history for t + A f is either entirely neglected, partially admitted, or entirely admitted.
486
J. PA~VAN and J. LACKNEY
10 kips
50 kips
Fig. 5. Response of truss roof structure to applied concentrated step loads of 10,50, and 75 kips using hierarchial implicit technique.
Hierarchial least-square solution algorithm
500 klps
Fig. 6. Response to truss roof structure to applied concentrated step loads of 100, 250, and 500 kips using hierarchial implicit technique.
J. PAIXJVAS and Table
1. Effects of load size on solution
J.
stability
LACKNEY for a truss
with a rise-to-run
Time
Step
I
Hierarchial Scheme
Newton Raphson Load (kips)
ratio of 0.
Reason For Lack of Convergence
Time Step
50
200
200
150
200
200
200
1
Out-of-Balance Loads
200
300
1
Out-of-Balance Loads
200
450
1
Out-of-Balance Loads
200
600
1
Out-of-Balance Loads
200
900
1
Out-of-Balance Loads
1200
1
O&t-of-Balance Loads
1500
1
Out-of-Balance Loads
2100
1
Out-of-Balance Loads
Reason For Lack of Convergence
II
I,
II
II
Since criterion functionare in general difficult freedom level. This enables the following: to formulate for nonlinear structural problems, for (i) localized and individualized solution conour purposes a wide variety of numerical experistraint; ments were performed to enable the “best” choice (ii) streamlined storage requirements: of Ai. These included such wide-ranging structures (iii) localized updating and assembly; and as truss roofs with a full range of rise ratios, spher- (iv) handling severe nonlinearities stably and effiical caps, arches, and so on. Based on such work, ciently. the following criteria for the choice of A were asNote the formulation and development of the certained, namely: (i) For degrees of freedom with scheme was undertaken so as to minimize changes no external load application, Aiis set to be 0.1. (ii) needed for retrofitting to standard general purpose For degrees of freedom with applied external loads, (GP) codes. For most available GP codes, the alXiwas set to 0.9. With this choice of Ai, the cross- .gorithmic architecture is of the Newton-Raphsonsection of problems handled by the hierarchial im- type. As such, the tangent stiffness requires global plicit scheme possesses stably converged solutions. updating and assembly. For the hierarchial scheme, Typically, four or less iteration cycles were re- as noted earlier, updating and assembly occurs quired even for large magnitude load excursions. solely at the degree-of-freedom level. In this conAs seen from Tables 1 and 2, such was not the case text, all the necessary updating and assembly profor the Newton-Raphson implicit formulation. cessors are already in place with other minor adjustments required in rearchitecturing storage. Currently, we are working on a scheme which 7. SUMMARY will enable self-adaptive adjustment of A. Work is In the preceding section a hierarchial type im- also proceeding on a hybrid methodology involving plicit formulation was developed and bench- use of NR-type scheme as prediction and hierarmarked. Due to the manner of formulation, the it- chial procedure as a corrector. Such work will be eration process occurs at the local degree of reported in future papers.
Hierarchial least-square
solution algorithm
489
Table 2. Effects of load size on solution stability for a truss with a rise-to-run ratio of 0.0 Hierarchial Scheme
Newton Raphson Reason For Lack of Convergence
Time Step
500
200
1000
3
Out-of-Balance Loads
200
1500
3
Out-of-Balance Loads
200
2000
2
Out-of-Balance Loads
200
2500
2
Out-of-Balance Loads
200
3000
2
Out-of-Balance Loads
3500
2
Out-of-Balance Loads
4000
2
Out-of-Balance Loads
4500
2
Out-of-Balance Loads
5000
2
Out-of-Balance Loads
REFERENCES
Reason For Lack of Convergence
200
Acknowledgment-The first author acknowledges the many fruitful discussions with Chris Chamis of NASA Lewis which stimulated this effort.
1. 0.
Time Step
Load (kips)
II
tural dynamics. J. Eng. Mech. Div.. ASCE 85(EM3). 67-94 (1959). 8. E. L. Wilson, A computer program for the dynamic stress analysis of underground structures. Report No. SESM 68-1, Department of Civil Engineering, University of California. Berkleley (1968). 9. G. A. Wempner, Discrete approximations related to nonlinear theories of solids. Int. J. Solids Struct. I,
C. Zienkieweiz, The Finite Element Method. 1581 (1971). McGraw-Hill. New York (19771. 10. E. Riks, An incremental approach to the solution of 2. K. J. Bathe, Finite Element Procedures in Engineersnapping and buckling problems. In?. J. Solids Struct. ing Analysis. Prentice-Hall, Englewood Cliffs, NJ 15, 529 (1979). (1982). 1I. M. A. Crisfield. A fast incremental/iterative procedure 3. T. J. R. Hughes, Implicit-explicit finite element techthat handles snapthrough. Compttt. Structures 13, 55 niques for a symmetric and nonsymmetric system, in (1981). Recent Advances in Nonlinear Computational Me- 12. J. Padovan and T. Arechaga, Formal convergence chanics (Edited by E. Hinton. D. R. J. Owen, and C. characteristics of elliptically constrained incremental Taylor), pp. 255-267. Pineridge Press, Swansea, U.K. Newton-Raphson algorithms. Int. J. Engng Sci. 20, (1982). 1077 (1982). 4. T. J. R. Hughes, J. Winget, and I. Levit, Element by 13. J. Padovan, S. Tovichaichaikul, and T. Arechaga, Opelement solution procedures. Nonlinear Structural erating characteristics of hyperbolically and elliptiAnalysis, NASA CP 2297 NASA Lewis (1984). cally constrained self adaptive incremental Newton5. J. Padovan and J. Lackney. On the development of Raohson aJaorithms. J. of the Franklin Institute 316. hierarchial solution strategies for nonlinear finite ele197 (1983). ment formulations. Compttt. Strttctures 19, 535-544 14. J. Meirovitch, Analytical Methods in Vibrations. (1984). McMillan, London (1967). 6. J. Padovan and S. Tovichakchaikul, On the solution IS. J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems. Comput. Meth. of elastic-plastic static and dynamic post-buckling colAppl. Mech. Engng 1, 1 (1972). lapse ofgeneral structure. Comput. Structures 16, 199 (1983). 16. Y. C. Fung, Foundations of Solid Mechanics. Pren7. N. M. Newmark. A method of computation for structice-Hall, Englewood Cliffs, NJ (1965).
HIERARCHIAL
Departments
1986
0045-79491.96 53.00 + .w S 1986 Pcrgamon Press Ltd.
IMPLICIT DYNAMIC LEAST-SQUARE SOLUTION ALGORITHM*
JOSEPH PADOVANand JOSEPH LACKNEY of Mechanical and Polymer Engineering, University of Akron, Akron, Ohio 44325. U.S.A. (Beceived
10 Ocrober 1984)
Abstract-This paper develops an implicit type transient solution strategy which possesses hierarchial levels of application. In particular, due to the manner of formulation, stiffness updating, assembly inversion, solution constraint, as well as iteration are all performed at a localized level. The level of iterative calculations depends on the type of hierarchial partitioning employed, namely degree of freedom, nodal, elemental, material/nonlinear group, substructural, and so on. Since the iterative solution process and application of constraints are applied at a local level, the resulting so-called hierarchial implicit solution algorithm possesses very stable and efficient numerical properties and is highly storage efficient. To demonstrate the scheme, the results of several benchmark examples are presented. These enable comparisons with the Newton-Raphson solved implicit transient solution method. Overall the comparisons illustrate the superior stability and efticiency of the hierarchial scheme. 1. INTRODUCTION
In recent years, the finite element (FE) solution of nonlinear transient structural problems has been receiving increased attention[ 1, 2). Overall, the solution procedures employed can be categorized into several basic overall methodologies namely: (i) explicit type formulations[2]; (ii) implicit type formulations[2]; (iii) implicit-explicit (mixed) schemes[3] and; (iv) element-to-element schemes involving hybrid integration operators[4]. Each of the foregoing approaches has various advantages and disadvantages. For instance, the explicit approach is algorithmically quite simple and relatively efficient in the context of storage requirements. However, with regard to stability, very small time step increments are required to yield converged solutions. In contrast, the implicit type approach is generally quite stable with regard to step size. Unfortunately, the scheme is relatively less efficient algorithmically. This follows from the fact that the Newton-Raphson scheme used to solve the nonlinear formulation requires stiffness matrix updating assembly and global level inversion. Hence, it is very much more storage intensive. For wide-ranging fluid solid interaction problems, hybrid implicit-explicit schemes are employed[3]. They were derived to deal with the different time characteristics associated with the governing fluid and structural field equations. In this context, such schemes have no special advantages for purely structural problems. Lastly. in recent investigations Hughes et a/.[41 have been developing element-to-element-type transient formulations. The thrust of this work is to * This work supported by NASA Lewis under grant NAG3-54.
reduce the storage requirements for large-scale simulations. To date, applications to heat conduction problems have shown improved architectural-type efficiencies. As yet no large deformation applications have been reported. In more recent studies, Padovan and Lackney[S] have developed a so-called hierarchial solution strategy for static nonlinear finite element formulations. Somewhat reminiscent of the element-toelement scheme, the hierarchial approach enables termination, updating, assembly and matrix inversion as well as solution termination to be pursued at a local level, namely degree-of-freedom, nodal, elemental, material group, nonlinear group or in fact any level of partitioning. In benchmarking comparisons with the Newton-Raphson (NR) scheme, the hierarchial scheme was found to have improved stability and storage utility characteristics. In the context of the foregoing, the current paper will develop a hierarchially solved version of the implicit formulation of dynamic nonlinear structural problems. The overall thrust of the work will be to develop an implicit-type solution strategy which possesses improved numerical stability efficiency as well as reduced storage requirements. This will involve an alternative structuring of the dynamic problem. For the present purposes, a dynamic least-square representation of the dynamic field equations will be used to establish the necessary implicit formulation. The resulting nonlinear formulation will then be used to derive the necessary hierarchial solution algorithm. The main emphasis will be to develop an implicit transient solution methodology which enables localized/hierarchial: (i) stiffness updating; (ii) assembly; (iii) inversion, and (iv) iteration. Additionally, the application of individualized local/hierarchial constraints will be possible so as to control/stabilize the solution de-
479
J. PADOVANand
480
J. LACKNEV
velopment. Due to the generality of the scheme a where [Al. [B], [Cl, and [D] are coefficient mawide variety of nonlinear dynamic FE formulations trices defined by the particular scheme employed can be solved. and In the sections which follow, detailed discussion will be given on: (i) an overview of the NR-type AY = Y(t + II) - Ylr). (2.3) solution of the implicit formulation: (ii) establishment of a dynamic least-square formulation: (iii) development of hierarchial solution strategy: and (iv) Based on eqn (2.3). eqn (2.2) reduces to overview of algorithmic structure. Additionally, a series of benchm~ks will be overviewed. These iI- IMlfAlAY + v, [S*]?3 dz*, /tiAr = F(t + If) I iustrate the operational characteristics of the hierarchial implicit solution methodology. Particular - [fMl{[BlY(f) + [C]Qt, + [D]*(t)I. emphasis will be given to dealing with highly non(2.5) linear response problems. t.
O~RV~W
OF CURREM
IMPLICIT APPROACHES
Without going into details, it is worthwhile overviewing the major algorithmic features of the currently widely used implicit dynamic formulation of the governing field equations of structural mechanics. Typically, such a formulation employs a Newton-Raphson scheme to solve the resulting nonlinear equations. 3ased on the usual displacement type formulation[ I], the typical nonlinear FE equations associated with large deformation problems take the form[l, 21 [M]Y +
Iv,fB*lTS dzr, =
F,
(2.2)
To solve eqn (2.2), generally ii is approximated by a wide variety of numerical integration schemes such as in the Newmark[7] and Wilson[8] methods. Such an approximation reduces eqn (2.2) from a second-order ordinary di~erential equation to a nonlinear equation in the variable Y(i + At). In particular, for numerical integration schemes such as the Newmark or Wilson, ii is approximated by the following linear type operator. namely Yfr f Ar) = [AlilY
+ [I?]Y(t) + [C]%‘(r) + [D]ii(t,,
[[M] + [K]]A Y = F(r f At) - [M]{[B]Y(t)
where
(2.i)
where [Ml is the mass matrix, Y the nodal deflection, ( ) the acceleration, S is the second PiolaKirchhoff stress vector, V, is the volume, F the external load vector, [ ]r defines matrix transposition, and [B* ] is a matrix which defines weighting factors redistributing the internal stress state as concentrated forces on the various element nodes. For implicit-type transient solution formulations, the time history is known up to t. To obtain the dependent fields at I f 31, the next time increment, eqn (2.1) is evaluated at this point in time, namely,
= F(t -i- A?).
Lastly, the integral expression appearing in eqn (2.5) is expanded in Taylor series to yield the following implicit Newton-Raphson-type solution algorithm, that is[2]
(2.3)
such that [D] is the material stiffness, [S] is the prestress matrix, and [G] are coefftcient matrices arising from the FE methodology. Noting the form of eqn (2.6), it follows that up dating, assembly inversion, and iteration must be undertaken on a globai level. Such features obviously lead to numerical and storage inefficiencies. Recently, to enhance the convergence and stability of the NR scheme, Wempner,[9] Riks,[lO], Crisfield,[ll]. and Padovan, Tovichaikchakul, and Arechagal6, 12, 131 introduced constraint concepts. Such schemes have greatly enhanced the range of the NR scheme. This is achieved by controlling the allowable load step size through the use of various constraint functions. Regardless of the enhancement to static situations, the application to implicit transient formulation is awkward. This fohows from the fact that the apptication of constraints would cause a continuous readjustment of time step size involving an increased number of stiffness UP date, assembly, and inversion steps. While the solution would be stabilized, in general the numerical and storage efficiencies would be substantially reduced. The foregoing shortcomings of the overall family of NR-type implicit algorithms points to the need for developing alternative schemes which posses both numerical stabilityieffciency and which are not storage intensive. This then motivates the formulation of the hierarchial implicit methodoIogy outlined in detail in the next series of sections.
Hierarchialleast-squaresolution algorithm 3. DYNAMIC
VIRTUAL
LEAST-SQUARE
FE FORML~ATION
Typically, the D’Aiembert principiefl4] is used to formulate the requisite implicit dynamical soiution aigo~thm for transient mechanics problems. To enable the development of the hierarchial scheme noted earlier, the Gaussian[ IS] least-square version of the governing equations of motion is employed. For the present purposes this is achieved through the use of the virtual work principle. In particular, assuming large deformation problems, the requisite virtual work expression takes the form]151
=
av,
Ltj = 1 (tttj + =
’
Q
IMI = S(Y)‘F =
I,
J
avo
s33.
s,:.
s23.
PJ~LVI
(3.8)
S,3)
(3.9)
dv,
GltiSja(Gik+ Ui.k)Ojds
(3. IO)
[B*l = WI + f&lfGf
(3.11)
L = (@I + 4 rB,llcl)Y
(3.12)
such that
+
tgl.il~i~~,
Se = {e(~))=SYS(Y)~e~~).
(3.14)
Equating (3.7) with eqn (3.14). we see that since SY is arbitrary, the dynamic error function E(t) is
T
(3.3) E(t) = e’e =
Xj.sXi.rUsrr
vu[B*]%dv, + [M]Y - F (J
“* (6LtjStj + PoSttttti - SttfSft) dv, 2 ~LiiSj~(~i~
+
ui,t)nj
ds
>
.
(3.4)
Assuming the use of a displacement type finite eiement formulation we have that[l] u = [N]Y,
(3.5)
where [Nl, Y are, respectively, the shape function matrix and the nodal displacement matrix with @I.
To obtain the dynamic least-square error function, we note that SE takes the following quadratic form, namely
(3.2)
{J
UT =
(3.13)
given by the expression Uj.i
where p, xi, and o,, are respectively the current density, the Euler coordinates, and Cauchy stress. Based on eqn (3.1), the dynamical virtual least square expression is given by
-
ST = (Sii, s22,
(3.1)
SlfjS,~(Si~+ Iti,,+)njdS,
where S( ), S;j, Ltj, Sijv po, it, nj, V,, aV, are, respectively, the variational operator, Kronecker delta, Lagrangian strain tensor, second Piola-Kirchhoff stress tensor, initial density, displacement vector, surface normal, initial volume, and associated surface area. Note that for large deformation probiems[l61
SE =
where
LT = (Lg,, L22, L33.LIZ, L23, L,,i.
J
Sij
481
k.
u3)
X
[B*]%
U
V”
du, + [M]Y - F
)
>
.
(3.f.5)
Note the instantaneous solution which minimizes eqn (3.15) represents the true physical solution. Based on eqn (3.15), either explicit or implicit type formulations can be derived. For the present purposes, the implicit type approach will be considered. In this context, given that Y and ail the other dependent variables are known for time t or less, the error function E is evaluated at t + At so as to estimate Y(t + At). The magnitude of the time increment At is dependent on the operator used to approximate the various accelerations and ve!ocities entering E(t + Ar). Evaluating E at t + &r yields the expression E(t + At) = ([x(r + Ar)lY(r + AI) f [M]%!(t+ Ar) - F(t + A# x ([i’?(r + Ar)]Y(r + Ar)
(3.6)
such that ( IT denotes matrix transposition. Employing eqn (3.5), eqn (3.4) reduces to the form
+ [M]P(t
+
At) - F(r + Ar)), (3.16)
where. SE=
[B*J%dv, + [M]Y - F ‘SYS(Y)7 > X
(J c*e
2213-R
v,, IB*17S dv, -t- [Mfii - F
>
[i?(r + At)] = j-p*(r
(3.7)
+
Ar#.%r + At)di;. (3.17)
482
J. PADOVAN and J. LACKNEY
To estimate Y, a variety of numerical operators can be employed. For present purposes we shall use the Newmark scheme, in particular,
ii(t + At) = --& (Y(t + At) - Y(r))
k(f f At) = Y(r) f Ar(1 - i;)ii(r) + CA&t -t At),
(3.19)
where a and t are parameters which can be so adjusted as to improve the stability/accuracy of the approximation. Based on eqn (3.18), eqn (3.16) reduces to the form
E(r + At) =
+
[Ml
[x(t + At)JY(t i- At) AY -&)
-
(a
- 1) y}
- F(t f At) T (x(t + At)]Y(r f At) + [M]
>(
- F(t -I- Ar)
>
The overall algorithm development involves three major phases, namely, (i) quadratic approximation of E(t + At): (ii) Taylor series approximation of E( f + A 1); and (iii) hierarchial structuring of the algorithm. Before starting the development it should be noted that E > 0 for all real Y. In this context for small neighborhoods of the minimizing solution, E( t + At) can be approximated by the following positive definite quadratic matrix polynomial expression, that is B,(Y) = ai + (Y - Yi)‘{Hi]-‘(S - Yi),
(4.1)
where idenotes the iteration count, aithe minimum of E(r + At) after i iterations, and 8i is the ith approximation of E(t + At). The matrix [H] appearing in eqn (4.1) represents the hyperslope of the dynamic error function space. The function E(t + At) can also be approximated through the use of Taylor-type expansions. To simplify such a development, we note from eqn (3.15) that E(r) is in essentially inner product form. In this context it can be recast as follows .v E(Y) = s f%Y)‘, I= I
(4.2)
(3.20)
where N denotes the total number of degrees of freedom. Expanding each of the various E,(Y) terms in truncated Taylor series yields the following approximation of E(t + At) namely
(3.21)
ei+l(Y) = 2 (E/(Yi) I
where AY = Y(r + Ar) - Y(r).
Through the use of the Newmark operator, eqn (3.16) was converted to a nonlinear equation in the variable Y(t + Ar). The next section wili explore the development of an efficient and stable solution strategy for eqn (3.20). Once Y(r + Ar) is obtained, the process is repeated at t + 2At, and so on until the entire range r,, h t s moatis evaluated. 4. HlERARCHlAL
SOLUTION STRATEGY
Earlier we derived an implicit type transient solution scheme involving a dynamic least-square type error function, namely eqn (3.16). By employing Newmark’s operator[7], eqn (3.16) yielded a nonlinear equation in Y(t t Aht), eqn (3.20). To avoid the difftculties associated with NR-type schemes, we shall develop a modified version of the hierarchial scheme of Padovan and Lackney[Sl to calculate Y(t + Ar). As noted earlier, the main thrust of the following development will be to formulate a new solution algorithm which provides for hierarchial partitioning of: fi) stiffness updating, (ii) stiffness inversion, (iii) iteration, and (iv) solution constraint/control.
+ (Y - YifTVEr(Yji))’-I- H.O.
(4.3)
such that V( ) is the multidimensional gradient op erator and as before &+ i(Y) is the ith approximation of E(t + At) with Yi defining the value of Y after i iterations. Equations (4.1) and (4.3) represent alternative functional estimates of E(t + At). They can be used to deveiop a wide variety of different solution ap proximations. For the present purposes we will consider the development a solution estimate with degree-of-freedom level hierarchy. Namely, the approximation process wiil be undertaken on a degree-of-freedom basis. In terms of (4.1) and (4.3). the (i c lfth degree-of-freedom level estimate is structured to include the effects of the ith element of the dynamic error function as well as a scaled amount of the ith approximation of E( t -I Af ). Specifically, e;*l(Y(t + Ar)t = (&(Y,(t + At)) + (Y(t f At) - Yi(t + At))rVEi(Yi(t + At)))’ + A; O;(Y(r + At)),
(4.4)
Hierarchial least-square solution algorithm
where i is cycled through ail the degrees of freedom until convergence is obtained. The parameter Xiappearing in eqn (4.4) is used to control the level of ith approximation history admitted to the solution. Employing eqn (4.1), eqn (4.4) takes the more convenient form
433
(i) Displacement Norm Check I!Yi+l - Yi I// II Yi 11<
tOl.
(4.11)
(ii) Error Function Check (aii+f
-
Cti)iiXi C tOl.
(4.121
Bi+‘(Y(t + At)) = (Ei(Yi(t + At) + (Yft + At) and (iii) Out-of-balance load check. Note the convergence criteria are monitored continuously. Hence, if convergence is ascertained + hjoi(f + At) + (Y(f + Af) during a cycle through the various degree-of-free- Yift + Ar))‘[H;(f + At)]-’ dom level algorithms noted by (4.7)-(4.10), the x (Y(r + Ar) - Y;(r i- At)). (4.5) overall iteration process is stopped and the next time increment is then considered. - Yi(t + Ar))‘rEi(Yi(f + Jr))’
To obtain an estimate of l’i+I from the foregoing the (i + 1)th version of eqn (4.1) is equated with (4.3, specifically @ic’(Yfr + At)) =I o;c’(t + At) + (Y(t + AZ) - Yi+‘(t + At))‘[Hi+r(t + At)]-‘(Y(t + Al) - Yi+,(t -t At)).
(4.6)
Note two polynomiais are said to be equivalent if they possess equal coefficients. In this context in order to ensure identical levels of accuracy between the solution approximations defined by eqns (4.5) and (4.6), equating coefficients yields the following estimates of oi+‘(t + At), [Hi,‘(t + Al)] and Yi+‘(r + At), namely ai+t(f + At) = Xi ai(t + At) 1 + (&(r + Ar))’ Yi+tft + At) = Yi(t f At) -
’
(4.7)
Yi(t + At) > 1 yi(r + AtI
X [Hi(t + At)]Ei(t + At)VEi(r + At)
(4.8)
X V&(t + Atf([Hg(t + At)JV&(t + At)),
(4.9)
where yi(r + At) is defined by the expression Yi(f
+
At)
=
Ai + (VEi(t + At))” X fHi(t
f
Af)JVEi(r + At)
(4.10)
The foregoing relations define a complete algorithmic set which can be used to estimate the governing field variables at r + At. Note the various equations defined by (4.7)-(4.10) must be cycled for i = I through i = N, the number of degrees of freedom. If convergence is not obtained, the cycle is repeated. For the present purposes, solution convergence is defined by three separate criteria, namely
5. ALGORITHMIC
STRUCTURE
The preceding section developed the degree-offreedom level hierar~hial solution strategy required to handle the dynamic virtual ieast square FE formulation. Before going on to discuss the overall architecture of the scheme, it is perhaps worthwhile to overview the solution flow. Noting the form of eqns (4,7)-(4.10), it follows that the basic flow of control in the hierar~hial implicit scheme proceed as fohows: (i) Commence r + Ar time step. (ii) Cycle through individual degrees of freedom. . Update Y,, [Hi], oi, . . . , etc. l Monitor convergence checks. (iii) Repeat cycle through individual degrees of freedom until appropriate convergence criteria are satisfied. l If satisfied proceed to next time step. l If not satisfied after designated number ot total cycles cancel iteration process. Outwardly, the foregoing solution process apRears similar to the NR approach to implicit integration. Upon closer scrutiny we see from eqns (4.7)-(4.10) that the iteration process is occurring at an element level. Hence, rather than updating an entire stiffness matrix as in the NR scheme, only one row is dealt with at a time, namely VEi(r + Ar). In addition to improving solution efficiency, such an updating process enables a greater improvement in solution accuracy for a given amount of computational effort. Furthermore, by the judicious resealing of hi, the local level constraint, greater control is possible during the solution process. Note the actual choice of hi will be outlined in the benchmarking section of the paper. In addition to the improved control and updating characteristics, the algorithm defined by eqns (4.7)-(4.10) possesses significantly improved storage requirements. Specifically, noting the manner of use and updating of I&], only small partitions need to be maintained in memory. For the degreeof-freedom level implicit algorithm developed herein, Fig. 1 illustrates the changing storage requirements of [Hi] as i is cycled in the interval [I,
J. PADOVANand J. LACKYEY
Fig. 1. Storage requirements of [H] at various iterative stages.
NI. For current purposes, in the manner of Padovan and LackneyfS] [NJ is reset to the identity matrix at the beginning of each fuil cycle through the N degrees of freedom. Such a procedure both maintains the storage efftciency for multicycle iteration situations as well as stabilizes the algorithm. Due to the st~cture of the overalf implicit hierarchial algorithm, the foregoing architectural advantages apply both to linear and nonlinear tran-
AE = 60000X
1 O3
Fig. 2. Material properties, geometry, and ditions of truss roof.
boundarycon-
sient problems. For linear problems only one cycle through all the degrees of freedom is necessary to obtain Y(t + 1t) as we11as the various other dependent variabies. 6. BENCHSIARKING To illustrate the numerical stability and efticiency of the implicit hierarchial solution strategy developed in the foregoing sections. several benchmark examples will be treated. The main emphasis of this work will be to illustrate the inherent stability of the scheme to time step size as well as the range/ magnitude of load excursions. Such behavior will be considered for structures which exhibit significant nonlinearity. In the context of the foregoing remarks, Fig. 2 illustrates the geometry. material properties, and boundary conditions of the truss roof support which serves as the benchmark example. Before considering the numerical efftciency and stability characteristics of the implicit hierarchial algorithm, it is worthwhile overviewing the nonlinear behavior of the truss. This is best achieved by defining the static force-deflection characteristics. Figures 3 and 4 illustrate the trusses static response behavior for sev-
Hierarchial least-square
/ 80--
60--
Fig. 3. Effects of positive rise run ratios on truss roofs vertical deflection.
era1 rise height ratios. As can be seen, for trusses with negative ratios, the response is highly nonlinear of the hardening type. For positive rise ratios, softening behavior is noted. Such properties lead to frequency spectra which are deformation dependent. In this context, as the magnitude of load excursions grow, the dynamic response behavior is either hardened, softened, or combines such effects. Because of this, the amount of nonlinearity excited in a given time step of an implicit type formulation is highly dependent on kinematic, kinetic, and material properties. Of particular importance is the fact that while time increment size may be ad-
Fig. 4. Effects of negative rise run ratios on truss roofs vertical deflection.
solution algorithm
485
equate to represent inertial effects. the level of nonlinearity generated may be sufficient to cause convergence problems for the equations solver. With this in mind, the following benchmarking will indicate the level of nonlinearity excited as well as the relative stability of the Newton-Raphson and hierarchial-type equation solvers. To demonstrate the potential effects of structural nonlinearity on the dynamic response, Figs. 5 and 6 illustrate the transient behavior of the truss to step loading at the apex point. Specifically, as the load amplitude is raised, increasingly higher frequencies are excited for trusses with negative rise ratios. This is clearly seen in the series of responses given in Figs. 5 and 6. Such behavior is a direct outgrowth of the hardening characteristic of the truss. The foregoing results were obtained both by the Newton-Raphson and hierarchially solved implicit dynamic formulation. Excellent correlation was achieved over the entire range of truss behavior of engineering interest. To illustrate the inherent stability of the hierarchial implicit scheme, comparisons were made with the numerical characteristics of the NewtonRaphson approach. Specifically, for a range of given time increment sizes, the magnitude of the step loading was progressively increased. This caused increasingly higher levels of structural nonlinearity to be excited. As can be seen from Tables I and 2, while the hierarchial scheme was able to converge for the full range of load excursions depicted, the Newton-Raphson version of the implicit formulation failed at very low load excursion levels. In fact as depicted in the tables, the hierarchial implicit formulation was able to handle load excursions several orders of magnitude greater with little increase in the total number of iterations. The foregoing benchmarks illustrate the inherent numerical stability of the hierarchial implicit formulation. Such enhanced operational characteristics are an outgrowth of the following features of the scheme, that is: (i) The continuous localized updating of the various field parameters namely El, VEi. Y, . . . , etc. (ii) The imposition of localized constraints through the scaling parameters hi. (iii) The capability of independently varying the constraint level for each localized degree of freedom. Before overviewing the choice of hi, it is worthwhile clarifying the use of the term “history.” For the present purposes two uses are implied, namely: (i) time history as in the succession of events which occur in time; and (ii) iteration history as in the succession of numerically iterated calculations at a given point in time. Based on this definition of usage, it should be noted that for hie[O, l] the iterated history for t + A f is either entirely neglected, partially admitted, or entirely admitted.
486
J. PA~VAN and J. LACKNEY
10 kips
50 kips
Fig. 5. Response of truss roof structure to applied concentrated step loads of 10,50, and 75 kips using hierarchial implicit technique.
Hierarchial least-square solution algorithm
500 klps
Fig. 6. Response to truss roof structure to applied concentrated step loads of 100, 250, and 500 kips using hierarchial implicit technique.
J. PAIXJVAS and Table
1. Effects of load size on solution
J.
stability
LACKNEY for a truss
with a rise-to-run
Time
Step
I
Hierarchial Scheme
Newton Raphson Load (kips)
ratio of 0.
Reason For Lack of Convergence
Time Step
50
200
200
150
200
200
200
1
Out-of-Balance Loads
200
300
1
Out-of-Balance Loads
200
450
1
Out-of-Balance Loads
200
600
1
Out-of-Balance Loads
200
900
1
Out-of-Balance Loads
1200
1
O&t-of-Balance Loads
1500
1
Out-of-Balance Loads
2100
1
Out-of-Balance Loads
Reason For Lack of Convergence
II
I,
II
II
Since criterion functionare in general difficult freedom level. This enables the following: to formulate for nonlinear structural problems, for (i) localized and individualized solution conour purposes a wide variety of numerical experistraint; ments were performed to enable the “best” choice (ii) streamlined storage requirements: of Ai. These included such wide-ranging structures (iii) localized updating and assembly; and as truss roofs with a full range of rise ratios, spher- (iv) handling severe nonlinearities stably and effiical caps, arches, and so on. Based on such work, ciently. the following criteria for the choice of A were asNote the formulation and development of the certained, namely: (i) For degrees of freedom with scheme was undertaken so as to minimize changes no external load application, Aiis set to be 0.1. (ii) needed for retrofitting to standard general purpose For degrees of freedom with applied external loads, (GP) codes. For most available GP codes, the alXiwas set to 0.9. With this choice of Ai, the cross- .gorithmic architecture is of the Newton-Raphsonsection of problems handled by the hierarchial im- type. As such, the tangent stiffness requires global plicit scheme possesses stably converged solutions. updating and assembly. For the hierarchial scheme, Typically, four or less iteration cycles were re- as noted earlier, updating and assembly occurs quired even for large magnitude load excursions. solely at the degree-of-freedom level. In this conAs seen from Tables 1 and 2, such was not the case text, all the necessary updating and assembly profor the Newton-Raphson implicit formulation. cessors are already in place with other minor adjustments required in rearchitecturing storage. Currently, we are working on a scheme which 7. SUMMARY will enable self-adaptive adjustment of A. Work is In the preceding section a hierarchial type im- also proceeding on a hybrid methodology involving plicit formulation was developed and bench- use of NR-type scheme as prediction and hierarmarked. Due to the manner of formulation, the it- chial procedure as a corrector. Such work will be eration process occurs at the local degree of reported in future papers.
Hierarchial least-square
solution algorithm
489
Table 2. Effects of load size on solution stability for a truss with a rise-to-run ratio of 0.0 Hierarchial Scheme
Newton Raphson Reason For Lack of Convergence
Time Step
500
200
1000
3
Out-of-Balance Loads
200
1500
3
Out-of-Balance Loads
200
2000
2
Out-of-Balance Loads
200
2500
2
Out-of-Balance Loads
200
3000
2
Out-of-Balance Loads
3500
2
Out-of-Balance Loads
4000
2
Out-of-Balance Loads
4500
2
Out-of-Balance Loads
5000
2
Out-of-Balance Loads
REFERENCES
Reason For Lack of Convergence
200
Acknowledgment-The first author acknowledges the many fruitful discussions with Chris Chamis of NASA Lewis which stimulated this effort.
1. 0.
Time Step
Load (kips)
II
tural dynamics. J. Eng. Mech. Div.. ASCE 85(EM3). 67-94 (1959). 8. E. L. Wilson, A computer program for the dynamic stress analysis of underground structures. Report No. SESM 68-1, Department of Civil Engineering, University of California. Berkleley (1968). 9. G. A. Wempner, Discrete approximations related to nonlinear theories of solids. Int. J. Solids Struct. I,
C. Zienkieweiz, The Finite Element Method. 1581 (1971). McGraw-Hill. New York (19771. 10. E. Riks, An incremental approach to the solution of 2. K. J. Bathe, Finite Element Procedures in Engineersnapping and buckling problems. In?. J. Solids Struct. ing Analysis. Prentice-Hall, Englewood Cliffs, NJ 15, 529 (1979). (1982). 1I. M. A. Crisfield. A fast incremental/iterative procedure 3. T. J. R. Hughes, Implicit-explicit finite element techthat handles snapthrough. Compttt. Structures 13, 55 niques for a symmetric and nonsymmetric system, in (1981). Recent Advances in Nonlinear Computational Me- 12. J. Padovan and T. Arechaga, Formal convergence chanics (Edited by E. Hinton. D. R. J. Owen, and C. characteristics of elliptically constrained incremental Taylor), pp. 255-267. Pineridge Press, Swansea, U.K. Newton-Raphson algorithms. Int. J. Engng Sci. 20, (1982). 1077 (1982). 4. T. J. R. Hughes, J. Winget, and I. Levit, Element by 13. J. Padovan, S. Tovichaichaikul, and T. Arechaga, Opelement solution procedures. Nonlinear Structural erating characteristics of hyperbolically and elliptiAnalysis, NASA CP 2297 NASA Lewis (1984). cally constrained self adaptive incremental Newton5. J. Padovan and J. Lackney. On the development of Raohson aJaorithms. J. of the Franklin Institute 316. hierarchial solution strategies for nonlinear finite ele197 (1983). ment formulations. Compttt. Strttctures 19, 535-544 14. J. Meirovitch, Analytical Methods in Vibrations. (1984). McMillan, London (1967). 6. J. Padovan and S. Tovichakchaikul, On the solution IS. J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems. Comput. Meth. of elastic-plastic static and dynamic post-buckling colAppl. Mech. Engng 1, 1 (1972). lapse ofgeneral structure. Comput. Structures 16, 199 (1983). 16. Y. C. Fung, Foundations of Solid Mechanics. Pren7. N. M. Newmark. A method of computation for structice-Hall, Englewood Cliffs, NJ (1965).