MultiPhysicsVault/.raw/FiniteElementsinPlasticityTheoryandPractice/FiniteElementsinPlasticityTheoryandPractice_028.md

CALL ZERO(ELOAD,MELEM,MEVAB,MPOIN,MTOTG,MTOTV,NDOFN,NELEM, PLAS 39
. NEVAB,NGAUS,NSTR1,NTOTG,EPSTN,EFFST, PLAS 40
. NTOTV,NVFIX,STRSG,TDISP,TFACT, PLAS 41
. TLOAD,TREAC,MVFIX) PLAS 42
C
C*** LOOP OVER EACH INCREMENT PLAS 43
C
DO 100 IINCS = 1,NINCS PLAS 44
C
C*** READ DATA FOR CURRENT INCREMENT PLAS 45
C
CALL INCREM(ELOAD,FIXED,IINCS,MELEM,MEVAB,MITER,MTOTV, PLAS 46
. MVFIX,NDOFN,NELEM,NEVAB,NOUTP,NOFIX,NTOTV, PLAS 47
. NVFIX,PRESC,RLOAD,TFACT,TLOAD,TOLER) PLAS 48
C
C*** LOOP OVER EACH ITERATION PLAS 50
C
DO 50 IITER = 1,MITER PLAS 51
C
C*** CALL ROUTINE WHICH SELECTS SOLUTION ALORITHM VARIABLE KRESL PLAS 52
C
CALL ALGOR(FIXED,IINCS,IITER,KRESL,MTOTV,NALGO, PLAS 53
. NTOTV) PLAS 54
C
C*** CHECK WHETHER A NEW EVALUATION OF THE STIFFNESS MATRIX IS REQUIRED PLAS 55
C
IF(KRESL.EQ.1) CALL STIFFP(COORD,EPSTN,IINCS,LNODS,MATNO, PLAS 56
. MEVAB,MMATS,MPOIN,MTOTV,NELEM,NEVAB,NGAUS,NNODE, PLAS 57
. NSTRE,NSTR1,POSGP,PROPS,WEIGP,MELEM,MTOTG, PLAS 58
. STRSG,NTYPE,NCRIT) PLAS 59
C
C
C*** SOLVE EQUATIONS PLAS 60
C
C
CALL FRONT(ASDIS,ELOAD,EQRHS,EQUAT,ESTIF,FIXED,IFFIX,IINCS,IITER, PLAS 61
. GLOAD,GSTIF,LOCEL,LNODS,KRESL,MBUFA,MELEM,MEVAB,MFRON, PLAS 62
. MSTIF,MTOTV,MVFIX,NACVA,NAMEV,NDEST,NDOFN,NELEM,NEVAB, PLAS 63
. NNODE,NOFIX,NPIVO,NPOIN,NTOTV,TDISP,TLOAD,TREAC, PLAS 64
. VECRV) PLAS 65
C
C*** CALCULATE RESIDUAL FORCES PLAS 66
C
CALL RESIDU(ASDIS,COORD,EFFST,ELOAD,FACTO,IITER,LNODS, PLAS 67
. LPROP,MATNO,MELEM,MMATS,MPOIN,MTOTG,MTOTV,NDOFN, PLAS 68
. NELEM,NEVAB,NGAUS,NNODE,NSTR1,NTYPE,POSGP,PROPS, PLAS 69
. NSTRE,NCRIT,STRSG,WEIGP,TDISP,EPSTN) PLAS 70
C
C*** CHECK FOR CONVERGENCE PLAS 71
C
CALL CONVER(ELOAD,IITER,LNODS,MELEM,MEVAB,MTOTV,NCHEK,NDOFN, PLAS 72
. NELEM,NEVAB,NNODE,NTOTV,PVALU,STFOR,TLOAD,TOFOR,TOLER) PLAS 73
C
C*** OUTPUT RESULTS IF REQUIRED PLAS 74
IF(IITER.EQ.1.AND.NOUTP(1).GT.0) PLAS 75
.CALL OUTPUT(IITER,MTOTG,MTOTV,MVFIX,NELEM,NGAUS,NOFIX,NOUTP, PLAS 76
. NPOIN,NVFIX,STRSG,TDISP,TREAC,EPSTN,NTYPE,NCHEK) PLAS 77
C
C*** IF SOLUTION HAS CONVERGED STOP ITERATING AND OUTPUT RESULTS PLAS 78
C
IF(NCHEK.EQ.0) GO TO 75 PLAS 79
50 CONTINUE PLAS 80
C
C*** PLAS 81
C
PLAS 82
PLAS 83
PLAS 84
PLAS 85
PLAS 86
PLAS 87
PLAS 88
PLAS 89
PLAS 90
PLAS 91
PLAS 92
PLAS 93
PLAS 94
PLAS 95
PLAS 96
PLAS 97
PLAS 98
PLAS 99
PLAS 100
PLAS 101
PLAS 102
PLAS 103 

IF(NALGO.EQ.2) GO TO 75	PLAS 104
STOP	PLAS 105
75 CALL OUTPUT(IITER,MTOTG,MTOTV,MVFIX,NELEM,NGAUS,NOFIX,NOUTP,NPOIN,NVFIX,STRSG,TDISP,TREAC,EPSTN,NTYPE,NCHEK)	PLAS 106
100 CONTINUE	PLAS 107
STOP	PLAS 108
END	PLAS 109
	PLAS 110

7.9 Numerical examples

The first numerical example considered is illustrated in Fig. 7.12(a). The problem studied is that of a thick cylinder subjected to a gradually increasing internal pressure, with plane strain conditions being assumed in the axial direction. A Von Mises yield criterion is assumed and the numerical solutions obtained compared with the theoretical results of Reference 14. The pressure/radial displacement characteristics are shown in Fig. 7.12(b) and good

radar

Angle (°)	Value
40	50
45	50
50	50
55	50
60	50
65	50
70	50
75	50
80	50
85	50
90	50
95	50
100	50
105	50
110	50
115	50
120	50
125	50
130	50
135	50
140	50
145	50
150	50
155	50
160	50
165	50
170	50
175	50
180	50
185	50
190	50
195	50
200	50
205	50
210	50
215	50
220	50
225	50
230	50
235	50
240	50
245	50
250	50
255	50
260	50
265	50
270	50
275	50
280	50
285	50
290	50
295	50
300	50
305	50
310	50
315	50
320	50
325	50
330	50
335	50
340	50
345	50
350	50
355	50
360	50
365	50
370	50
375	50
380	50
385	50
390	50
395	50
400	50
405	50
410	50
415	50
420	50
425	50
430	50
435	50
440	50
445	50
450	50
455	50
460	50
465	50
470	50
475	50
480	50
485	50
490	50
495	50
500	50
505	50
510	50
515	50
520	50
525	50
530	50
535	50
540	50
545	50
550	50
555	50
560	50
565	50
570	50
575	50
580	50
585	50
590	50
595	50
600	50
605	50
610	50
615	50
620	50
625	50
630	50
635	50
640	50
645	50
650	50
655	50
660	50
665	50
670	50
675	50
680	50
685	50
690	50
695	50
700	50
705	50
710	50
715	50
720	50
725	50
730	50
735	50
740	50
745	50
750	50
755	50
760	50
765	50
770	50
775	50
780	50
785	50
790	50
795	50
800	50
805	50
810	50
815	50
820	50
825	50
830	50
835	50
840	50
845	50
850	50
855	50
860	50
865	50
870	50
875	50
880	50
885	50
890	50
895	50
900	50
905	50
910	50
915	50
920	50
925	50
930	50
935	50
940	50
945	50
950	50
955	50
960	50
965	50
970	50
975	50
980	50
985	50
990	50
995	50
1000	50

Elastic modulus, E = 2.1 \times 10^{4} \, dN/mm^{2}

Poissons ratio, \nu = 0.3

Uniaxial yield stress, \sigma_{r}=24.0\ dN/mm^{2}

Strain hardening parameter, \mathbf{H}' = 00

Von mises yield criterion

(a)

line

displacement of inner face (× 10⁵ mm)	3 point gauss rule	2 point gauss rule	theory (Ref. 14)
0	0	0	0
8	8	8	8
12	12	12	12
16	16	16	16
20	18	18	18
24	19	19	19
28	19	19	19

Fig. 7.12 (a) Mesh and material properties employed in the elasto-plastic analysis of an internally pressurised thick cylinder under plane strain conditions. (b) Displacement of the inner surface with increasing pressure for the problem of Fig. 7.12(a).

agreement between the numerical and analytical solutions is evident. In the numerical studies, collapse was deemed to have occurred if the iterative procedure diverged for an incremental load increase.

line

x	σθ
100	13
110	12
120	11
130	10
140	9
150	8
160	7
170	6.5
180	6
190	5.5
200	5

line

x	σₙ
100	16
120	15
140	13
160	11
180	9
200	8

line

x	σ_e
100	13
110	16
120	19
130	18
140	16
150	14
160	12
170	11
180	10
190	9
200	8

line

x	σθ
100	10
120	15
140	20
160	22
180	20
200	18

text_image

Two rows of circular symbols, likely representing data points or indicators, with dots inside each row.

text_image

P 2 3 4 5 6 7 8 9 10

Fig. 7.13 Hoop stress distributions at various pressure values for the problem of Fig. 7.12(a).

Fig. 7.13 shows the circumferential (hoop) stress distributions for specified pressure values. Again a good agreement is evident. In solution both a two-point and three-point Gaussian integration rule was considered. Whilst the nodal displacements obtained by use of both rules are practically identical, it is seen from Fig. 7.13 that use of a 2 \times 2 integrating rule gives superior stress values to a 3 \times 3 rule. This is a general result for elasto-plastic problems and therefore use of a two-point rule is recommended. This phenomenon is an example of the benefit of a reduced integration order for parabolic isoparametric elements. ^{(15)}

line

Central deflection, w	Intensity of UDL, P
0.00	0
0.05	100
0.10	150
0.15	200
0.20	250
0.25	275
0.30	280

Fig. 7.14 Load/central deflection response for a uniformly loaded simply supported circular plate.

The second example considered is the simply supported circular plate shown in Fig. 7.14.

The plate is modelled by five axisymmetric elements and the loading takes the form of a progressively increasing uniformly distributed load. The growth in central deflection with increasing load is shown in Fig. 7.14. A converged solution was obtained for P = 270 but the numerical process diverged for P = 280 and consequently the collapse load is taken to be 270. This is in good agreement with the value of 260 quoted in Ref. 16, particularly in view of the coarse mesh employed in the present study. Fig. 7.15 shows the deflection profile with increasing applied load.

line

r/a	w/h (P = 100)	w/h (P = 200)	w/h (P = 260 lb/in²)
0.0	0.05	0.15	0.20
0.2	0.05	0.15	0.20
0.4	0.05	0.15	0.20
0.6	0.05	0.15	0.20
0.8	0.05	0.15	0.20
1.0	0.05	0.15	0.20

Fig. 7.15 Deflection profiles for the problem of Fig. 7.14 at various applied load values.

7.10 Problems

7.1 In Section 7.2.1 it was stated that the Von Mises law implies that yielding begins when the (recoverable) elastic energy of distortion, D , reaches a critical value. Prove this by showing that J_2' is proportional to D , since D can be written as

D = \frac {1}{2} \sigma_ {i j} \epsilon_ {i j} - \frac {(1 - 2 \nu)}{1 2 \mu (1 + \nu)} (\sigma_ {i i}) ^ {2}. \tag {7.98}

text_image

σ₃ S √3P₁ 0 σ₂ σ₁ √3P₀

Fig. 7.16 Geometric representation of the Berg yield criterion—Problem 7.2.

7.2 A yield criterion has been proposed by Berg ^{(17)} which attempts to account for the tensile failure of a material due to the formation of voids at a sufficiently high strain level. The yield surface is illustrated in Fig. 7.16 and can be seen to be made up of two distinct portions. For stress levels below a mean hydrostatic tension of P_{I} the material yields according to the Von Mises cylinder of radius S. The yield surface in the tensile range is terminated by an elliptic cap whose extremity is defined by P_{0} . The three constants S, P_{I} and P_{0} are material constants and must be experimentally determined. The two distinct portions of the yield surface can be expressed as

\sqrt {2} (J _ {2} ^ {\prime}) ^ {\ddagger} = S \quad \text { for } \sigma_ {m} \leqslant P _ {I}

[ 2 J _ {2} ^ {\prime} + H (\sigma_ {m} - P _ {I}) ^ {2} ] ^ {\frac {1}{2}} = S \quad P _ {I} \leqslant \sigma_ {m} \leqslant P _ {0}, \tag {7.99}

where H = S^2 / (P_I - P_0)^2 and \sigma_m is the mean hydrostatic pressure.

Show that this yield criterion can be expressed in the form of three constants C_{1} , C_{2} and C_{3} as indicated in Section 7.4 where

C _ {1} = 0, \quad C _ {2} = \sqrt {2}, \quad C _ {3} = 0 \quad \text { for } \quad \sigma_ {m} \leqslant P _ {I}

C _ {1} = H \left(\sigma_ {m} - P _ {I}\right) / S, \quad C _ {2} = 2 J _ {2} ^ {\prime} / S, \quad C _ {3} = 0 \quad P _ {I} \leqslant \sigma_ {m} \leqslant P _ {0}.

7.3 A certain material yields when the maximum principal stress reaches a critical value, Y. Assuming identical behaviour in tension and compression, determine the geometrical form of the yield surface. The solution is given in Fig. 7.17.

text_image

σ₃ 30° √2/3 Y σ₁ σ₂

Fig. 7.17 \pi plane representation of a yield criterion based on maximum principal stress values—Problem 7.3.

7.4 The assumption of a linear strain hardening material law may prove to be inadequate for certain situations. If the uniaxial stress/strain test curve for the material is known, then it is possible to represent the stress-plastic strain relationship in a piecewise linear fashion as shown in Fig. 7.18 and the instantaneous yield stress can be written in the form \sigma_{Y} = \sigma_{Y^{0}} + S(\bar{\epsilon}_{p}) where S(\bar{\epsilon}_{p}) is the piecewise linear function describing the increase (or decrease) in the initial yield stress \sigma_{Y^{0}} with the increase of effective plastic strain \bar{\epsilon}_{p} . The program modifications required to describe this behaviour will all be included in subroutine RESIDU, except for changes in material property specification which will need to be made in subroutine INPUT. Carry out all necessary modifications.

line

Effective plastic strain ε̅p	Stress, σ
Δε̅p (1)	σ⁰_γ
Δε̅p (2)	σ²_γ
n-1 (3)	σ²_γ
n (n)	σ²_γ

Fig. 7.18 Piecewise-linear representation of material strain hardening—Problem 7.4.

7.5 By using the mesh of Fig. 7.12(a) and solving as an axisymmetric problem, use program PLANET (documented in Appendix II, Section A2.1) to determine the elasto-plastic stress and displacement distributions in a sphere when it is loaded by an incrementally applied internal pressure. The dimensions and material properties of the sphere are given by reference to Fig. 7.12. Assume a Tresca yield criterion for solution and compare your results with the solution given in Ref. 1.
7.6 Use program PLANET to solve the problem illustrated in Fig. 1.2, Chapter 1. Use both a Tresca and Von Mises yield criterion and compare the plastic zone distributions obtained with those of Fig. 1.2.
7.7 Subroutine CONVER, described in Section 6.5.4, bases convergence of the nonlinear solution process on the global norm of the residual force vector. Modify subroutine CONVER so that convergence is based on expression (3.27) in which the summation signs are absent; so that convergence is monitored locally at each of the nodes 1 to N in turn.
7.8 Modify subroutine CONVER, Section 6.5.4 so that convergence is monitored locally at each node according to the displacement changes that occur during a particular iteration, r, as follows.

\frac {| \Delta d ^ {r} |}{| d ^ {1} |} \times 1 0 0 \leqslant \text { TOLER }, \tag {7.100}

where d^{1} is the elastic displacement occurring upon application of the load increment and \Delta d^{r} is the change in nodal displacement during the r^{th} iteration.

7.9 Modify program PLANET to undertake the elasto-plastic solution of three-dimensional solids. To simplify the task consider only the Von Mises yield criterion and assume that the solid is loaded by nodal point loads only.
7.10 The yield criterion to be employed in program PLANET is specified by means of control parameter NCRIT in subroutine INPUT described in Section 6.5.1. In some applications, such as steel-concrete composites, it is necessary to employ a different yield surface for different parts of the structure. Modify program PLANET so that the yield criterion governing elasto-plastic behaviour is separately specified for each element in the solid.

7.11 References

1. HILL, R., The Mathematical Theory of Plasticity, Oxford University Press, 1950.
2. PRAGER, W., An Introduction to Plasticity, Addison-Wesley, Amsterdam and London, 1959.
3. HOFFMAN, O. and SACHS, G., Introduction to the Theory of Plasticity for Engineers, McGraw-Hill, 1953.
4. BRIDGMAN, P. W., Studies in Large Plastic Flow and Fracture, McGraw-Hill, New York, 1952.
5. BISHOP, A. W., The strength of soils as engineering materials, Geotechnique, 16, 89–130 (1966).
6. DAVIS, E. M., Theories of plasticity and the failure of soil masses, Ch. 6 Soil Mechanics, Ed. I. K. Lee, Butterworths, London, 1968.
7. ZIENKIEWICZ, O. C., VALLIAPPAN, S. and KING, I. P., Elasto-plastic solutions of engineering problems; Initial stress finite element approach, Int. J. Num. Meth. Engng. 1, 75–100 (1969).
8. YAMADA, Y., YOSHIMURA, N. and SAKURAI, T., Plastic stress-strain matrix and its application for the solution of elastic-plastic problems by Finite Element Method, Int. J. Mech. Sci. 10, 343–354 (1968).
9. BLAND, D. R., The associated flow rule of plasticity, J. Mech. Phy. of Solids, 6, 71–78 (1957).
10. NAYAK, G. C. and ZIENKIEWICZ, O. C., Convenient form of stress invariants for Plasticity, Journ. of the Struct. Div. Proc. of A.S.C.E., 949-953, April 1972.
11. FREDERICK, D. and CHANG, T. S., Continuum Mechanics, Allyn and Bacon, 1965.
12. KOITER, W. T., Stress-strain relations, uniqueness and variational theorems for elastic-plastic materials with singular yield surface, Quart. Appl. Math., 11, 350–354 (1953).
13. HINTON, E. and OWEN, D. R. J., Finite Element Programming, Academic Press, London, 1977.
14. HODGE, P. G. and WHITE, G. N., A quantitative comparison of flow and deformation theories of plasticity' J. Appl. Mech. 17, 180–184 (1950).
15. ZIENKIEWICZ, O. C. and HINTON, E., Reduced integration, function smoothing and non-conformity in finite element analysis (with special reference to thick plates), J. of the Franklin Institute, 302, Nos. 5 and 6, Nov./Dec. 1976.
16. ARMEN, H., Jr., PIFKO, A. and LEVINE, H. S., Finite element analysis of structures in the plastic range, N.A.S.A. Contractor Report, CR.1649 (1971).

17. BERG, C. A., Plastic dilation and void interaction, Battelle Inst. Material Science Colloquia, Sept. 1969. Inelastic Behaviour of Solids, Ed. M. F. Kanninen et al., McGraw-Hill, 171–210, 1970.
18. KRIEG, R. D. and KRIEG, D. B., Accuracies of numerical solution methods for the elastic-perfectly plastic mode, Trans. ASME, J. Pressure Vessel Technology, 99, 4, 510–515 (1977).
19. SCHREYER, H. L., KULAK, R. F. and KRAMER, J. M., Accurate numerical solutions for elastic-plastic models, Trans. ASME, J. Pressure Vessel Technology, 101, 3, 226–234 (1979).
20. SAMUELSSON, A. and FRÖIER, M., Finite elements in plasticity—A variational inequality approach, in The Mathematics of Finite Elements and Applications III, Ed. J. R. Whiteman, Academic Press, 1979.