Mixed FEM for Shells of Revolution Based on Flow Theory and its Modifications

Full Text

1. Introduction For the majority of deformable materials, Hook’s law is only valid at loading levels, at which the stresses do not exceed the yield stress of the material. Usually, plastic deformations emerge in stress concentration zones already at insignificant levels of loading. Hence, structural analysis with account of elastoplastic deformation zones is an important engineering problem. Two elastoplastic deformation theories are most commonly used for solid bodies: flow plasticity theory and the theory of incremental elastoplastic deformation11[1-3]. Displacement-based finite element method (FEM) has been widely used for elastoplastic deformation analysis2 [4-7]. This method was applied to thermoplastic and contact problems of continuum mechanics [8-12]. FEM was also effectively employed in finite strain cases of elastoplastic deformation processes [13-16]. Mixed finite element method has been extensively applied to elastoplastic deformation problems [17-21]. In this study, a prismatic finite element with triangular bases has been developed in mixed FEM formulation. Three versions of governing equations are used as constitutive relations. The first version uses the flow theory equations. The second version employs the governing equations obtained from the authors’ hypothesis, that the components of the incremental plastic strain tensor are proportional to the components of the combined stress deviator tensor. The third version does not separate strains into the elastic and plastic parts. For determining the relationships between the strain increments and the stress increments, the condition of proportionality between the components of the incremental strain deviator tensor to the components of the incremental stress deviator tensor was used. 2. Methods 2.1. Shell Geometry Arbitrary point М0t of the shell, which is located at distance t from the middle surface, is defined by the following position vector: R R ar0t r0 tr0, (1) s r r r where R0 xi r sinθ j r cosθk is the position vector of the corresponding point r М0 of the middle r r surface; r is the radius of curvature of the middle surface point; i j k, , are the unit vectors of the Cartesian coordinate system; x,θ are the axial and angular coordinates of point М0; ar0 a ar10 r10 is the normal line to the middle surface at point М0; a ar10 r10 are the unit basis vectors at point М0. The basis vectors of arbitrary point М0t are determined by differentiating position vector (1): r r r r r r g R g R10 ,0st; 02 ,0rtθ; g R30 ,0tt ar0, (2) and by following [17], the matrix expressions of the derivatives of the basis vectors of an arbitrary point in the basis of this point are formed: gr,0s m gr0 ; gr,0r n gr0 ; gr,t0 l gr0 , (3) 3 1 3 3 3 1 3 1 3 3 3 1 3 1 3 3 3 1 where gr,0s T g gr r r,10s ,20s g,0s3 ; gr,0rθ T gr,0r1θ gr,0r 2θ gr,0r 3θ ; gr,0t T gr,t01 gr,t02 gr,t03 ; 1 3 1 3 1 3 gr0 T gr10 g gr20r30 are the row matrices of the derivatives of the basis vectors of М0t. 1 3 Under gradually applied load, the incremental displacement vector at a load step is represented by components in the basis of point М0t: r 1 0r v g2 0r2 v g3 0r3 gr0 T v , (4) V v g1 1 3 3 1 where v T v1 v2 v3 is the row matrix of displacements of point М0t. 1 3 . 2024;20(1):27-39 The derivatives of the displacement vector are also expressed in terms of the basis vectors of point М0t: 1 1g f1 2 VVVrrr ,,,s ff1 0f1 0gr1 0rgr f f2 02 0g2 0rgrgr ff3 31 33 0f3 02 3g ,r3 0g ;rg ;r (5) r 2 1 2 2 t 3 1 3 2 where f11 v,1s v m1 11 v m2 21 v m3 31;L f33 v,3t v l1 13 v l2 23 v l3 33; mij,nij,lij are the elements of matrices [m], [n] and [l]. Under specified load, an arbitrary point of the shell will displace to position Мt, which is determined by position vector r r r R R Vt 0t . (6) The strain increments for a load step are governed by relations [3] in a geometrically linear definition ij 1 gri0 vr,j gr0j vr,i . (7) 2 Considering (5), strains (7) can be expressed in matrix form as L v , (8) 6 1 6 3 3 1 Т ss tt 2 s 2 st 2 t is the row matrix of strain increments; L is the matrix where 1 6 of differentiation operators. 2.2. Relations of Flow Plasticity Theory Full strain increments ij are combinations of elastic strains ije and plastic strains ijp : ij ije ijp . (9) The relationships between the elastic strain increments and the stress increments are defined by expressions[4] ije 1 1 ij c ij , (10) E where Е is the material Young’s modulus; is the Poisson’s ratio; c is the mean value of the normal stress increments; ij is the Kronecker delta. In the flow theory, the plastic strain increments are defined by relations[5] ijp k ij c ij , (11) 6 where k is the coefficient of proportionality, which is defined according to expression k 3 i E1k E1н i . (12) 2 Here: i is the stress intensity; Eн is the modulus of the initial segment of the stress-strain intensity diagram; Eк is the tangent modulus at the considered point on the stress-strain intensity diagram; i mn is the stress intensity increment. i mn By combining (10) and (11) and taking into account (12), the matrix expression for the constitutive equations of the flow theory is formed: ε С1П σ . (13) 6 1 6 6 6 1 The second version of the post-yield constitutive equations uses the hypothesis of proportionality between the components of the incremental plastic strain tensor and the components of the incremental stress deviator tensor: ijp 1 ij c ij . (14) Proportionality coefficient 1 is defined according to[6] expression i 32 Е1k Е1н . (15) By combining (10) and (14), the matrix expression for the second version of constitutive relations is obtained: ε С2П σ . (16) 6 1 6 6 6 1 The third version of constitutive relations is based on the hypothesis of proportionality between the incremental strain deviator tensor and the incremental stress deviator tensor components: ij ij c 2( ij ij c) , (17) where ψ2 3 εi 3 1 , and the volumetric strain increment is determined as in the case of elastic 2 σi 2 Еk deformation, c c 1 2 . E Based on (17), the third version of the constitutive relations is formed: ε С3П σ . (18) . 2024;20(1):27-39 2.3. Finite Element Stiffness Matrix A prismatic finite element with triangular bases is considered. The nodal unknowns are the displacement and stress increments. Coordinates s, ,θ t of an arbitrary point of the shell are defined in terms of nodal coordinates using linear functions ξ, ,ηζ with ranges 0≤ξ,η≤ 1; -1≤ζ≤1, λ={f (ξηζ, , )}T { }λy , (19) 1 6× 6 1× where { }λ y T = λ λ λ λ λ λ{ i j k m n p} is the row of nodal coordinate s,θ or t ; {f (ξ1 6×,η,ζ)}T = (1- -ξ η)1-ζMξ1-2ζMη1-2ζM(1- -ξ η)1+2ζMξ1+2ζMη1+2ζ 2 By using linear approximating functions (19), the interpolation expressions for Δv components and the components of the incremental stress tensor are formed: { }Δ =v [ ]А {Δvy}; { }Δσ = [ ]S {Δσy}, (20) 31× 318× 181× 6 1× 6 36× 361× where { }Δvy Т = Δ Δ Δ Δ Δ Δ{ v1i v1j v1k v2i v2 j v2k... Δ Δ Δv3m v3n v3p} is the row-matrix of the nodal displace- 118× ment increments; { }Δσ = Δσ Δσ Δσ Δσ Δσ ΔσТ { ss θθ tt sθ st θt} is the row of stress increments at a point; 1 6× {Δσy}Т = Δσ { ssy} {T Δσθθy}T L{Δσθty}T is the row of stress increments at the nodes of the finite 1 36× 1 6× 1 6× 1 6× element. { }Δε =[ ]L { }Δ =v [ ] [ ]L A Δvy = [ ]B Δvy . 6 1× 6 3× 31× 6 3 318× × 181× 6 18× 181× The nonlinear mixed functional for a load step, obtained in [17], is expressed as (21) Considering (20), strain increments (8) can be represented in matrix form: { } { } Ф V 1 6× 6 3× 3 1× V 1 6× 6 6× 3 1× - 12 S{Δ1 3×v} {T Δ3 1×q dS} - Δ S{1 3×v} { }T 3 1q dS× +V { } {σ ε1 6× T Δ6 1× }dV; (μ=1,2,3 .) (22) 2 Taking into account matrix relations (18) and (21), functional (22) for the prismatic finite element becomes HEORY OF Ф у T S T B dV v y 12 y T S T C П S dV y 1 36 V 36 6 6 18 18 1 1 36 V 36 6 6 6 6 36 36 1 1 vy T A T q dS vy T A T q dS vy T B T dV. (23) 2 1 18 S 18 3 3 1 1 18 S 18 1 3 1 1 18 V 18 6 6 1 By varying functional (23) with respect to nodal unknowns y T and vy T , the following systems of 1 36 1 18 equations are obtained: Q v y T 36 18 18 1y 36 H 36 36 1y 0; vy T 18 Q 36T 36 1 18f q1 18 R 1 0, where Q S T B dV ; H S T С П S dV; fq A T q dS; 36 18 V 36 6 6 18 36 36 V 36 6 6 6 6 36 18 1 S18 3 3 1 R A T q dS B T dV is the Raphson residual. 18 1 S 18 3 3 1 V 18 6 6 1 Systems (24) can be combined into one (24) K Zу Fy , (25) Ф Ф 54 54 54 1 54 1 H Q where K Q36 36T 36 0 18 - is the matrix of the stress-strain state of the hybrid finite element at a load 54 54 18 36 18 18 step; Zу Т y T vy T - is the vector of nodal unknowns; Fу Т 0 T M fq T R T - 1 54 1 36 1 18 1 54 1 36 1 18 1 18 is the vector of nodal loads with residuals. 3.1. Sample Calculation 1 The shell of revolution depicted in Figure 1 with the middle surface in the shape of a truncated ellipsoid was analyzed. The following input values were specified: а = 0.21 m; в = 0.15 m; h = 0.01 m; lk 0.2 m; Е = 2×105 MPa; = 0.3. The height of truncation of the elliptical shell is zk в 1 lak22 0.15 1 0.200.2122 0.0457 m. The stress-strain curve for the shell material was assumed to be in the form of Figure 2, where σТ 200 MPa is the yield stress of the material; εТ 0.001 is the yield strain; εk 0.02 is the final strain; σk 400 MPa is the final stress. Figure 1. Truncated elliptical shell Sou rce: made by R.Z. Kiseleva The stress-strain intensity curve was constructed using formulas[7] σi 1 σ11 σ22 2 σ22 σ33 2 σ33 σ11 2 1 σ2 0 σ2 σ; εi 3 ε11 ε22 2 ε22 ε33 2 ε33 ε11 2 3 ε νε 2 0 νε ε 2 2 1 3 ν ε. Figure 2. Stress-strain curve of the elliptical shell material Figure 3. Stress-strain intensity curve of the shell material Sou rce: made by R.Z. Kiseleva Sou rce: made by R.Z. Kiseleva Values of the parameters of the stress-strain intensity curve: iТ Т 200 MPa is the stress intensity at yield point; εiТ 2 1 ν ε Т 2 1 0.3 0.001 = 0.866667×10-3 is the strain intensity at yield point; 3 3 ik 2 1 k 2 1 0.3 0.01 = 0.866667×10-2 is the final strain intensity; 3 3 iк к = 300 MPa is the final stress intensity. The stress-strain intensity curve is assumed to be defined by function σi f εi in the form of a parabola σi aεi2 bεi c (when i iТ ), where а = -6612835.5282 MPa; b = 242231.47902 MPa; с = 1795.0330258 MPa. The presented shell of revolution was analyzed for the case of elastic deformation (q = 18.0 MPa). The normal stress values at the fixed support are presented in Table 1, where the first column contains the number of discretization nodes of the shell along its axis (NM) and along its thickness (NT). The other columns contain normal stresses of the internal fibers along the axis int and circumference 11 int . For the external fibers, these variables are denotes as σ11ext and σθθext respectively. The results presented in Table 1 demonstrate convergence of the computational process with respect to normal stresses of the shell at the fixed support. 3.2. Sample Calculation 2 The analysis of the shell from the previous section was performed under internal pressure q = 27.65 MPa. The specified load value was achieved in 16 steps and in 32 steps, and the results of the analysis using the three versions of constitutive equations were found to be virtually identical. The values of meridional stresses ss and circumferential stresses after 32 load steps are presented in Table 2. The stress values were calculated along the shell thickness h in the left section using the third version of the constitutive equations. Table 2 Numerical values of meridional and circumferential stresses after 32 load steps along the shell thickness h in the left section The of results from Table 2 are used to plot the distributions of meridional stresses (Figure 4) and circumferential stresses (Figure 5). In order to control the accuracy of computation of meridional stresses, the check of х 0 is performed. The check gives an acceptable discrepancy in the values of the resultant external and internal forces: Qext-Qint δ= 100% 2.4% , Qext where Qext is the resultant external force; Qint is the resultant internal force. Figure 4. Distribution of meridional stresses ss Figure 5. Distribution of circumferential stresses σθθ along the section thickness Sou rce: made by R.Z. Kiseleva Sou rce: made by R.Z. Kiseleva As seen from Figure 5, the circumferential stresses exceed the elastic limit significantly. Table 3 provides the values of meridional and circumferential stresses in the external fibers along the meridian arc length. Table 3 Numerical values of meridional and circumferential stresses in external fibers along the meridian arc length The results from Table 3 were used to plot the distributions of meridional stresses ss and circumferential stresses (Figure 6). Figure 6. Distributions of meridional and circumferential stresses in external fibers along the meridian arc length Sou rce: made by R.Z. Kiseleva HEORY OF The values of the meridional stresses in the end section are almost zero, which complies with the loading condition. The circumferential stresses vary insignificantly along the meridian. The analysis of the results in Tables 1-3 indicates correctness of the developed algorithm and shows adequate convergence of the computational process. 4. Conclusion The 3D stress-strain state of a shell is studied without using the straight-normal hypothesis for elastoplastic deformation. 1. The constitutive relations beyond the elastic limit are implemented in three versions. The first version uses the relationships of the flow theory. The second version employs the governing equations, where the authors’ hypothesis is used for determining the plastic strain increments. The hypothesis assumes that the components of the incremental plastic strain tensor are proportional to the components of the incremental stress deviator tensor. The third version of equations is based on the hypothesis of proportionality between the components of the incremental strain deviator tensor and the components of the incremental stress deviator tensor without separating the strains into elastic and plastic. 2. The analysis of the shell is performed using mixed FEM. For this purpose, the authors developed a 6node solid prismatic finite element with triangular bases. The nodal unknowns are the displacement vector components and the nodal stress tensor components. The target variables are approximated by the nodal unknowns using bilinear shape functions. 3. The presented study shows that all three versions of the governing equations for plastic deformation produce identical results. The analysis of the constitutive equations shows that the most physically reasonable version is the third one. This version does not separate the strain increments into elastic and plastic parts, and is based on the hypothesis of proportionality between the components of the incremental strain deviator tensor to the components of the incremental stress deviator tensor. The proposed governing equations, without the strain separation, correspond to the physical meaning of the process of deformation and have great potential for analyzing reservoirs, submersibles and other engineering structures containing shells of revolution.

About the authors

Rumia Z. Kiseleva

Volgоgrad State Agrarian University

Author for correspondence.
Email: rumia1970@yandex.ru
ORCID iD: 0000-0002-3047-5256

Candidate of Technical Sciences, Associate Professor of the Department of Applied Geodesy, Environmental Management and Water Management

Volgograd, Russia

Natalia A. Kirsanova

Financial University under the Government of the Russian Federation

Email: nagureeve@fa.ru
ORCID iD: 0000-0003-3496-2008

Doctor of Physical and Mathematical Sciences, Professor of the Department of Mathematics

Moscow, Russia

Anatoliy P. Nikolaev

Volgоgrad State Agrarian University

Email: anpetr40@yandex.ru
ORCID iD: 0000-0002-7098-5998

Doctor of Technical Sciences, Professor of the Department of Mechanics

Volgograd, Russia

Yuriy V. Klochkov

Volgоgrad State Agrarian University

Email: klotchkov@bk.ru
ORCID iD: 0000-0002-1027-1811

Doctor of Technical Sciences, Professor, Head of the Department of Higher Mathematics

Volgograd, Russia

Vitaliy V. Ryabukha

Volgоgrad State Agrarian University

Email: vitalik30090@mail.ru
ORCID iD: 0000-0002-7394-8885

Postgraduate student of the Department of Mechanics

Volgograd, Russia

References