- Email: [email protected]
An efficient and accurate iterative method, allowing large incremental steps, to solve elasto-plastic problems
AN EFFICIENT AND ACCURATE ITERATIVE METHOD, ALLOWING LARGE INCREMENTAL STEPS, TO SOLVE ELASTO-PLASTIC PROBLEMS C. NYSSEN Novatome, 20 Avenue Edouard-Herriot, 92350 Le Plessis Robinson, France (Received11 Ma-v 1980) AWraet--In non liear structural anaIysis, an economical computation algorithm should be able to compute for large load increments and the number of iterations per step must remain little sensitive to the increment size. In case of elasto-plasticity, several difficulties are encountered in deriving an efficient solution scheme. A state determination algorithm is proposed which, combined with slight adaptations of the classical Newton-Raphson method, allows to obtain the required property and to compute for accurate solutions. Evenmore, the solutions obtained are reasonably independent on the chosen calculation strategy. Several examples, including combined geometric and efasto-plastic non linear&s, ante of the derived algorithm.
illustrate the pefform-
to include temperature dependence of mechanical properties and kinematic hardening if necessary. If small strains are assumed, the increments of total Green strains may be split into an elastic part and a plastic part :
A large number of theoretical and computation advances in the analysis of nonlinear structures have been made in recent years. Many computer analyses have been carried out and a number of general purpose computer programs for nonlinear elastic and inelastic analyses have been developed. Nevertheless, the choice for an efficient, accurate and economical calculation algorithm to solve elastoplastic problems with a great number of degrees of freedom is still a tricky task [l-3]. It seems to be a general agreement that a good algorithm should be able to deal with large load increments, this requirement being obtained without alteration of the solution accuracy. The main dithculties in case of elastoplasticity arise from the irreversibility of the plastic strains and from the dis~ntinuous change of the material properties. From this, two distinct phases must be considered when deriving a computation algorithm in elastoplasticity : the numerical iterative method to solve the system of nonlinear equations and the numerical algorithm for scourate state dete~ination of the material. As will be shown, these two fundamental aspects are closely bounded. We propose here a particular state determination algorithm which allows to maintain good convergence properties to the Newto~~ph~n method, provided it is slightly modified. It will be shown that the numerical solution computed is quite independent of the particular calculation strategy adopted, which is also an essential requirement.
i,, =g;+ it. The elastic strains are related to the stresses by the classical linear elastic constitutive relations :
where the CiJU_are the elastic moduli. The stresses should always satisfy the plasticity criterium: F(a, f&J=&-X2=0
(3)
with ii2=f(Cij);X=XO+
H’dg,,
where f is a wnvex function which defines the plasticity criterium used. X, is the initial elasticity limit, while dB is an equivalent plastic strain. H’ is the coefficient of strain hardening. During the plastic flow, the stresses must be plastically admissible and satisfy eqn (3). Following the Drucker stability postulate [Sj. the plastic strains are given by
2. SYSTEMOF RQUATIONSTO SOLVE In this analysis, it will be assumed that the elastoplastic materials obey the Prandtl-Reuss-von Mises plastic flow rules [4]. Only isotropic materiab with isotropic hardening and at uniform temperature will be considered in the formulation for brevity. All considerations can be straightforwards generalized
The exact solution of the plasticity problem should satisfy relations (l)-(3) and (5), the equilibrium equations and the appropriate boundary conditions. During plastic flow, the increments should satisfy the consistency condition : a.Ii6.Ii =~xHBP=~~~H~=A~. 63
(6)
64
C. NYSSEN
From (1 ), (2), (5) and (61, /i may always be expressed in terms of strain increments to obtain the classical elasto-plastic incremental stress-strain relations : (7) with
one obtains K’(q)Aq=g-
(8). the I parameter takes a value of unity during plastic flow, i.e. iff(a,)= X2 and f>O: it takes a zero value otherwise. In a finite element static analysis based on the displacement formulation, the equation of equilibrium may be stated as
s
B’(q)a d I’=g, Y where a are the ~rchho~-Trefftz stresses, B is the incremental strain matrix which generally depends on the generalized displacements q, and g is a vector of generalized loads. Since the strains are assumed to be small, the stresses are related to the total strains through the constitutive equations (7), i.e. u=: r D(a)ds,
(10)
Je
where the integration is taken along the deformation path since pIasticity is essentially path dependent. The nonlinear equations to solve in an eiasto-plastic problem may be written
3. REQUIRED PROPERTKES OF AN EFFICIENT SOLUTION METHOD
3.1 Solution proceiiWes It is now generally admitted that an iterative procedure must be used to achieve sufficient accuracy in solving the nonlinear system (11) [ l-3,6-8]. The main iterative methods can be considered as particular cases of the Newton-Raphson method. Suppose that we know an approximation of the solution at iteration k, noted (a, a+41,
(12)
which satisfy the compatibility requirements and the constitutive relations (7). The exact solution may always be written as q=q+Aq
E=E-+AE
~=u+Ae,
(131
where Aq, As and Aa denote the corrections to be done on the solution (12). which are not necessary inlinitesimal. Expressing that the solution (13) must satisfy eqns (ll), it is then possible to compute for the corrections A, after linearization in these corrections. Noting that
B’fq)crdY=r,
115)
where
K’(q)=
In
s 8’
s c
[B’(q)%P3~q)+~~qw1
dV
(161
is the tangential stiffness matrix computed using the tangent material properties evaluated at the approximate state (12). In (15). r states for the out-of-balance forces. Solving (15) for Aq, a new estimation of the displacements is obtained. The corresponding state of the material is then obtained by computing the corresponding strain increments to which the stress increments are related through the elasto-plastic constitutive relations :
From the new stress state, the new out-of-balance forces may be eomputed and the procedure repeated until an appropriate norm of r becomes less than a user fixed tolerance E. It is worth noting that a solution is then reached whatever the stiffness matrix used for the iterations. Constant stiffness method. If the stiffness matrix is formed and decomposed oniy once, the constant stifEness method [9, lo] is recovered. In case of strong nonlinear&s, the method often fails to converge even if an acceleration scheme is used and this method cannot be retained solely for a general purpose computing code. Strict and modified Newton-Raphson method. In the strict Newton-Raphson method, the tangential stiffness is changed at every iteration. A disadvantage of this procedure is that a large amount of compu~~o~l effort may be required to form and decompose the stiffness matrix. But, if the convergence of the method is sufficiently fast, it can be compensated for by saving time in the state determination phase. In fact, when evaluating the computational efficiency of a solution scheme, it should be noted that the computational cost of this phase can be significant as compared with the solution phase cost [l-3]. Depending on the degree of no~in~~ty, computing time can sometimes be saved by reforming the stiffness matrix on each few iterations (modified Newtow Raphson method). It is also interesting to note that the factorized tangent stiffness matrix may also be gradually constructed at low cost while iterating using quasi Newton methods [ 11, 121. These methods seem very promising but were not investigated in this work. When only geometric nonlinear%& are present and in absence of instability, the convergence of the Newto~~p~on method is quadratic and in structural-analysis, convergence is almost always obtained in a few iterations. Due to the ability of the method to converge even in case of strong nonlinearities, the method is widely used in nonlinear structural analysis. When applied to elasto-plastic problems, some difficulties are still encountered due to the irreversibility of the plastic stress and the discontinuous variation of material properties due to eiasto-plastic state transrtion. Its use and efficiency is sometime controverted in
65
An efficrentand accurate iterative method to solve elasto-plastic problems this case [13]. The main problems of the method may be summarized as follows. The first difficulty is inherent to each type of iterative method. To each computing strategy corresponds a particular path in the load space given by B’adV
9= s”
(18)
where the u are the stresses computed during the iterations. This path can be quite sophisticated. Since an elasto-plastic solution is essentially path dependent, the value of the final solution must be discussed. Convergence difficulties may be encountered if the method is applied without special care. In case of combined material and geometric nonlinearities, the convergence of the Newton-Raphson method can excessively slow down and even be lost in particular cases [13,14]. Due to these difficulties, some modifications of the classical Newton-Raphson method should be introduced to retain the decisive advantages of reliability of the method. 3.2 Defmition of an efficient iterative scheme A solution method to be included in a general purpose nonlinear computer program should have a great reliability and offer a great flexibility in the choice of computing strategy since the best one vary with each particular problem. The computed solution should be reasonably independent of the particular strategy adopted. Finally, the method should be economic to use which implies two aspects. Firstly, in order to save engineer time, convergence should be reached whatever the reasonable choice of the load incrementation may be. Secondly, in order to save computing time, the rate of convergence should be sticiently fast to insure that the total computing cost, including the stiffness matrix reformations and the state determinations, should be acceptable. It is now well recognized [3,14,13 that these requirements can only be achieved if the computation algorithm is able to deal with large load increments, provided the number of iterations per step remains quite independent of the increment size. The state determination algorithm should consequently be very accurate, which can significantly increase the computation cost of the corresponding phase and leads to limit the number of iterations per step. It will be shown that the required properties can be reached using a Newton-Raphson type method with an appropriate state determination algorithm.
4. STATE DETERMINATION ALGORITHM
Computing of strain increments from displacement increments involves only kinematics. The problem of computing a stress increment from a given strain increment involves the material constitutive relations at two levels: i.e. how to determine the increment of stress Aa for a given strain increment ALEand how to determine the total increment of stress over a whole increment. 4.1 Stress increment associated to a given strain increment l%e stress increment during plastic flow must be
such that the plastic flow rule and the plasticity criterium should be satisfied with about the same degree of accuracy. Due to the nonlinearity of (3) and (5), the true stress increment must be computed using, t+bZ A#= Ic
D(a) ds.
(19)
Various algorithms have been designed for this purpose. 4.1.1 Single step methods. A 5rst approximation of (19) is obtained using a Euler one step forward integration scheme : AaPz D(a)Ac.
(20)
The fact that the direction of the plastic flow is only correct in the beginning of the increment can lead to important error in the final orientation of the stress vector in the stress space [15, 161. Moreover the final stress a“ will not lie on the yield surface. An additional correction is then needed and generally, the stresses are brought radially back to the yield surface. The plastic 5ow which takes place during this correction is usually completely neglected. So, the radial return seems not very consistent. To eliminate this difficulty, a single step Euler backward method can be used. In this case the flow rule is satisfied only at the end of the intaval considered while the elastic stresses are again projected radially to the yield locus. But in this implicit method, the additional plastic 5ow occuring during this return can easily be taken into account if the von Mises criterium is used [lW7]. Schreyer et al. [15] have demonstrated in a parametric study that for a single step method, the error associated with the implicit scheme is usually much less than the error obtained with the Euler forward method, but if the increments are large, the error can still be unacceptable. Moreover, the convergence properties of the iterative method singularly degrades when using the implicit scheme if the stress state is far from uniaxial. 4.1.2 Subincrementation schemes. Improvements of the accuracy of integration (19) can be achieved by multi-stepping [3, 6, 14, 16, 181. In this case, both implicit and Euler schemes lead to results that may be considered quite satisfactory for most engineering problems, for approximately the same number of subincrements [15], the errors being even less for the explicit Euler scheme in the case considered in [ 151. The algorithm implemented is based on an Euler forward multistep integration scheme since it is the most easy to use in plane stress problem and for multilinear hardening law. It is coupled with a consistent correction to bring the stresses back to the yield surface at each subincrement and is formulated as follows. Since higher derivatives are not available, a constant strain rate is assumed like in most structural programs This assumption may be violated in the true physical loading path. Strictly speaking, small enough load steps should be used. However the proposed algorithm performs very well even for quite large increments of load. When plastic loading occurs. the increment of strain AE is then divided into m equal subincrements. To estimate the parameter m, the tronmture error of a one step integration is first estimated. If a, designate
66
C. NYSSEN
the difference between the stress state reached m a single step and in a double step Euler forward integration, a measure of the troncature error for one step is taken as 2J’ ?crJ/X.
(31 I
It is then assumed that the total error is roughly I,‘m of the single step error if a m substep procedure is used. The number of subincrements is then taken equal to m =2~f’%J,)/XE,
(22)
where E is roughly equal to the requtred relative error level. Sufficient accuracy was always obtained using E=O.OS and the number m was almost always less than 30. More sophisticated algorithms could be used to determine m, based on the variation of the plastic flow direction during the integration process [14, 181, but the simple formula (22) has always proven to be sufficient in the applications. If the strain hardening law is multilinear, a strain subincrement &,,,=W,
(23)
is divided into the required number of substeps so that no discontinuity occurs within as subinterval. For each subinterval, the associated stress increment is computed by Ab,m)=W)Aa,m,,
(24)
where e states for the stresses at the beginning of the subinterval. At the end of each subincrement an additional correction is made to bring the stresses back to the yield surface. To do this, the following scheme, which is consistent with all plasticity laws to the first order, is proposed. In a displacement iterative scheme, this correction denoted 6 must be done at fixed total strains, which implies 6Ee+6&P=o.
(25)
On the other hand, we would like to have f(o+W=(X+W2,
(26)
$7=X’-/(e)+&XH”Si.
(27)
or linearizing
If due account is taken from plastic flow rule (S), it comes that (28) where H is the Hooke matrix. The corresponding correction for the stresses is given by 6a= -62Ha.
4.2 Total stress increment evaluation over an increment From strict application of the Newton-Raphson method the total strain increment at iteration i of increment n should be evaluated by
(391
As demonstrated by (28). plastic flow always occurs during this correction.
130) with at = a:, where a starred value is a value at the beginning of the considered increment. Nevertheless, if (30) is applied without care, the solution obtained may be very path dependent, and depends strongly of the particular calculation strategy adopted [3]. Moreover, fictitious numerical unloadings may occur during the iterative process. If they are considered as irreversible, these can lead to erroneous results especially in cases of combined geometric and material nonlinearities. Evenmore, convergence may be lost due to loadingunloading cycles. TO avoid this difficulty and obtain a reasonable path independent algorithm for a given load increment. many authors [l. 7, 13, 181 have proposed to compute the stress increment by integration on the total strain increment at each iteration. i.e. for iteration i of increment n, on As;= i dAs, .4=,
(31)
where dAslr is the new strain increment at iteration k. Nevertheless, this technique has a serious drawback: the cost of the state determination does not decrease during the iterative procedure since the whole strain increment is always integrated. Moreover. the strain path used to compute for (AC?): differs from the strain path used in the Newton-Raphson iteration for the approximated stress : (Ac“):,=(Aa’):-‘+D[(c”)2
dAs;.
32)
Consequently, the new state of stress may differ from the approximated value by a quantity which is not quadratic in dAe:. This technique can lead to serious convergence difficulties especially in case of combined geometric and material nonlinearities. This fact has led Bushnell [ 141 to separate in this case the geometric and material iterations to maintain the convergence of the iterative scheme. But this procedure may be costly since each material iteration requires several geometric iterations. All these drawbacks may be eliminated in another way. The path dependence of the state determination may be eliminated using an assumption of incremental reversibility. “A point which deforms plastically during one increment is assumed to unload plastically until the plastic work done become again equal to its value at the beginning of the considered increment. Then only will elastic unloading take place”. The corresponding stress-strain relation is given on Fig. 1 for a one dimensional case. Due to this assumption, the stress increment may be computed using (30). The number of integration subincrements decrease rapidly as the iterations converge, which save computing time for the state determination. Moreover, the procedure is more consistent with the Newton-Raphson method and good convergence properties are maintained by this procedure. The solution becomes also reasonably
An efficient and accurate iterative method, allowing large incremental steps, to solve elasto-plastic problems
Fig. 1. incremental reversibility assumption.
independent of the particular adopted for the increment.
computing
strategy
5. IMPROVEMENTS OF I-HE NEWTON-RAPHSON METHOD FOR ELASTGPLASTIC CASES
Two additional adaptations have been introduced in the classical Newton-Raphson method. When the loads are decreased, the plastic points undergo too large unloading when the first iteration of the corresponding increment is made using the tangent stiffness matrix. The convergence of the Newton-Raphson method may then excessively be slowed down. To avoid this difficulty, the possibility to restore the elastic stiffness matrix at the first iteration of each increment of loads was introduced. In the case of continuous loading, this additional iteration is not very time consuming since the elastic sti&ss matrix triangulation, made at the first iteration of the problem, may be saved and does not need to be rerun. Moreover, when the incremental reversibility assumption is introduced, this iteration does not generally change the solution. Secondly, non convergent loading-unloading cycles
Fig. 2. Unloading-loading cycles in a uniaxial case
67
may be observed when a plastic point unloads as far as to reload plastically during the same iteration. This phenomenon was often observed numerically in case of combined material and geometric nonlinearities. To avoid these spurious cycles, the elastic stiffness is automatically restored at this point during the increment until a new elastic state is reached during the iterative process (Fig. 2). In this way, the convergence of the Newton-Raphson method has always been restored in all the cases where the loss of convergence was due to the considered phenomenon. The final algorithm obtained has been tested on numerous examples including thermoplasticity with isotropic or kinematic hardening and combined geometric and material nonlinearities. The convergence rate reached has been always satisfactory even for very large load increments. Moreover. the number of iterations needed per step was nearly independent of the increment amplitude and the solution computed reasonably insensitive to the particular strategy adopted on the linear parts of the loading path. A few examples will now be presented to illustrate the accuracy and the efficiency of the proposed algorithm.
6. EXAMPLES
The proposed algorithm was implemented in general purpose computer programs [ 19, 201. The results presented here were obtained using either of these programs. 6.1 Tension-torsion tube Problem description. The physical problem being solved is that of a thin walled cylindrical tube which is subjected to axial tension and torsion. Since the stress distribution in the tube, remote from ends, is wnstant everywhere, the only requirements for the finite element model is that it represents material subjected to uniform tension and shear. This could be accomplished by a single plane constant stress isoparametric quadrilateral element. The material wnstants are given in Fig. 3. The non proportional loading path is specified in the same figure. The elastolplastic’nonlinearities are very severe for this case, which was retained by Clinard et al. [21] as a benchmarking problem for testing nonlinear computing codes. Results and discussion. To compute the solution, a systematic incremental history was chosen, as if the final solution was unknown. Each linear segment of the loading path was equally divided in n increments. To obtain a reference solution, n was tirstly chosen equal to 80. The solution obtamed agree with the averaged computed solution using ADINA. ANSYS, CREEPLAST, CREEPABSA and PLACRE as given in [21], within 57;. Then n was chosen equal to 5 which corresponds to quite large increments. In this case, two different computing strategies were adopted: pure Newton-Raphson iteration and modified NewtonRaphson method with a first elastic iteration at each step. Axial and shear stress strain results are plotted in Fig. 4. The results agree closely with the reference solution. When the elastic stiffness matrix is used at the tirst iteration, the shear stress-strain curves differ slightly, but the maximum difference is less than 9:$,. Figure 5 gives the convergence curves of the algorithm for a few increments and demonstrates its quality.
c. NYSSEN
Fig. 4. Stress-stun
response III the tension-torsxon
6.2 E&at&dar pixrte under cimcentrratetf lcrad droller ~~crip~io~. A simply supported ~um~ium
circular plate is submitted to a concentrated load at its center. Geometricand material data are given in Fig. 5. Isotropic hardening was adopted using a piecewise linear law with eight segments. The plate was tested by Levine et al. [22], which present aiso computed results in good agreement with the test results. but they use about 100 load steps to reach the finaf load of 1000 Ibs. Computed results using BOSOR 5 are also avaikbke in [ 141 with separate material and geometric iterations. As reported by Bushnell {i4], the stress path differs significantly from radial path near the center of the
test
plate and since b&h geometriai and material nonIinearities are s~gnifi~nt, this is a good figuration an which to verify the analysis and the strategy. Resz&s obtuined. The totai loads of 1000 Ibs was applied in ten, two and one increments respectively. The deflections of the center of the plate is given in Fig. 5, where they are compared to the other available results. Good agreement with the experimental results should be noted. Computed results and CPU times are compared in Tabie f . It is a striking fact that the differences between the central defkctions, computed in the different load incrementation cases, are less than 1”;. The differences on the final stresses computed in one or ten steps may
An efficient and accurate iterative method to solve clssto-plastic problems
Fig. 5. Convergence curves obtained in the tension-torsion
test applying large increments of load.
69
6.3 Eva&don of the burst pressureof a rocker motor Problem dewiption. The last example comes from an industrial study. It consists of a rocket motor case for which experimental results are available. The complete structure shown in Fig. 4 has three different material properties. In the first part, denoted composite material, we have anisotropic material which behaves elastically. The second part, denoted rubber joint, ensures the bounding between the vessel and the polar boss. The third part is made of aluminium alloy and undergoes important plastic deformations during loading. Internal pressure is raised up to the burst failure of the rocket motor. Geometric nonhnearities in the vessel and in the rubber joint, which undergoes very large strains (shear up to 2000/,), are very important in the behaviour of the structure 13,221. Experimental results are available for strains on the surface of the polar boss.
Fig. 7. Typical rocket motor. t
I
I
1
0.1
W,
I 0.2
I
I 0.3
.
Results. The internal pressure up to 63.35 bars was applied in only four increments, as given in Table 2. The finite element mesh used leads to 1965 DOF for a total of 271 quadrangular isoparametric axisymetric elements. Since the geometric nonhnearities were ant, the pure Ne~o~~p~ method was used ~ou~out the analysis. More details about the cahxdation may be found in [XZ]. The characteristics of the solution are summarized in Table 2. As can be seen the convergence properties obtained were very
w&tes
Fig. 6. Central deflection of an ~uminium plate.
more important and reach 5-10x. The comparison of the CPU times shows clearly that computing cost can be saved using large incremental steps. It is also interesting to compare the proposed algorithm to the metbod developsa by Bushnell in Bosor 5. For the first load step of 500 lbs, this Iast required 12 cycles of material iterations, for each of which two to three Ne~~~~n geometric iterations were needed, leading to a total of 30 geometric iterations. Up to 45 subincrements were needed in the material iterations to ensure that the stresses remain on the yield locus. With the proposed algorithm, this number was always less than 30 and decreases rapidly for successive iterations. The total number of iterations was only 7, which emphasize the economy that can result by s~ul~~usly iterating on material and geometric nonlinear&s. bc locally
0
-1-
Fig. 8. Evolution of the skin strains of the polar boss.
70
C. NYSSEN
Table 1. Aiuminium plate under concentrated load. Companson of computed solutions
izia-
Increments
..._(lb)
number .-________________
,---_s.
10
I I
2
:entra 1
deflcctfol (lo-2 in.) .____w-----_______
1
IJumber of
-3
PO times
(set
lterationl IBM 370/158 .-_______. L _____________ 4 69.
i .,,_
100.
1.849
200.
4.979
5
74.
300.
9.624
5
76.
4co.
14.16
5
77.
500.
18.01
7
98.
600.
21.52
4
61.
7cQ.
24.74
4
62.
800.
27.81
4
60.
900.
30.76
4
60.
looo.
33.66
4
60.
MO. 1000.
17.95 33.79
7 6
114. 99.
lcoo.
33.70
12.
--
L
204.
Table 2. Nonlinear analysis of a rocket motor 1ncremcnt 1 2 ________________________.___________.,___________..-----------..-----------PressJre (hrs)
Number
of
Number of elements
3
4
14
34
56
4
4
6
5
0
58
154
167
1 684
1 422
63.35
plastic
times (ICC) xnm37o/rsa
CPW
1 213
good and the total CFW time is about ten times the time of the linear analysis of the same structure. Figure 8 shows the comparison of the measured and computed strains for the polar boss. The agreement is quite satisfactory.
7. CONCLUSIONS
The Newton-Raphson method has been adapted to solve for dasto-plastic calculations. Combined with the proposed stress state determination algorithm it leads to a very ethcient iterative method. Fast convergence and accurate solutions have always been achieved for stable structures even if large incremental loads were applied. The computed solution is quite independent of the solution strategy. This allows to save both computing time and analysts time, since the loading history may be defined more easily even for an unknown structural response. Additional computational e&iency could certainly be gamed by using an automatic decision for reformation or not of the stiffness matrix at each iteration, depending on the value of the residual loads norm and its rate of decay.
REFERENCES 1. D. P. Mondkar and G. H. Powell, Evaluation of solution schemes for nonlinear structures. Compul. Srrucrures 9,
1 130
223-236 (1978). 2. J. A. Stricklin and W. E. Haisler, Formulations and solution procedures for nonlinear structural analysis. Cotnpur. Structures 7, 125-136 (1977). 3. C. Nyssen, Modelisation par elements finis du comportement non IinQire de structures airospatiales. Ph.D. dissertation, Liege (1979). 4. R. Hill, The mathematical theory of plastictty Oxford (1950). 5. D. C. Drucker, A more fundamental approach to plastic stress-strain relations. Proc.. 1st U.S. Nat. Cong. ofAppl. Mech., 487-491, Chicago (1951). 6. G. C. Navak and 0. C. Zienkiewicz, Elasto-plastic stress analysis. A generalization for various constitutive relations including strain softening. Int. J. Nwn. Meth. Engng $113-135 (1972). 7. K. J. Bathe, E. Ramm and E. L. Wilson, Finite element formulations for large deformation dynamic analysis. Inc. J. Num. Merh. Engng 9,353-386 (1975). 8. H. Jr. Armen, A. B. Pifko, H. S. Levine and G. Isakson, Plasticity. In Finite Element Techniques in S~rucl~rol Mechanics (Edited by H. Tottenham and C. rebbia), Southampton (1972). 9. 0. C. Zienkiewia, S. Valliapan and I. P. King, Elastoplastic solutions of engineering problems. “Initial stress” tinite element approach. Inr. J. Num. Merh. Engng 1,75-100 (1%9). 10. J. H. Argyris and D. W. Scharpf, Methods of elastoplastic analysis. Proc. of Symp. on Finite Element Technology at ISD, Stuttgart, Germany (1969). 11. J. Dennis and J. More, Quasi-Newton methods. motivation and theory. SIAM Rev. 19 (1977). 12. H. Matthies and G. Strang. The solution of nonhnear
An efficient and accurate iterative method to solve elasto-plastic problems
13.
14.
15.
16.
17.
18.
finite element equations. Int. J. Num. Mesh. Engng. To appear. J. A. Stricklin and W. E. Haisler. Evaluation of solution procedures for material and for geometrically nonlinear structural analysis by the direct stiffness method. AIAAJASME 13th Struct.. Structure Dyn. Material Conf. San Antonio (1972). D. Bushnell, A strategy for the solution of problems involving large ddlections, plasticity and creep. Int. J. Num. Meth. Engng 10,1343-1356 (1976). H. L. Schreyer, R. F. Kulak and J. M. Kramer. Accurate numerical solution for elasto-plastic models. J. Pressure Vessel Technol. 101.226-234 (1979). R. D. Krieg and D: B. Krieg,‘ Accuracies of numerical solution methods for the elastic-perfectly plastic model. J. Pressure Vessel Technol. 512-515 (1977). Q. S. Nguyen and J. Zarka. Quelques methodes de resolution numtrique en Clasto-plasticite classique et en Clastoviscoplasticite. Siminaire Plastic@ et fiscoplasticiti,Ecole Polytechnique, Paris (1972). D. P. Mondkar and G. H. Powell, Evaluation of state
19. 20.
21.
22.
23.
71
determination calculation in nonlinear analysis. 4th Int. Co& on Structure Mech. in React. Technol., L2/5. San Francisco (1977). NOVNL, Fiche signalttique. NOVATOME technical note MSD 78.057. S.A.M.C.E.F., Sysdme d’Analyse des Milieux Continus par Elements Finis. Theoretical Manuel. University of Liege, Belgium. J. A. Clinard, D. A. McKinley, W. C. Kroencke and W. K. Sartory, Verification by comparison of independent computer program solutions. Pressure vessels and piping computer program evaluation and qualification, ASME PVB-PB-024,27-50 (1977). H. S. Levine, H. Armen Jr., R. Winter and A. Pinko, Nonlinear behavior of shells of revolution under cyclic loading. Grumman Research Department Report RE426J (1972). C. Nyssen, Nonlinear incremental analysis up to failure of aeronautical structures. Proc. of Int. Conf. on Engng Appl. of the Finite Element Methods, p. 22. Hovik, Norway (1979).
illustrate the pefform-
to include temperature dependence of mechanical properties and kinematic hardening if necessary. If small strains are assumed, the increments of total Green strains may be split into an elastic part and a plastic part :
A large number of theoretical and computation advances in the analysis of nonlinear structures have been made in recent years. Many computer analyses have been carried out and a number of general purpose computer programs for nonlinear elastic and inelastic analyses have been developed. Nevertheless, the choice for an efficient, accurate and economical calculation algorithm to solve elastoplastic problems with a great number of degrees of freedom is still a tricky task [l-3]. It seems to be a general agreement that a good algorithm should be able to deal with large load increments, this requirement being obtained without alteration of the solution accuracy. The main dithculties in case of elastoplasticity arise from the irreversibility of the plastic strains and from the dis~ntinuous change of the material properties. From this, two distinct phases must be considered when deriving a computation algorithm in elastoplasticity : the numerical iterative method to solve the system of nonlinear equations and the numerical algorithm for scourate state dete~ination of the material. As will be shown, these two fundamental aspects are closely bounded. We propose here a particular state determination algorithm which allows to maintain good convergence properties to the Newto~~ph~n method, provided it is slightly modified. It will be shown that the numerical solution computed is quite independent of the particular calculation strategy adopted, which is also an essential requirement.
i,, =g;+ it. The elastic strains are related to the stresses by the classical linear elastic constitutive relations :
where the CiJU_are the elastic moduli. The stresses should always satisfy the plasticity criterium: F(a, f&J=&-X2=0
(3)
with ii2=f(Cij);X=XO+
H’dg,,
where f is a wnvex function which defines the plasticity criterium used. X, is the initial elasticity limit, while dB is an equivalent plastic strain. H’ is the coefficient of strain hardening. During the plastic flow, the stresses must be plastically admissible and satisfy eqn (3). Following the Drucker stability postulate [Sj. the plastic strains are given by
2. SYSTEMOF RQUATIONSTO SOLVE In this analysis, it will be assumed that the elastoplastic materials obey the Prandtl-Reuss-von Mises plastic flow rules [4]. Only isotropic materiab with isotropic hardening and at uniform temperature will be considered in the formulation for brevity. All considerations can be straightforwards generalized
The exact solution of the plasticity problem should satisfy relations (l)-(3) and (5), the equilibrium equations and the appropriate boundary conditions. During plastic flow, the increments should satisfy the consistency condition : a.Ii6.Ii =~xHBP=~~~H~=A~. 63
(6)
64
C. NYSSEN
From (1 ), (2), (5) and (61, /i may always be expressed in terms of strain increments to obtain the classical elasto-plastic incremental stress-strain relations : (7) with
one obtains K’(q)Aq=g-
(8). the I parameter takes a value of unity during plastic flow, i.e. iff(a,)= X2 and f>O: it takes a zero value otherwise. In a finite element static analysis based on the displacement formulation, the equation of equilibrium may be stated as
s
B’(q)a d I’=g, Y where a are the ~rchho~-Trefftz stresses, B is the incremental strain matrix which generally depends on the generalized displacements q, and g is a vector of generalized loads. Since the strains are assumed to be small, the stresses are related to the total strains through the constitutive equations (7), i.e. u=: r D(a)ds,
(10)
Je
where the integration is taken along the deformation path since pIasticity is essentially path dependent. The nonlinear equations to solve in an eiasto-plastic problem may be written
3. REQUIRED PROPERTKES OF AN EFFICIENT SOLUTION METHOD
3.1 Solution proceiiWes It is now generally admitted that an iterative procedure must be used to achieve sufficient accuracy in solving the nonlinear system (11) [ l-3,6-8]. The main iterative methods can be considered as particular cases of the Newton-Raphson method. Suppose that we know an approximation of the solution at iteration k, noted (a, a+41,
(12)
which satisfy the compatibility requirements and the constitutive relations (7). The exact solution may always be written as q=q+Aq
E=E-+AE
~=u+Ae,
(131
where Aq, As and Aa denote the corrections to be done on the solution (12). which are not necessary inlinitesimal. Expressing that the solution (13) must satisfy eqns (ll), it is then possible to compute for the corrections A, after linearization in these corrections. Noting that
B’fq)crdY=r,
115)
where
K’(q)=
In
s 8’
s c
[B’(q)%P3~q)+~~qw1
dV
(161
is the tangential stiffness matrix computed using the tangent material properties evaluated at the approximate state (12). In (15). r states for the out-of-balance forces. Solving (15) for Aq, a new estimation of the displacements is obtained. The corresponding state of the material is then obtained by computing the corresponding strain increments to which the stress increments are related through the elasto-plastic constitutive relations :
From the new stress state, the new out-of-balance forces may be eomputed and the procedure repeated until an appropriate norm of r becomes less than a user fixed tolerance E. It is worth noting that a solution is then reached whatever the stiffness matrix used for the iterations. Constant stiffness method. If the stiffness matrix is formed and decomposed oniy once, the constant stifEness method [9, lo] is recovered. In case of strong nonlinear&s, the method often fails to converge even if an acceleration scheme is used and this method cannot be retained solely for a general purpose computing code. Strict and modified Newton-Raphson method. In the strict Newton-Raphson method, the tangential stiffness is changed at every iteration. A disadvantage of this procedure is that a large amount of compu~~o~l effort may be required to form and decompose the stiffness matrix. But, if the convergence of the method is sufficiently fast, it can be compensated for by saving time in the state determination phase. In fact, when evaluating the computational efficiency of a solution scheme, it should be noted that the computational cost of this phase can be significant as compared with the solution phase cost [l-3]. Depending on the degree of no~in~~ty, computing time can sometimes be saved by reforming the stiffness matrix on each few iterations (modified Newtow Raphson method). It is also interesting to note that the factorized tangent stiffness matrix may also be gradually constructed at low cost while iterating using quasi Newton methods [ 11, 121. These methods seem very promising but were not investigated in this work. When only geometric nonlinear%& are present and in absence of instability, the convergence of the Newto~~p~on method is quadratic and in structural-analysis, convergence is almost always obtained in a few iterations. Due to the ability of the method to converge even in case of strong nonlinearities, the method is widely used in nonlinear structural analysis. When applied to elasto-plastic problems, some difficulties are still encountered due to the irreversibility of the plastic stress and the discontinuous variation of material properties due to eiasto-plastic state transrtion. Its use and efficiency is sometime controverted in
65
An efficrentand accurate iterative method to solve elasto-plastic problems this case [13]. The main problems of the method may be summarized as follows. The first difficulty is inherent to each type of iterative method. To each computing strategy corresponds a particular path in the load space given by B’adV
9= s”
(18)
where the u are the stresses computed during the iterations. This path can be quite sophisticated. Since an elasto-plastic solution is essentially path dependent, the value of the final solution must be discussed. Convergence difficulties may be encountered if the method is applied without special care. In case of combined material and geometric nonlinearities, the convergence of the Newton-Raphson method can excessively slow down and even be lost in particular cases [13,14]. Due to these difficulties, some modifications of the classical Newton-Raphson method should be introduced to retain the decisive advantages of reliability of the method. 3.2 Defmition of an efficient iterative scheme A solution method to be included in a general purpose nonlinear computer program should have a great reliability and offer a great flexibility in the choice of computing strategy since the best one vary with each particular problem. The computed solution should be reasonably independent of the particular strategy adopted. Finally, the method should be economic to use which implies two aspects. Firstly, in order to save engineer time, convergence should be reached whatever the reasonable choice of the load incrementation may be. Secondly, in order to save computing time, the rate of convergence should be sticiently fast to insure that the total computing cost, including the stiffness matrix reformations and the state determinations, should be acceptable. It is now well recognized [3,14,13 that these requirements can only be achieved if the computation algorithm is able to deal with large load increments, provided the number of iterations per step remains quite independent of the increment size. The state determination algorithm should consequently be very accurate, which can significantly increase the computation cost of the corresponding phase and leads to limit the number of iterations per step. It will be shown that the required properties can be reached using a Newton-Raphson type method with an appropriate state determination algorithm.
4. STATE DETERMINATION ALGORITHM
Computing of strain increments from displacement increments involves only kinematics. The problem of computing a stress increment from a given strain increment involves the material constitutive relations at two levels: i.e. how to determine the increment of stress Aa for a given strain increment ALEand how to determine the total increment of stress over a whole increment. 4.1 Stress increment associated to a given strain increment l%e stress increment during plastic flow must be
such that the plastic flow rule and the plasticity criterium should be satisfied with about the same degree of accuracy. Due to the nonlinearity of (3) and (5), the true stress increment must be computed using, t+bZ A#= Ic
D(a) ds.
(19)
Various algorithms have been designed for this purpose. 4.1.1 Single step methods. A 5rst approximation of (19) is obtained using a Euler one step forward integration scheme : AaPz D(a)Ac.
(20)
The fact that the direction of the plastic flow is only correct in the beginning of the increment can lead to important error in the final orientation of the stress vector in the stress space [15, 161. Moreover the final stress a“ will not lie on the yield surface. An additional correction is then needed and generally, the stresses are brought radially back to the yield surface. The plastic 5ow which takes place during this correction is usually completely neglected. So, the radial return seems not very consistent. To eliminate this difficulty, a single step Euler backward method can be used. In this case the flow rule is satisfied only at the end of the intaval considered while the elastic stresses are again projected radially to the yield locus. But in this implicit method, the additional plastic 5ow occuring during this return can easily be taken into account if the von Mises criterium is used [lW7]. Schreyer et al. [15] have demonstrated in a parametric study that for a single step method, the error associated with the implicit scheme is usually much less than the error obtained with the Euler forward method, but if the increments are large, the error can still be unacceptable. Moreover, the convergence properties of the iterative method singularly degrades when using the implicit scheme if the stress state is far from uniaxial. 4.1.2 Subincrementation schemes. Improvements of the accuracy of integration (19) can be achieved by multi-stepping [3, 6, 14, 16, 181. In this case, both implicit and Euler schemes lead to results that may be considered quite satisfactory for most engineering problems, for approximately the same number of subincrements [15], the errors being even less for the explicit Euler scheme in the case considered in [ 151. The algorithm implemented is based on an Euler forward multistep integration scheme since it is the most easy to use in plane stress problem and for multilinear hardening law. It is coupled with a consistent correction to bring the stresses back to the yield surface at each subincrement and is formulated as follows. Since higher derivatives are not available, a constant strain rate is assumed like in most structural programs This assumption may be violated in the true physical loading path. Strictly speaking, small enough load steps should be used. However the proposed algorithm performs very well even for quite large increments of load. When plastic loading occurs. the increment of strain AE is then divided into m equal subincrements. To estimate the parameter m, the tronmture error of a one step integration is first estimated. If a, designate
66
C. NYSSEN
the difference between the stress state reached m a single step and in a double step Euler forward integration, a measure of the troncature error for one step is taken as 2J’ ?crJ/X.
(31 I
It is then assumed that the total error is roughly I,‘m of the single step error if a m substep procedure is used. The number of subincrements is then taken equal to m =2~f’%J,)/XE,
(22)
where E is roughly equal to the requtred relative error level. Sufficient accuracy was always obtained using E=O.OS and the number m was almost always less than 30. More sophisticated algorithms could be used to determine m, based on the variation of the plastic flow direction during the integration process [14, 181, but the simple formula (22) has always proven to be sufficient in the applications. If the strain hardening law is multilinear, a strain subincrement &,,,=W,
(23)
is divided into the required number of substeps so that no discontinuity occurs within as subinterval. For each subinterval, the associated stress increment is computed by Ab,m)=W)Aa,m,,
(24)
where e states for the stresses at the beginning of the subinterval. At the end of each subincrement an additional correction is made to bring the stresses back to the yield surface. To do this, the following scheme, which is consistent with all plasticity laws to the first order, is proposed. In a displacement iterative scheme, this correction denoted 6 must be done at fixed total strains, which implies 6Ee+6&P=o.
(25)
On the other hand, we would like to have f(o+W=(X+W2,
(26)
$7=X’-/(e)+&XH”Si.
(27)
or linearizing
If due account is taken from plastic flow rule (S), it comes that (28) where H is the Hooke matrix. The corresponding correction for the stresses is given by 6a= -62Ha.
4.2 Total stress increment evaluation over an increment From strict application of the Newton-Raphson method the total strain increment at iteration i of increment n should be evaluated by
(391
As demonstrated by (28). plastic flow always occurs during this correction.
130) with at = a:, where a starred value is a value at the beginning of the considered increment. Nevertheless, if (30) is applied without care, the solution obtained may be very path dependent, and depends strongly of the particular calculation strategy adopted [3]. Moreover, fictitious numerical unloadings may occur during the iterative process. If they are considered as irreversible, these can lead to erroneous results especially in cases of combined geometric and material nonlinearities. Evenmore, convergence may be lost due to loadingunloading cycles. TO avoid this difficulty and obtain a reasonable path independent algorithm for a given load increment. many authors [l. 7, 13, 181 have proposed to compute the stress increment by integration on the total strain increment at each iteration. i.e. for iteration i of increment n, on As;= i dAs, .4=,
(31)
where dAslr is the new strain increment at iteration k. Nevertheless, this technique has a serious drawback: the cost of the state determination does not decrease during the iterative procedure since the whole strain increment is always integrated. Moreover. the strain path used to compute for (AC?): differs from the strain path used in the Newton-Raphson iteration for the approximated stress : (Ac“):,=(Aa’):-‘+D[(c”)2
dAs;.
32)
Consequently, the new state of stress may differ from the approximated value by a quantity which is not quadratic in dAe:. This technique can lead to serious convergence difficulties especially in case of combined geometric and material nonlinearities. This fact has led Bushnell [ 141 to separate in this case the geometric and material iterations to maintain the convergence of the iterative scheme. But this procedure may be costly since each material iteration requires several geometric iterations. All these drawbacks may be eliminated in another way. The path dependence of the state determination may be eliminated using an assumption of incremental reversibility. “A point which deforms plastically during one increment is assumed to unload plastically until the plastic work done become again equal to its value at the beginning of the considered increment. Then only will elastic unloading take place”. The corresponding stress-strain relation is given on Fig. 1 for a one dimensional case. Due to this assumption, the stress increment may be computed using (30). The number of integration subincrements decrease rapidly as the iterations converge, which save computing time for the state determination. Moreover, the procedure is more consistent with the Newton-Raphson method and good convergence properties are maintained by this procedure. The solution becomes also reasonably
An efficient and accurate iterative method, allowing large incremental steps, to solve elasto-plastic problems
Fig. 1. incremental reversibility assumption.
independent of the particular adopted for the increment.
computing
strategy
5. IMPROVEMENTS OF I-HE NEWTON-RAPHSON METHOD FOR ELASTGPLASTIC CASES
Two additional adaptations have been introduced in the classical Newton-Raphson method. When the loads are decreased, the plastic points undergo too large unloading when the first iteration of the corresponding increment is made using the tangent stiffness matrix. The convergence of the Newton-Raphson method may then excessively be slowed down. To avoid this difficulty, the possibility to restore the elastic stiffness matrix at the first iteration of each increment of loads was introduced. In the case of continuous loading, this additional iteration is not very time consuming since the elastic sti&ss matrix triangulation, made at the first iteration of the problem, may be saved and does not need to be rerun. Moreover, when the incremental reversibility assumption is introduced, this iteration does not generally change the solution. Secondly, non convergent loading-unloading cycles
Fig. 2. Unloading-loading cycles in a uniaxial case
67
may be observed when a plastic point unloads as far as to reload plastically during the same iteration. This phenomenon was often observed numerically in case of combined material and geometric nonlinearities. To avoid these spurious cycles, the elastic stiffness is automatically restored at this point during the increment until a new elastic state is reached during the iterative process (Fig. 2). In this way, the convergence of the Newton-Raphson method has always been restored in all the cases where the loss of convergence was due to the considered phenomenon. The final algorithm obtained has been tested on numerous examples including thermoplasticity with isotropic or kinematic hardening and combined geometric and material nonlinearities. The convergence rate reached has been always satisfactory even for very large load increments. Moreover. the number of iterations needed per step was nearly independent of the increment amplitude and the solution computed reasonably insensitive to the particular strategy adopted on the linear parts of the loading path. A few examples will now be presented to illustrate the accuracy and the efficiency of the proposed algorithm.
6. EXAMPLES
The proposed algorithm was implemented in general purpose computer programs [ 19, 201. The results presented here were obtained using either of these programs. 6.1 Tension-torsion tube Problem description. The physical problem being solved is that of a thin walled cylindrical tube which is subjected to axial tension and torsion. Since the stress distribution in the tube, remote from ends, is wnstant everywhere, the only requirements for the finite element model is that it represents material subjected to uniform tension and shear. This could be accomplished by a single plane constant stress isoparametric quadrilateral element. The material wnstants are given in Fig. 3. The non proportional loading path is specified in the same figure. The elastolplastic’nonlinearities are very severe for this case, which was retained by Clinard et al. [21] as a benchmarking problem for testing nonlinear computing codes. Results and discussion. To compute the solution, a systematic incremental history was chosen, as if the final solution was unknown. Each linear segment of the loading path was equally divided in n increments. To obtain a reference solution, n was tirstly chosen equal to 80. The solution obtamed agree with the averaged computed solution using ADINA. ANSYS, CREEPLAST, CREEPABSA and PLACRE as given in [21], within 57;. Then n was chosen equal to 5 which corresponds to quite large increments. In this case, two different computing strategies were adopted: pure Newton-Raphson iteration and modified NewtonRaphson method with a first elastic iteration at each step. Axial and shear stress strain results are plotted in Fig. 4. The results agree closely with the reference solution. When the elastic stiffness matrix is used at the tirst iteration, the shear stress-strain curves differ slightly, but the maximum difference is less than 9:$,. Figure 5 gives the convergence curves of the algorithm for a few increments and demonstrates its quality.
c. NYSSEN
Fig. 4. Stress-stun
response III the tension-torsxon
6.2 E&at&dar pixrte under cimcentrratetf lcrad droller ~~crip~io~. A simply supported ~um~ium
circular plate is submitted to a concentrated load at its center. Geometricand material data are given in Fig. 5. Isotropic hardening was adopted using a piecewise linear law with eight segments. The plate was tested by Levine et al. [22], which present aiso computed results in good agreement with the test results. but they use about 100 load steps to reach the finaf load of 1000 Ibs. Computed results using BOSOR 5 are also avaikbke in [ 141 with separate material and geometric iterations. As reported by Bushnell {i4], the stress path differs significantly from radial path near the center of the
test
plate and since b&h geometriai and material nonIinearities are s~gnifi~nt, this is a good figuration an which to verify the analysis and the strategy. Resz&s obtuined. The totai loads of 1000 Ibs was applied in ten, two and one increments respectively. The deflections of the center of the plate is given in Fig. 5, where they are compared to the other available results. Good agreement with the experimental results should be noted. Computed results and CPU times are compared in Tabie f . It is a striking fact that the differences between the central defkctions, computed in the different load incrementation cases, are less than 1”;. The differences on the final stresses computed in one or ten steps may
An efficient and accurate iterative method to solve clssto-plastic problems
Fig. 5. Convergence curves obtained in the tension-torsion
test applying large increments of load.
69
6.3 Eva&don of the burst pressureof a rocker motor Problem dewiption. The last example comes from an industrial study. It consists of a rocket motor case for which experimental results are available. The complete structure shown in Fig. 4 has three different material properties. In the first part, denoted composite material, we have anisotropic material which behaves elastically. The second part, denoted rubber joint, ensures the bounding between the vessel and the polar boss. The third part is made of aluminium alloy and undergoes important plastic deformations during loading. Internal pressure is raised up to the burst failure of the rocket motor. Geometric nonhnearities in the vessel and in the rubber joint, which undergoes very large strains (shear up to 2000/,), are very important in the behaviour of the structure 13,221. Experimental results are available for strains on the surface of the polar boss.
Fig. 7. Typical rocket motor. t
I
I
1
0.1
W,
I 0.2
I
I 0.3
.
Results. The internal pressure up to 63.35 bars was applied in only four increments, as given in Table 2. The finite element mesh used leads to 1965 DOF for a total of 271 quadrangular isoparametric axisymetric elements. Since the geometric nonhnearities were ant, the pure Ne~o~~p~ method was used ~ou~out the analysis. More details about the cahxdation may be found in [XZ]. The characteristics of the solution are summarized in Table 2. As can be seen the convergence properties obtained were very
w&tes
Fig. 6. Central deflection of an ~uminium plate.
more important and reach 5-10x. The comparison of the CPU times shows clearly that computing cost can be saved using large incremental steps. It is also interesting to compare the proposed algorithm to the metbod developsa by Bushnell in Bosor 5. For the first load step of 500 lbs, this Iast required 12 cycles of material iterations, for each of which two to three Ne~~~~n geometric iterations were needed, leading to a total of 30 geometric iterations. Up to 45 subincrements were needed in the material iterations to ensure that the stresses remain on the yield locus. With the proposed algorithm, this number was always less than 30 and decreases rapidly for successive iterations. The total number of iterations was only 7, which emphasize the economy that can result by s~ul~~usly iterating on material and geometric nonlinear&s. bc locally
0
-1-
Fig. 8. Evolution of the skin strains of the polar boss.
70
C. NYSSEN
Table 1. Aiuminium plate under concentrated load. Companson of computed solutions
izia-
Increments
..._(lb)
number .-________________
,---_s.
10
I I
2
:entra 1
deflcctfol (lo-2 in.) .____w-----_______
1
IJumber of
-3
PO times
(set
lterationl IBM 370/158 .-_______. L _____________ 4 69.
i .,,_
100.
1.849
200.
4.979
5
74.
300.
9.624
5
76.
4co.
14.16
5
77.
500.
18.01
7
98.
600.
21.52
4
61.
7cQ.
24.74
4
62.
800.
27.81
4
60.
900.
30.76
4
60.
looo.
33.66
4
60.
MO. 1000.
17.95 33.79
7 6
114. 99.
lcoo.
33.70
12.
--
L
204.
Table 2. Nonlinear analysis of a rocket motor 1ncremcnt 1 2 ________________________.___________.,___________..-----------..-----------PressJre (hrs)
Number
of
Number of elements
3
4
14
34
56
4
4
6
5
0
58
154
167
1 684
1 422
63.35
plastic
times (ICC) xnm37o/rsa
CPW
1 213
good and the total CFW time is about ten times the time of the linear analysis of the same structure. Figure 8 shows the comparison of the measured and computed strains for the polar boss. The agreement is quite satisfactory.
7. CONCLUSIONS
The Newton-Raphson method has been adapted to solve for dasto-plastic calculations. Combined with the proposed stress state determination algorithm it leads to a very ethcient iterative method. Fast convergence and accurate solutions have always been achieved for stable structures even if large incremental loads were applied. The computed solution is quite independent of the solution strategy. This allows to save both computing time and analysts time, since the loading history may be defined more easily even for an unknown structural response. Additional computational e&iency could certainly be gamed by using an automatic decision for reformation or not of the stiffness matrix at each iteration, depending on the value of the residual loads norm and its rate of decay.
REFERENCES 1. D. P. Mondkar and G. H. Powell, Evaluation of solution schemes for nonlinear structures. Compul. Srrucrures 9,
1 130
223-236 (1978). 2. J. A. Stricklin and W. E. Haisler, Formulations and solution procedures for nonlinear structural analysis. Cotnpur. Structures 7, 125-136 (1977). 3. C. Nyssen, Modelisation par elements finis du comportement non IinQire de structures airospatiales. Ph.D. dissertation, Liege (1979). 4. R. Hill, The mathematical theory of plastictty Oxford (1950). 5. D. C. Drucker, A more fundamental approach to plastic stress-strain relations. Proc.. 1st U.S. Nat. Cong. ofAppl. Mech., 487-491, Chicago (1951). 6. G. C. Navak and 0. C. Zienkiewicz, Elasto-plastic stress analysis. A generalization for various constitutive relations including strain softening. Int. J. Nwn. Meth. Engng $113-135 (1972). 7. K. J. Bathe, E. Ramm and E. L. Wilson, Finite element formulations for large deformation dynamic analysis. Inc. J. Num. Merh. Engng 9,353-386 (1975). 8. H. Jr. Armen, A. B. Pifko, H. S. Levine and G. Isakson, Plasticity. In Finite Element Techniques in S~rucl~rol Mechanics (Edited by H. Tottenham and C. rebbia), Southampton (1972). 9. 0. C. Zienkiewia, S. Valliapan and I. P. King, Elastoplastic solutions of engineering problems. “Initial stress” tinite element approach. Inr. J. Num. Merh. Engng 1,75-100 (1%9). 10. J. H. Argyris and D. W. Scharpf, Methods of elastoplastic analysis. Proc. of Symp. on Finite Element Technology at ISD, Stuttgart, Germany (1969). 11. J. Dennis and J. More, Quasi-Newton methods. motivation and theory. SIAM Rev. 19 (1977). 12. H. Matthies and G. Strang. The solution of nonhnear
An efficient and accurate iterative method to solve elasto-plastic problems
13.
14.
15.
16.
17.
18.
finite element equations. Int. J. Num. Mesh. Engng. To appear. J. A. Stricklin and W. E. Haisler. Evaluation of solution procedures for material and for geometrically nonlinear structural analysis by the direct stiffness method. AIAAJASME 13th Struct.. Structure Dyn. Material Conf. San Antonio (1972). D. Bushnell, A strategy for the solution of problems involving large ddlections, plasticity and creep. Int. J. Num. Meth. Engng 10,1343-1356 (1976). H. L. Schreyer, R. F. Kulak and J. M. Kramer. Accurate numerical solution for elasto-plastic models. J. Pressure Vessel Technol. 101.226-234 (1979). R. D. Krieg and D: B. Krieg,‘ Accuracies of numerical solution methods for the elastic-perfectly plastic model. J. Pressure Vessel Technol. 512-515 (1977). Q. S. Nguyen and J. Zarka. Quelques methodes de resolution numtrique en Clasto-plasticite classique et en Clastoviscoplasticite. Siminaire Plastic@ et fiscoplasticiti,Ecole Polytechnique, Paris (1972). D. P. Mondkar and G. H. Powell, Evaluation of state
19. 20.
21.
22.
23.
71
determination calculation in nonlinear analysis. 4th Int. Co& on Structure Mech. in React. Technol., L2/5. San Francisco (1977). NOVNL, Fiche signalttique. NOVATOME technical note MSD 78.057. S.A.M.C.E.F., Sysdme d’Analyse des Milieux Continus par Elements Finis. Theoretical Manuel. University of Liege, Belgium. J. A. Clinard, D. A. McKinley, W. C. Kroencke and W. K. Sartory, Verification by comparison of independent computer program solutions. Pressure vessels and piping computer program evaluation and qualification, ASME PVB-PB-024,27-50 (1977). H. S. Levine, H. Armen Jr., R. Winter and A. Pinko, Nonlinear behavior of shells of revolution under cyclic loading. Grumman Research Department Report RE426J (1972). C. Nyssen, Nonlinear incremental analysis up to failure of aeronautical structures. Proc. of Int. Conf. on Engng Appl. of the Finite Element Methods, p. 22. Hovik, Norway (1979).