Efficient Transient Analysis for Large Locally Nonlinear Structures

gardless of the size of the (linear) portion of the model, the motivation exists for the development of a solution method which isolates the nonlinearities, thereby preserving the linearity of the bulk of the model for an efficient linear solution prior to (repeated) nonlinear design analyses. It is also desirable that such a solution procedure not require all model degrees-of-freedom (DOF) to be retained in the analysis; the solution procedure should require only those DOF directly associated with the nonlinear portions of the model.


Introduction
In the design analysis of structural systems subjected to transient loading, the calculation of dynamic response is a time consuming computational task, greatly exacerbated by the presence of structural nonlinearities in the model.With this computational burden, the repeated analyses required in the optimal design of large nonlinear systems are prohibitive.The optimal design of nonlinear shock isolation is an example of such a design analysis.
Given that, in a traditional direct integration solution, a finite element (FE) model is rendered nonlinear by the inclusion of a single nonlinear element, re-gardless of the size of the (linear) portion of the model, the motivation exists for the development of a solution method which isolates the nonlinearities, thereby preserving the linearity of the bulk of the model for an efficient linear solution prior to (repeated) nonlinear design analyses.It is also desirable that such a solution procedure not require all model degrees-of-freedom (DOF) to be retained in the analysis; the solution procedure should require only those DOF directly associated with the nonlinear portions of the model.
Given these motivations, this work develops a new and general formulation for transient structural dynamics with localized sources of structural nonlinearities.The formulation for nonlinear transient structural synthesis provide fast nonlinear (or linear) transient analysis (or re-analysis) with arbitrary loading for large linear structural FE models which can have localized, but generally nonlinear components of arbitrary magnitude, number, and spatial distribution.The nonlinear formulation is an extension of the general linear formulation [1] for transient structural synthesis.
The governing equation to be developed is a nonlinear Volterra integral equation.Integral equations have appeared in mechanics, and specifically, in the structural dynamics literature in a variety of applications.Integral equations in mechanics are most closely associated with Green's function methods [11,13], and boundary integral/element methods [8], and have been applied both to continuous and discrete problems.The use of integral equations for basic problems in continuous-system structural dynamics is outlined in Hurty and Rubinstein [6].This approach recently was applied to structural system identification [10].The use of Green's function and boundary integral methods has appeared in soil-structure interaction problems [7], and in discrete form (impulse response), in [12].

Overview of formulation
The theory defines a transient analysis that is independent of model size, in that the theory is cast in phys-Shock and Vibration 6 (1999) 1-9 ISSN 1070-9622 / $8.00 © 1999, IOS Press.All rights reserved ical coordinates (i.e., non-modal), and the transient analysis is done using only those structural DOF of interest.These DOF must include, as a minimum, those DOF associated with the nonlinear elements, which are treated independently of the (linear) model.Additionally, other DOF for which synthesized response information is desired can be included as needed.Therefore, it is possible to synthesize the transient response for an arbitrarily large model using an arbitrarily small number of DOF, the minimum number defined only by the number of nonlinear elements in the model.This unique feature of the theory has been demonstrated for the linear formulation [2], and in the frequency domain in [3][4][5].Functioning as a re-analysis procedure, the formulation directly calculates the new transient response for a system resulting from structural changes and/or coupling with other structures, without a re-assembly or full re-analysis.These features will be more fully described in what follows.
The synthesis can be formulated as a substructure coupling approach.This feature allows nonlinear elements to be isolated by division of the system into substructures, which provides additional computational advantages in that the synthesis can exploit inherent physical boundaries in a problem, and linear substructures can remain linear.An example of this is a current submarine design concept: modular isolated internal deck structures connected by nonlinear isolators to the hull.The substructure models (each internal deck module and the hull itself) remain linear, and are solved once, prior to optimization.The synthesis would be used to connect the substructures through the nonlinear isolators, and to calculate the combined system nonlinear transient response.Its important to note for this example that prior to optimization, the baseline transient response is calculated only for those linear substructures subjected to excitation, e.g., the hull.Only impulse response functions need be generated for the remaining substructures.Note that the combined system is rendered nonlinear due to the interconnecting isolation, and this nonlinear transient response is described exactly by the governing equations of synthesis, to be developed in what follows.
As mentioned, impulse response functions (IRF) are used to represent all structures.The IRF are most efficiently calculated using modal superposition.The form of the modal IRF assures rapid convergence due to the presence of the (damped or undamped) natural frequency in the denominator.Therefore, convergence of the structure representation is established once, prior to synthesis.However, the degree of convergence required must be determined based on the use of the IRF in synthesis.

Theory
In the context of the physical coordinate synthesis formulation to be developed, a structural system is defined to consist of one or more uncoupled substructures.A single governing equation for nonlinear transient synthesis will be derived and this equation will address each of the following three general analysis categories: (1) Structural modification -the addition and/or removal of linear and/or nonlinear structural elements; (2) Prescribed base motion -application of base motion to structure through linear and/or nonlinear elements; (3) Substructure coupling -the joining of substructures (a linear analysis).
Each of the above analysis categories defines a set of DOF.Referring to Fig. 1, a structural system is shown, comprised of two substructures, each of arbitrary number of DOF.
The subset of DOF associated with structural modifications is denoted the "m-set" and may include DOF from all substructures.The subset of DOF associated with prescribed base motion excitation is denoted the "b-set" and again, may include DOF from all substructures.The subset of DOF associated with substructure coupling is denoted the "c-set" and may include DOF from all substructures.Given that all nonlinear elements are installed in the synthesis, the substructures are linear, and hence the coupling is a linear synthesis [2].The subset of DOF denoted the "i-set" refers to those system DOF about which synthesized response information is required, but are not directly involved with the synthesis, either modification, base excitation, or coupling.The DOF sets are: All three synthesis categories are governed by an equation (to be derived below) of the following general form: This is a nonlinear Volterra integro-differential equation, where F represents a differential operator.The first issue is the question of existence and uniqueness of a solution, and the proof of existence centers around the Banach fixed point theorem.To summarize, the integral operator in Eq. ( 1) provides a contractive mapping on a complete metric space, which implies the existence of a unique solution.If the kernel function {F (t, τ , {x(t)} * )} is C 0 -continuous on a closed interval [0, t] (as in a numerical solution over finite time), and is Lipschitz, the above integral operator is a contractive mapping, with respect to a suitable metric, hence assuring Cauchy convergence and the existence of a solution.While a contractive mapping assures Cauchy convergence, a contractive mapping on a complete metric space assures that convergence (to a fixed point) is implied by Cauchy convergence (with respect to a suitable metric) [12].The contraction of the integral operators is exploited in the development of efficient iterative numerical solution methods.
We now summarize the several classes of governing integral equations, and provide an example of the synthesis of the transient response for a "large" system with nonlinear springs to ground, subjected to a ground shock.This example will demonstrate the convergence of the iterative solution scheme, a result of the contractive nature of the integral operators.

Derivation of governing equation of nonlinear transient synthesis
The total solution for (linear) transient response can be written in terms of the convolution integral: Here, the vectors {x(t)} of physical responses and the vector {f (t)} of excitations are partitioned according to the previously defined sets of DOF, i.e., This partitioning will be implicit in all matrix equations which follow, unless otherwise indicated.Each element of these vectors, Eq. ( 3), in general are vectors as well.
The matrix [H(t)] is the impulse response function (IRF) matrix, any element of which can be written as: where φ p i is the i-th element of the p-th eigenvector of the substructure prior to synthesis, ω r and ω dr are the p-th undamped and damped natural frequencies, respectively, ζ p is the p-th modal damping ratio, r is the number of rigid body modes, and n is the total number of modes.Note that the IRF matrix [H] contains elements from all substructures involved in the synthesis.We can decompose each excitation component into an externally applied portion and a component due to synthesis, as follows.The i-set DOF are by definition subject only to externally applied excitations, so where the superscript "e" indicates an externally applied excitation.The c-set DOF are subject to externally applied excitations as well as coupling reactions hence, where the overstrike ∼ indicates the Coupling reaction, [R] is a Boolean matrix reflecting the equilibrium which exists between the coupled DOF [1,[3][4][5], and the superscript * indicates a synthesized (unknown) quantity.The m-set DOF are subject to externally applied excitations as well as reactions due to the modifications hence where the reactions due to the modifications are an arbitrary nonlinear function f * m (t) of the synthesized displacements, velocities, accelerations, and time.The bset DOF are subject to externally applied excitations as well as excitations due to prescribed base motion acting through an arbitrarily nonlinear structural element, typically involving displacement-and velocitydependent forces only, i.e., Therefore, Eq. ( 2) for the synthesized response becomes, (6) where {x(t)} contains both the initial displacement and response due to externally applied excitations, {x(t)} = {x(t)} t=0 + t=τ t=0 [H(t − τ )]{f e (τ )} dτ , (7) and {f (τ )} * are the synthesized reactions acting on all DOF sets, Equation ( 6) is a nonlinear Volterra integro-differential equation of the second kind.Direct solution is possible for linear problems; for nonlinear problems iterative solutions are required, and these exploit the contractive nature of the integral operators in achieving good convergence properties.

Example
The following example will demonstrate the nonlinear synthesis procedure in the solution of a pre- scribed transient base displacement excitation applied to an idealized "isolated deck" through four nonlinear springs/linear dampers located at the four corner nodes, as shown in Fig. 2. A piece of equipment is mounted elastically on the deck, modeled using a single lumped mass and linear spring.The deck is modeled using four-noded quadrilateral elements.The deck is a square steel plate, 10 inch on a side, and the plate thickness is 0.125 inch.The lumped mass is taken as 5.5% of the total plate mass.The linear spring stiffness is 1000 lbf/in.
The deck is of approximately 51,500 DOF.The prescribed transient base motion y(t) is shown in Fig. 3.The results will be compared to a nonlinear direct transient analysis performed using a widely-used commercial FE code.Solution times for the synthesis and the commercial FE code ("Direct FE") are provided.The purpose of this comparison is to demonstrate the computational efficiency and accuracy of the synthesis as compared with traditional direct integration schemes, as found in common commercial FE codes.The synthesis makes use of the mode frequencies and shapes calculated using the commercial FE code.The direct FE solution is taken as the "reference" in that the synthesis results are demonstrated to be more efficiently obtained, and without the numerical damping found in the direct FE solution.Both solutions were obtained from a SGI Octane TM 195 MHz R10000 computer.The nonlinear spring force-deflection curve is shown in Fig. 4.This curve is used by both the synthesis and the Direct FE (for comparison) in the form of an interpolated table-lookup.
A normal modes solution was performed once in the Direct FE code in order to generate the undamped modes for the synthesis.These modes were not used in the Direct FE solution.The solution of the governing equation for nonlinear synthesis, Eq. ( 6), is outlined in the following steps: (1) As the nonlinear isolators are "installed" by the synthesis, the free-free modes of the platespring-mass system are calculated, and used to generate impulse response functions for the system, at the b-set coordinates (the plate corner points, where the isolators will be installed via synthesis), and at the single i-set coordinate, the vertical displacement of the lumped mass.For this problem, modes up to 6 KHz were retained.(2) An initial guess at the nonlinear isolator force time histories is made (a vector of unit values) (3) The isolator force time histories are convolved with the IRF functions, producing the (nonlinear) b-set response time histories.
The time step used in the synthesis was 8.0e−5 s.The synthesis took 4 min 15 s.For the Direct FE, a time step of 1.0e−4 s was used, requiring 43 min 41 s.Both analyses were run to a final time of 0.04 s.The solution times are summarized in Table 1.For completeness, the nonlinear springs were replaced by linear springs in the Direct FE required, which required 25 min 12 s.
The Direct FE analysis by definition has no modal damping.In the plots which follow, the Direct FE re-  sults are labeled ζ = 0 to emphasize this fact, as a comparison of these Direct FE results with the synthesis results with varying levels of modal damping in the synthesis is made.Again, the isolation damping is identical for the synthesis and Direct FE.
The nonlinear transient out-of-plane response of the deck corner (at isolator) as calculated using the synthesis procedure (ζ = 0) is shown in Fig. 5, compared with the solution obtained from Direct FE.This comparison reveals that the synthesis provides a similar response as Direct FE, with the significant difference in the apparent level of damping.The higher modal content is evident in the synthesis response, while the Direct FE solution exhibits a significant level damping.Given that the isolation damping is identical for the two solutions, one might attribute this difference to the effect of numerical damping of the solution procedure in the Direct FE.
The comparison of the response for the isolated lumped mass is shown in Fig. 6.Again, the synthesis provides a very close response to that of Direct FE, with the difference being a delay, or period elongation in the Direct FE solution, a characteristic attributable to numerical integration procedures such as the Beta-Newmark [1], or simply slight differences in higher mode frequencies.The nonlinear force at the same deck corner node as calculated by the synthesis is shown in Fig. 7.
The synthesis is repeated, but with a modal damping in the deck of 7% (ζ = 0.07), and compared with the same Direct FE results as before (ζ = 0).Fig. 8 shows the nonlinear transient out-of-plane response of the deck corner (at isolator) for both analyses.It is seen that increasing the modal damping in the synthesis brings the synthesized response closer to that of Direct FE, which again leads to the conclusion that the   Direct FE possesses a degree of numerical damping which underestimates peak transient response.Fig. 9 shows the comparison for the isolated mass, with 0% and 7% modal damping in the deck in the synthesis.The amplitudes are closer, and the period differences are slightly reduced.

Discussion of results
The Direct FE model was run at half the time step used above to verify the degree of convergence.The synthesis was run at a quarter of the time step used above and for a greater number of modes to ensure its convergence.
The transient response as calculated by the synthesis in general compares very well with the Direct FE results.The difference is seen to be a slight period elongation and a higher (numerical) damping level in the Direct FE results relative to the synthesis results.The period elongation might be due to differences in higher mode frequency estimates.The conclusion about the numerical damping is based on the fact that the isolation damper (i.e., C) values are the same in both the Direct FE and the synthesis, and the Direct FE deck model itself is undamped.In order to bring the synthesis results closer to the Direct FE results with respect to the sharpness of the curves (damping), it was necessary to introduce approximately 7% modal damping in the synthesis deck model in order for the synthesized response to exhibit similar damping as the (undamped deck) Direct FE response.Hence the conclusion that the traditional direct integration introduces a noticeable numerical damping.To explain this conclusion further, consider Fig. 10. Figure 10 compares the synthesis result with Direct FE for zero isolator and modal damping, i.e., the completely undamped case.Here again, the sharpness of the nonlinear transient peaks is not present in the Direct FE result.
The synthesis solution does not appear to exhibit significant numerical damping and period elongation and this can be justified from the nature of the integral operator defined in the formulation which has all eigenvalues (for the linear solution) equal to one [1].

Summary
A new formulation for locally nonlinear transient synthesis has been presented.The formulation is cast in physical coordinates, and is exact regardless of the form of the nonlinear elements included.Only those DOF directly associated with nonlinear elements need be included in the synthesis, and hence the synthesis process is independent of model size.The synthesis solution has been shown to be very efficient as compared with traditional direct integration procedures.Further work remains regarding verification of the accuracy of the method relative to tradition integration algorithms.A practical advantage of the synthesis is its ability to easily handle locally nonlinear modal transient analyses, as it does not require any "large mass" techniques, nor does it require any linear springs to base.The example demonstrated does not include the static component of the isolator force, and this is a subject of current work.The fast re-analysis offered by this new method suggests its use in the optimization of locally nonlinear dynamics problems, such as optimal shock isolation.

Fig. 1 .
Fig. 1.General structural system for synthesis comprised of two substructures.
} , where * indicates convolution.(4) If the response time histories are converged, i.e., {x b (t)} k+1 ∼ = {x b (t)} k , go to Step 5. Otherwise, go to Step 3. (5) Using the converged b-set responses and the IRF between the b-set coordinates and the i-set coordinates, the i-set responses are calculated using convolution,

Fig. 10 .
Fig. 10.Comparison of synthesis with Direct FE no damping: Plate corner (at isolator) vertical response vs time.

Table 1
Summary of solution times required for one analysis