On the static analysis of structures with strong geometric nonlinearity

capum & Strucrurcs Vol. I I, pp. 137-145 Pergamon Press Ltd., 1980. Printed in Great Britain ON THE STATIC ANALYSIS OF STRUCTURES WITH STRONG GEOMETRIC NONLINEARITY R. L. WEBSTER Staff Scientist, Thiokol Corporation/Wasatch Division, Brigham City, UT 84302, U.S.A. (Received 13 June 1979) Abstrret_--Highly flexible structures, such as cable systems or intlatable shells derive their ability to support their own weight and other loads loom some secondary loading (e.g., pretensioning or pressurization). These structures can have very strong geometric nonlinearities when this secondary loading is small. This makes the determination of the initial static shape of such structurea very difficult. Various methods for calculating tile initial static shape are reviewed and their deficiencies noted A novel iterative solution method is then presented which is based on a viscous relaxation concept. The procedure is shown to be highly stable with strong geometric nonlinear&s. Its performana is demonstrated by comparisons with other methods on a simply -able system. INTRODUCIION Structural systems which respond with large deflections are said to exhibit a geometric nonlinearity. The deformations must be included in the expressions of static equilibrium for such structures Typically, the dominant terms in the deformations are rotations; however, problems involving large stretching are also encountered. Geometrically nonlinear problems re- quire that the kinematic relations between strain and displacements contain second order terms [I]. This can be accomplished in a variety of ways. Two popular approaches are the total Lagrangian (TL) and the updated Lagrangian (UL) methods [2]. The TL method uses the undeformed initial state of the structure as a reference for all physical quantities. Internal mechanical energy is described using the product of the 2nd Piola-Kirchhoff stresses and Green’s strains. The UL method is an incremental procedure which refers all physical quantities to the last configuration calculated. Each solution step is followed by a configuration update transformation which changes the reference state. The internal mechanical energy at the reference state is described using the product of the Cauchy stresses and Almansi strains. A Lagrangian solution form is still used for each step, but the reference state is the one at the beginning of the step. bmzrmediate states could be used for the reference, but this is seldom done. The TL form of the equation can be solved using iterations, increments, or combined methods. Incremental or incremental/iterative methods are used with the UL form. The dominant form of the solution iterations is the modified Newton-Raphson method. (See [3,4] for comprehensive descriptions of various nonlinear static solution methods.) One of the most popular incremental solution forms is the so-called incremental self-correcting method. It is essentially a series of simple linearized steps with a single Newton-Raphson iteration after each one. Purely incremental solutions are seldom used because of their tendency to produce growing errors. All of these methods rely on an estimate of the initial stiffness of the structure. This implies that a stable initial state can be defined. In the most common nonlinear structural problems, the initial reference state is quite easy to define. This is so because most structures have sufficient stifhtess to support their own weight. When the weight loading is small compared to the design loads, the dead loaded state is not very different from the Yweightkss” state. Such structures have an intrinsic stili&ss. An ideal flat plate falls in this category. Even though it may be thin and allow dead loaded deformations considerably more than its initial thickness, it is not too d&cult to get to the dead load state. One simply uses the idealixed flat con@uration with no preloads as the initial state and proceeds with the solution paying appropriate attention to non- linearities. When a very thin plate is idealized as a membrane by ignoring the small bending stiffness, the flat initial state has a nonsingular stiffness matrix only if the membrane is pretensioned. Attempts to calculate the dead loaded shape of such a membrane without preloads are likely to be frustrating. Even when a small bending stiffness is retained, the solution may be difficult to obtain. Examples of structures which are likely to present difficulties in obtaining an initial state suitable for structural calculations are Inflatable Shells-covers over sports facilities, portable habitats, etc. Deployable Structures-large, lightweight struc- tural systems for terrestrial and space appli- cations. Gable Structuressuspended roofs, moorings, large ocean structures. Mechanisms and linkages can also be included in this class in some unusual situations. The common denominator in all of these is the fact that in an unloaded, weightless condition they are essentially mechanisms with negligible intrinsic stiffness. This general type of problem has strong geometric nonlinearities even if the initial singularity is overcome. This paper explores some of the numerical methods that have been tried or proposed for this type ofproblem and offers yet another method. For the sake of completeness, some of the descriptive material on 137 138 R. L. WEBSTER the various methods which is availavle from other sources is repeated here. This facilitates the discussion of the deficiencies of the methods when applied to the initial configuration problem. Newton-Raphson method. Its general form is PROBLEM CHARACTERISTICS where The general discrete form of the static equations can be written in residual form as where ‘{R} = ‘{fi - ‘{g} = 0 (1) {R} represents the force residual {f} represents the external loads (g} represents the internal loads bjq}” are the displacements from the initial reference state as estimated by the rtth iteration. ‘,,jq):+’ are the adjustments made in the displace- ments between two successive estimates of the solution. ‘,, [R]n is the tangential stiffness matrix (Jacobian) or some approximation of it for the nth iteration. The procedure begins with an estimate of i(q) and a The left indices are used to identify the contiguration of the system used in the terms. The upper index designates the deformed state. Later, a lower index will be introduced to identify the reference state for displacements and terms which depend on spatial derivatives. This means an Eulerian quantity has both indices the same, while a total Lagrangian quantity has a zero for the reference index. Right indices will also be introduced on iterative solution forms. These indices will be used to designate iteration number. In geometrically nonlinear problems the internal loads depend on the displacements. When the external loads depend on the displacements (as in pressure loading), the loading is nonconservative. Combinations of the two effectss are common. calculation of S[g]“. Equation (2) is then applied recursively untrl both the residual and the displace- ment adjustments are smaller than a preselected level. A divergent condition is assumed when repeated increases in the residual and displacement increments occur. The usual form of the Newton-Raphson procedure is obtained by recalculation of the stiffness matrix at the beginning of each iteration. Modified schemes use approximations for the stiffness matrix and/or hold the matrix constant for a number of iterations. A common approximation is to apply an under- or over-relaxation factor as a scalar multiplier of the matrix. Quite often, geometrically nonlinear problems include a secondary nonlinearity. With the possibility of large displacements comes the possibility of displacement dependent constraints. Typical of these are the contact problems. The most common ones are situations where gaps close or portions of the body move to meet limits on displacements. Problems involving the determination of the contact region between deformable bodies or against a &ted constraint generally require considerable modeling detail and expense to get accurate solutions. The complementary problem where contact is broken and the constraint is removed is also encountered. These situations are approximated by checking the contact constraint condition and modifying it appropriately at each increment or iteration. Although eqn (1) represents the equation of equilibrium, it does not describe a method for finding the values of the displacements which satisfy equilibrium, nor does it give explicit information on how to calculate the internal and external forces. Calculations of the components of the residual could follow various analytical techniques; e.g., finite element methods, finite differences, method of lines, weighted residual methods such as collocation and Gale&in methods. The examples presented later use tinite element models, but any d&r&&ion approach which leads to eqn (1) and allows for convenient calculation of the Jacobian matrix (tangent stiffness matrix) associated with eqn (1) could be used. Associated with the basic procedure is a variety of convergence acceleration schemes which attempt to use results from successive iterations to improve the solution efliciency. One such procedure takes advan- tage of the fact that in some situations the modified Newton-Raphson method produces displacement estimates which oscillate about the correct solution. This is detected by tracking the largest component of 6141. . “+ 1 When oscillation is detected, a new estimate is made using a weighted average of the last two estimates This procedure has proved successful even when the two estimates indicate divergent behavior PI. Two major problems detract from the use of these Newton methods for the class of problems under discussion. First is the conditionally stable nature of these methods. There is an interval of convergence around the correct solution, within which a particular method will converge. When the starting estimate is outside this interval, convergence cannot be obtained without introducing special procedures. These con- vergence enhancement schemes are usually limited in scope and highly specialized so that their effectiveness is quite restricted (e.g., the one mentioned above). The interval of convergence is highly problem dependent and sensitive to variations in the solution method. In most situations this interval is not explicitly calculated. NEWTON-TYPE METHODS By far the most popular method for solving eqn (1) in structural applications is the modified The second problem is the need for making an initial estimate of the solution. This boils down to requiring the very thing that is being sought. A common means of starting a Newton-Raphson solution is to assume ;{q}O is zero and use the initial geometry to calculate a linear step. Quite often the stiffness matrix in the initial configuration is singular or very ill-conditioned. This means that the first step cannot be solved or the solution will be very far from the starting configuration. It is also quite likely to be far from the On the static analysis of structures with strong geometric nonlinearity 139 correct solution. This ill-conditioning can be com- pensated for to some extent by introducing an artificial stiffness. This can be done easily by a technique known as numerical damping [6]. Even with this artificial stiffness there is no guarantee that a solution starting from a very poor quality initial configuration will lead to a first step result that is anywhere near the correct solution. Without that guarantee, there is a strong likelihood that additional iterations will diverge. INCREMENTAL AND INCREMENTAL/TIXRATIVE METHODS An incremental operator applied to eqn (I), such as the first or&r terms of a Taylor Series expansion about a known or assumed configuration, leads to the equation form Here, i [KT] is the tangent stiffness matrix for the displacements i { q} . When nonconservative loads are encountered, the right hand side of the first of eqn (3) is linearized by evaluating ‘+ *{f} at the configuration represented by;(q).+ Repeated applications of eqn (3) with the loads being represented as a collection of load increments gives highly stable solutions as long as the stiffness matrix does not become singular. Unfortunately, accuracy is not a necessary companion of stability. Errors can grow dramatically with this method, and it contains no mechanism to control them. This characteristic in combination with the need for a starting guess which is assumed accurate (i.e., in equilibrium with the initial loading) makes the technique totally unsuitable for getting good initial configurations. An incremental, self-correcting method can be obtained by applying the incremental operator only to the internal forces in eqn (1). The result is (4) Substituting eqn (4) into eqn (1) gives :PGl ‘+*:~4 = ““‘Lfl -ad I+ A&?l =M + ‘+ A:{4). (5) Incremental application of this set of equations amounts to taking a linear step, calculating the residual and adding that residual to the force increment of the next step. Mathematically, this is equivalent to using the incremental form as a predictor for a step and then performing a single Newton-Raphson iteration. This approach couples the stability of the incremental method with the corrective features of the Newton-Raphson method. Since the iterations are not carried to convergence, accuracy can only be assured by small steps. In practice, the adequacy of the step size is usually 7 To be consistent in this case the stiffness matrix should contain a nonsymmetric term representing the change in loading with a change in displacement. This is usually neglected evaluated by solving the same problem with two different step sizes to see if significant changes result in the solutions. Although this method has been used with some success on structures with moderate geometric nonlinear&s, it still does not work well with strong geometric nonlinear&s. The circular problem of needing a good starting estimate to get a good solution is still present. The procedure used in the moderate situations consists of starting from some auxiliary load for which the solution is easily obtained. This loading state is then incrementally relaxed while the real loading is incrementally applied When strong geometric nonlinearities are present, it may be as difficult to get the auxiliary solution as it is to get the real one. Even if an auxiliary solution is obtained, it may take an unrealistically large number of steps to get accurate results; or (more seriouslyA the path between the auxiliary solution and the real solution may traverse regions where the stiffness matrix is singular or very ill-conditioned [7]. Another option available for the form given in eqn (5) is to apply multiple iterations at each load level. This leads to fully combined incremental/iterative methods. Usually the incremental predictor is replaced after the first step with an extrapolation scheme. Linear and quadratic extrapolation forms are used [8]. The rate of convergence is used to adjust the step size for improved performance [9]. These features used in conjunction with the special convergence enhance- ments have been responsible for the success of the Newton-type solutions. Unfortunately, they do little to remove the second objection of the previous section. Initial guesses are still required, and even a very small step size does not guarantee that the first step will lead to convergent iterations. ENERGY SEARCH MRTHODS Recasting the structural problem as a functional minimization and performing that minimization numerically is an attractive technique. The approach usually relies on a minimum energy principle and employs nonlinear programming methods to search for the problem parameters which correspond to the minimum [lo]. An important aspect of this form is that one is dealing with a scalar function of vector quantities. This means that the calculation and storage of large arrays are not required in the evaluation of the functional: a simple accumulation of dot products is used. Unfortunately, the more effective solution procedures do use matrix operations with large arrays in generating trial solutions, so the storage advantage fades away. The field of nonlinear programming is one of the glamour activities in numerical analysis, and it is quite easy to become enamored with its techniques. They are indispensable in optimum design, operations research, pattern recognition, and parameter identification; however, only little scratches have been made on the surface of structural response calculations using these methods. The bulk of nonlinear structural solutions are obtained by other methods. The pioneers in this field have been Schmit and his associates [l 11. Buchholt [12] and Murray [13] have used energy search methods on cable networks. Very little has been 140 R. L. WEBSTER done with the highly nonlinear problems considered here. The major problem encountered in using energy search methods to calculate structural response is the number of generalized coordinates involved. Typical cable and shell problems can use hundreds of degrees of freedom. The more successful applications of nonlinear programming deal with a small number of degrees of freedom. On large problems the procedures can be time-consuming. Choices must be made between the e&iency of the algorithm (which means large array manipulations) or a large number of steps. These choices are not unlike those required in the use of other methods. Because of the potential cost of application and the need for possibly large software developments coupled with the uncertainty of a successful outcome, the energy search methods have not been pursued in this effort. Still the siren call remains, and the author wistfully longs to have the freedom to romp in the open fields and explore what might be there. DYNAMIC RELAXATION A promising procedure for moving towards a correct solution from a nonequilibrium starting guess involves the use of dynamic equations. The starting guess with zero velocities is taken as the initial conditions. The static loads are then applied and held constant as the system is allowed to move dynamically until motion dies out. This effectively avoids the singularity problem if the damping and mass terms are appropriately defined. The promise of success soon fades as one tries to apply the method. If the physical characteristics are used to model the mass and damping, one can expect very large transients which persist for long times. This means excessive solution costs and difficulties jn controlling spurious oscillations in the solution. An obvious remedy is to assume artificial properties which assure strongly damped responses. Use of a highly stable integration scheme is desirable. The question then shifts to how to estimate those artificial properties. The necessity of obtaining contrived terms for both mass and damping properties detracts further from the method. A NEW ITERATIVE METHOD-VISCOUS RELAXATION Neglect of the inertia term in dynamic relaxation leads to a viscoelastic equation form. The specific form chosen is simply eqn (5) with a viscous term added. The result is R1[clr+*~4 +RIvGI ‘+:(4) = ““If) - ‘Id (7) where ;[C] is a damping matrix and the over-dot means differentiation with respect to time. Again, the initial conditions are taken as the guessed configuration at rest. The external forces are held constant at the &sired static values. Consideration of nonconservative forces requires new loading calcu- lations on each step even though the load intensity does not change. The choice of the damping matrix is somewhat arbitrary. To simplify the calculations, a scalar matrix was chosen with the form d[Cl =R1c [II (8) where [I] is the identity matrix, Thus the damping definition is reduced to a single term. The definition of this term still offers some latitude. Three methods were considered in this effort. All of them have the form AC = ‘/l Rf/j. (9) Here ‘p is an arbitrary damping factor which is usually selected to be in the neighborhood of 1.0. It may range from near zero to somewhat greater than 1.0. Its value is adaptively controlled by a process described below to help the solution move efficiently. The,# parameter selection options investigated are (1) i8 = '~R~T,'[K~]'(R}~~R}~~~R) (10) this is one ofthe parameters used by Felippa [6] for numerical damping of the modified Newton-Raphson procedure; (2) Rf/? is the average of the diagonal terms of; [K,. 1; (3) $I is some representative “stretching” stiffness; AE e.g. L for a typical rod or truss member. Values obtained from the first two procedures were found to behave poorly when large configuration errors were present. Bad configurations tend to produce large nonlinear contributions to the stiffness matrix and residual. This results in large damping terms causing sluggishness where rapid movement is desired. The third option tends to give more uniform values of,$ which leaves the task of controlhng the damping to the selection of ‘p. In the vicinity of a convergent solution, the first two options tend to give smaller damping than the third one. This reduced damping behavior near convergence is desirable and may be encouraged further through the ‘P factor. The beneficial nature of the reduced damping can be seen from the structure of eqn (7). As the damping approaches zero with the loading held constant, the equation degenerates to the Newton-Raphson form. This means that near the correct solution the convergence behavior is like that of the Newton-Raphson method; i.e. quadratic. An implicit integration scheme was selected to be used with eqn (7) to give high stability. This is important since the ability to get a solution at any step is at a premium and the accuracy of any one step is secondary. The form of eqn (7) assures that any inaccuracies can be overcome as long as a solution can be obtained on any step and a scheme can be found to force the steps to move in a direction to reduce the error. The particular form chosen is a single step method which uses no information prior to the current state. This is also important since it reduces the amount of storage required (results of previous steps need not be held). Also, there is no reason to believe that in the early stages successive states will have On the static analysis of structures with strong geometric nonlinearity 141 enough in common to be used for extrapolation. The particular form used is based on one of the general formulas given by Smith [14]. The starting point was the following difference formula: ‘+&{q} = ‘{q} + (1 - cc)Att(q} + aAt”*{q} (11) where a is an independently selected parameter between 0 and 1.0. A simple explicit integrator is obtained with a = 0, and a fully implicit form is obtained with a = 1.0. The familiar Crank-Nicolson form results from a = i. Restricting a > 0 and solving eqn (11) for’+*(fj} gives ““(4) =&(‘+‘{q} -‘(q)g-+{4] (12) Substitution of eqn (12) into eqn (7) yields Approximating’ + &‘(f) with the geometry at time t for nonconservative forces linearizes the equation and allows it to be used as an implicit integrator. Those familiar with the Newmark /I integrator may recognixe this equation as a first order parallel to that second order form. The procedures for applying eqn (13) are essentially the same as those for the incremental self-correcting procedure described earlier. When the load level is held over a series of time steps, iterations are produced. Increments are obtained by varying the loading with time. Thus time plays the role of a load level and iteration parameter. The paper by Smith [ 141 and some preliminary stability evaluations suggest a near 1 .O is a good choice for accuracy and stability. Experience with numerical test cases tends to confirm this choice. Detailed stability and convergence studies involving At, a, and ‘P need to be performed before explicit formulas for estimating their appropriate values are obtained. Since this appeared to be a formidable task, it was decided to select a particular a and use simple tests on the velocities at each step to determine if changes in At and ‘P were needed. These tests were devised after noting that with small values of *JJ successive velocity extimates changed signs when the solution estimate had large errors. This is generally a signal of divergence when velocity magnitudes are increasing. Acting on this signal, the value of’p is increased by a factor of ten and the iterations restarted. However, if these alternating signs occur with decreasing amplitudes, it is taken as a signal to increase At. A factor of two was used in this work. When the velocity pattern shows monotonic decline in magnitudes, the solution is assumed to be well behaved and approaching the solution. If the velocity decay is too slow, then ‘~1 is decreased to speed it up. The amount of the decrease depends on the previous history of decreases and increases to avoid cyclic changes. It will generally be by a factor of from 2.0 to 8.0. These procedures are heuristic and by no means optimal, but they have proven to be effective. Convergence to the static solution is accompanied by simultaneous occurrence of small velocities, position change, and residual. The sixeofeach of these is estimated using the folowing norms [15]: 112 (14) Ir N \ 112 (15) Ad = [iizl (‘+*‘4i -‘4i)‘) A, = f ,i (‘+“A - ‘qi)’ 112 (16) I 1 where N is the total number of degrees of freedom. The primary test for convergence is on the residual. When AR is less than a small tolerance value, convergence is assumed. Often the values of 4 and A,, will be small when A, is not. In this situation the damping is reduced by a factor of ten and the iterations continued. This results in the solution procedure becoming more like the Newton-Raphson method as noted above. The overall procedure is summarized in the flow chart in Fig. 1. DEMONSTRATION PROBLEMS A set of numerical experiments is .presented here to demonstrate the effectiveness of the viscous relaxation method. The test article is a single span cable with a large sag to span ratio. This particular problem has an analytical solution when cable elongation is negligible. This is the well-known catenary equation. The interest here is not in modeling catenaries but in the use of a simple problem with very large geometric non- lineaiities. In most of the test cases, the existence of the analytical solution was ignored and somewhat irrational choices for the starting “guess” were used. This was done to be more representative of the quality of initial guesses for more practical situations where the geometric characteristics are significantly more complicated. Three categories of tests were run Comparative tests-The same problem was run usmg Newton-Raphson, modified Newton-Raphson, dynamic relaxation, and vis- cous relaxation methods. Stability tests-Various values of a were tested to investigate convergence and stability behavior of the viscous relaxation method. Initial guess tests-Various poor quality initial guesses were tried to demonstrate the effectiveness of the viscous relaxation method. The characteristics ofthe problem and the final state are represented in Fig. 2. Only ten straight-line elements were used to avoid computation expense. The detail used here is considered to be suEicient for demonstration purposes. The comparison tests were run from an unstressed initial configuration with the line deployed as a straight horizontal line. One end was fixed in space and the other was constrained against vertical motion. 142 R. L. WEBSTER Fig. 1. Flow chart of viscous relaxation solution. Unstretched length* 200 ft W 3 ht per unit Ien th.O.10 Ib/ft Stif ness (ELI)-1.0~ 80 ’ lb Horizontul load, HsS.7735 lb span, Max x.152.2 ft sog, Max Y=SB.O ft Fig. 2. Characteristics of test problem. On the static analysis of structures with strong geometric nonlinearity 143 Table 1. Solution comparisons Solution form Special features Results Newton-Raphson Modified Newton-Raphson Incremental self- Correcting followed by full Newton-Raphson Dynamic relaxation Viscous relaxation Single step update [&I every iteration /1= 1.0 /J = 0.01 jl= 0.001 p = o.ooo1 jJ = 0.00001 single step update [&I after every 5th iteration /l = 1.0 /l = 0.01 /l = 0.001 /I = 0.6001 10 Increments N-Rp = 1.0 Natural properties (Sea water) Viscosity multiplied by 100 Viscosity multiplied by 10,000 O/l = o.m1 Initial At = l.Osec a= 1.0 a = 0.75 a = 0.50 a = 0.25 a = 0.10 Failed with increasing norms, small movement, many iterations. Failed with increasing norms, small movement, many iterations. Failed with increasin g norm% small movement, many iterations. Failed with in creasing norms, reasonable movement. Failed with negative stiksses on 3rd iteration due to large movement and irrational tensions on 2nd iteration. Not converged after 50 iterations, small movement. Not converged after 50 iterationa small movement Not converged after 50 iterations, small movement Failed due to large erratic displacements, after first stiffness update iteration. Incremental over shoots or under shoots. depending on initial conditions. Inconsistent tension patterns. N-R fails with increasing norms. Not converged after 660 time steps, Oscillating velocities, 11.58 CPU sec. Not converged after 660 time steps, Oscillating v&cities, 11.47 CPU sec. Not converged after 660 time steps, Oscillating velocities, 12.23 CPU set, small movement. Converged in 28 iterations 0.55 CPU sec. 34 Iterations 0.62 CPU sec. 40 Iterations 0.75 CPU sec. 34 Iterations 0.67 CPU sec. 34 Iterations 0.70 CPU sec. Note: CPU Set Represents the Central Processor Unit Execution Time on a CDC 7600 Computer. Table 1 summarizes the comparison tests and the results obtained. The superiority of the viscous relaxation is clear. It was apparent for this particular problem that strategies could be devised to get convergent Newton-Raphson solutions with (and possibly without) increments. This would require some ingenuity in setting up the initial conditiona damping controls, and load variations. Experience has shown this approach to be highly problem dependent. The difficulty with the Newton-Raphson solution was the sensitivity of the particular algorithm to non- monotonic norm behavior. It is assumed that repeated norm increases above a specitic size means divergence. Effective strategies for identifying divergent norm changes and reacting appropriately would be needed. This is precisely what is done in the viscous relaxation form. Some typical results when the Newton-Raphson solutions were aborted are shown in Fig. 3. In general, the full Newton-Raphson procedure performed better than the modified form, but with approximately 50% higher costs. Each of the dynamic relaxation runs had similar behavior in that all three trials weremoving toward the solution with oscillatory behavior. The increased viscosity slowed the progress of the solution, All three were stopped when the solution became significantly more costly than the viscous relaxation solution. Fig. 4 compares the results after 660 time steps. The same problem was used to test the sensitivity to the choice of the integration parameter, a. No dramatic sensitivity was found (Table 1). There was a tendency to adjust the damping more with the smaller values of a. It appears that more efficient damping adjustment schemes could be devised which may reduce solution iterations with low a and perhaps improve the performance with high a. The best choice for a in the present formulation appears to be 1.0. IA.44 R. L. WEBSTER -30 -40 -50 -60 -70 -80 2Dr 30 40 50 60 70 80 90 102 110 120 130 140 151~ 160__170 180 11 11 11 11 190 ZOO(FTI / 1 / j ---_. if- ---.- ---.__--.- 41STlTERATlON WITH 0.001 DAMPING FACTOR (FAILED ON INCREASING \\/ \ i CONVERGED SOLUTION / “L~OTH ITERATION WITH O.ooOl DAMPING FACTORIFAILED ON INCREASING NORMS1 , On the static analysis of structures with strong geometric nonlinearity ,z ITERPTION NUMBER SLACK ELEMENT Fig. 5. Viscous rdaxation progress from very poor start. 145 and manipulated based on premises not rigorously justified. Such justification should be sought with the aim of defining a more optimal procedure. It is hoped this paper will stimulate that effort. This viscous relaxation procedure may have applications beyond the initial configuration problem it was developed for. Some large displacement problems start from a reasonably stable initial state but soon develop unstable behavior and a subsequent re-stiffening (e.g. post-buckling response of shell structures). When the response after very large displacements is of interest, it is very costly to follow the entire response accurately. It may be effective to use the incremental or incremental self-correcting method to get out to the region of interest and then use viscous relaxation to correct for equilibrium errors. The usual Newton-Raphson forms often diverge when this approach is taken. The concepts embodied in this viscous relaxation procedure are certainly not new. It is not unlikely that others have developed similar formulations; however, to the author’s knowledge the particular form and the way it is applied are unique. It is interesting to note that with a = 1.0 the integration equation (13) degenerates to the initial value marching operator used extensively in heat conduction analysis. Perhaps not too coincidentally, the fully implicit (a = 1.0) form is equivalent to the damped form of the Newton-Raphson method used by Felippa [6]. The main distinction is in the step size dependence of the damping and the use of the velocity vector to adapt the solution to the characteristics of the problem. The similarities between this viscous relaxation form and the damped Newton-Raphson form were not apparent until the analytical developments noted here were completed. In the final analysis the important results are not the similarities but the differences between the two methods. Early efforts to find a stable solution for the initial configuration problem focused on the damped Newton-Raphson approach. No reliable scheme for selecting the damping to give stability in extreme cases was apparent. With the introduction of the velocity terms (a natural product of the viscous form), the information needed for the control of the damping was found. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. REFEUENCES Y. C. Fun& Foundations qf Solid Mechanics. Prentice- Hall, New Jersey, (1965). K. J. Bathe, E. Ramm, and E. L. Wilson, Finite element formulation for large deformation dynamic analy.sia Inc. J. Numerical Methods in Engineering, 9,353-386 (1975). J. R. Tillerson, J. A. Stricklin, and W. C. Ha&r, Numerical methods for the solution of nonlinear problems in structural analysis. Numerical Solution @ Nonlinear Structural Probkms. ASME. AMD Vol 6. 67-101 (1973). P. G. Bergen and T. tireide, A comparative study of different numerical solution techniaues as amlied to the nonlinear structural problem. C&p. Me%. in Appl. Me& Engng 2, 185-201 (1973). R. L. Webster, An application of the 6nite clement method to the determination of nonlinear static and dynamic responses of underwater cable structures. Ph.D. Thesis, Cornell University (Jan. 1976). C. A. Felippa, Finite element analysis of three dimensional cable structures. Proc. of the 1st International Co@ewnce on Computational Methods in Nonlinear Mechanics, Austin, TX, pp. 311-324 (!kp 1974). R. L. Webster, Finite element analysis of deep sea moors and cable systems. ASCE Fall Convention, Preprint 3033 (Ott 1977). B. 0. Almroth and F. A. Brogan, The STAGS Computer Code. Lockheed Missiles and Space Company, Inc., Sunnyvale, CA, Report No. LMSC-D558853 (1977). B. 0. Ahnroth, P. Stem, and F. A. Brogan, Future trends in nonlinear structural analysis. Comput. Stitures 10, 369-374 (1979). J. T. Oden, Finite Elements of Nonlinear Continua. McGraw-Hill, New York (1972). R. H. Mallett and L. A. Schmit, Nonlinear structural analysis by energy search. J. Structural Div. ASCE 93(ST3), Paper 5285, pp. 221-234 (Jun 1967). H. A. Buchholt, The configuration of prestressed cable nets. ACTA Polytechnica Scandanavia, Civil Eng and Building Construction Series No. 54, Trondheim, Norway (PB-182951) (1968). T. M. Murray, Application of direct energy mi&nization to the static analysis of cable supported structures. Ph.D. Thesis, Univ of Kansas Lawrence, KA (1970). J. M. Smith, Recent developments in numerical integration. ASME Paper No. 73-WA/ALIT-6, presented at the Winter Annual Meeting, Detroit, MI (Nov 1973). P. G. Bergen and R. W. Clough, Convergence criteria for iterative processes. AIAA J. lO(n 8). 1107-1108, (Aug 1972).