On the basis of physical equations of deformation theory of plasticity using Kirchhoff–Lava hypothesis, matrix dependences between columns of forces and moments and columns of deformations and curvatures of the shell midface are determined at the loading step. As a finite element, a quadrilateral fragment of the shell midface with nodal unknowns in the form of: increments of displacements and their derivatives; increments of deformations and increments of curvatures; increments of forces and increments of moments were used. To approximate the required quantities, the following expressions are adopted: bicubic functions with elements of Hermite polynomials of the third degree for displacements; bilinear functions for deformation and force parameters. To obtain the stiffness matrix of the finite element, the nonlinear Lagrangian functional on the loading step was used with an additional condition: the real work of the difference of forces determined using their direct approximation and using approximating expressions for displacements, on deformations and curvatures of the loading step must be equal to zero. Minimisation of the functional by nodal unknowns provides three systems of equations, the solution of which determines the stiffness matrix of the finite element used to calculate the displacement fields. The force and deformation parameters at the discretisation nodes of the shell are determined from the displacements found. Case studies show the effectiveness of using a three-field finite element method (FEM) technique compared to using FEM in the displacement method formulation (single-field technique).
1.Introduction
Definition of the object of study. Structures consisting of shells and their fragments are now increasingly used in various fields. These include hangars, warehouses, domes and slabs, tanks, bunkers, cisterns, pipelines and others. The worldwide trend aimed at significant reduction of material intensity of systems and objects actualises the problem of determining the stress-strain state (SSS) of structures in the form of shells and their fragments in a physically nonlinear formulation.
Assumption within the regulated limits of the plastic stage of the applied material of structures allows to reduce the overall material intensity of structures, including those consisting of shells.
Literature review. The developed theory of solid body deformation [1–4] turned out to be theoretically unrealisable in the practice of engineering calculations of real structures, which led to the development of numerical methods for solving the equations of solid body deformation mechanics [5–9]. At the present stage of development of structural mechanics and computer science, the main tool for investigating the SSS of shell structures beyond the elastic limit is numerical finite element method (FEM). In spite of the considerable volume of publications on this subject [10–15] and the availability of foreign (ANSYS, ABAQUS, NASTRAN, etc.) and domestic (PRIIS, LIRA, etc.) finite element computational systems, the problem of finding the most optimal formulations of the FEM for the calculation of shell structures in a physically nonlinear setting remains quite relevant. It is known that the FEM in the form of the displacement method is the most widespread at present, but it does not lack the disadvantages associated with the use of only one field of unknowns – the displacement field. The mixed variant of the FEM with kinematic and force fields of unknowns [16–23], as well as the FEM variant with three fields of unknowns [24, 25], are becoming the most promising.
Purpose and objectives of the study. In this paper, the three-field variant of the mixed FEM is applied to the calculation of shells under elastic-plastic deformation and the finite element solutions obtained on its basis are compared with the solutions obtained using the FEM in the formulation of the displacement method in a physically nonlinear setting.
2.Materials and Methods
In order to realise the three-field variant of the mixed FEM in a physically nonlinear formulation, it is first of all necessary to formulate the corresponding nonlinear three-field functional. For this purpose, we introduce the following matrix notations: is matrix-string of longitudinal forces and bending moments and their increments at the loading step, respectively; is matrix-string of deformation increments and curvature increments at the midpoint of the shell structure surface at the loading step; is matrix-string of step increments of displacement vector components; are matrix-string of external surface load vector components and their step increments, respectively.
The increments of longitudinal forces and increments of bending moments at the loading step can be expressed through the increments of deformation and increments of curvature using the following relations:
(1)
(2)
where and are matrices obtained by numerically finding definite integrals; are the increments of the contravariant components of the stress tensor at the loading step.
The plasticity matrix included in (1) and (2) is composed on the basis of the relationships of the deformation theory of plasticity [4], represented in a curvilinear coordinate system by the expression:
(3)
where are components of the deformation deviator; are intensity of deformations and stresses; is stress deviator; are first invariants of strain and stress tensors; are covariant and contravariant components of the metric tensor.
The increments of deformations in an arbitrary layer at a loading step are determined by differentiating (3) in the following general form:
(4)
The partial derivatives included in (4) were determined by the relations:
(5)
where is modulus of elasticity; is coefficient of transverse deformation.
The derivatives of the ratio of strain and stress intensities included in (5) are determined by the expressions:
(6)
where are tangent and secant modules of the deformation diagram;
Using (4), (5), (6) a matrix dependence is formed:
(7)
where
Matrix determined by matrix inversion
Taking into account the above, the matrix relationship is written:
(8)
where
Taking into account (1) and (2), the following matrix relation can be composed as follows:
(9)
where
If we use a four-node fragment of its midface [24] with nodes as a discretisation element of the shell structure, the increments of longitudinal forces and increments of bending moments at the loading step can be expressed through their nodal values by means of bilinear relationships:
(10)
where are local coordinates used to organise the procedure of numerical integration by Gauss quadrature;
Here, is understood as or On the basis of (10), a matrix dependence can be compiled:
(11)
where is a quasi-diagonal matrix, on the main diagonal of which the row matrices
Taking into account the notations introduced above and (1)–(11), the nonlinear mixed functional can be written in the following form:
(12)
The column of increments of deformations and curvatures of the medial surface of the shell structure included in (12) can be expressed similarly to (11) through their nodal values:
(13)
and can be represented by the Cauchy relations for thin shells [26] by a matrix product:
(14)
The column of step increments of the components of the displacement vector of a point of the centre surface of the shell structure can be interpolated through the nodal values of the component increments by means of products of Hermite polynomials of the third degree:
(15)
where is a quasi-diagonal matrix containing matrix-rows of polynomial Hermite functions; and are columns of kinematic nodal unknowns at the loading step in local and global coordinate systems, respectively, and is the transition matrix from column to column
The relation (14) taking into account (15) will take the form:
(16)
where
Here, means or and and are curvilinear global coordinates.
The functional (12) taking into account (8)–(11) and (13)–(16) can be transformed to the form:
(17)
For convenience of further calculations, we introduce the following matrix notations:
(18)
The functional (17) taking into account (18) can be written in the following form:
(19)
By successively minimising (19) by and we can obtain the following system of matrix equations:
(20)
To solve the system (20), we can use the substitution method. For this purpose, from the first and second equations (20), it is necessary to express the force and deformation step nodal unknowns:
(21)
By substituting (21) into the third equation (20), the following matrix equation can be obtained:
(22)
or in a more convenient form:
(23)
where is the stiffness matrix of the used four-node discretisation element of the three-field variant of the mixed FEM at one of the successive loading steps; is the column of step forces; is the Newton–Raphson correction at the loading step.
To obtain the stiffness matrix of the finite element at the loading step, a numerical integration procedure was applied using the Gauss method for the area of the mid-surface and the Simpson formula when integrating over the shell thickness.
The global stiffness matrix of the shell structure is composed of a four-node discretisation element by means of an index matrix formed according to the accepted boundary conditions of the calculated shell [27].
In order to verify the above algorithm for the calculation of shell structures in a physically nonlinear formulation, a comparative analysis of finite element solutions obtained using the developed variant of the mixed FEM with the solutions obtained on the basis of the FEM in the formulation of the displacement method was performed.
3. Results and Discussion
Example 1. As an example, a fragment of an elliptical cylinder made of duralumin alloy with the ratio of ellipse parameters of cross-section = 5, loaded with internal pressure of intensity = 6·10–3 MPa, was calculated. The calculation scheme of the shell is shown in Fig. 1. The following initial data were used: = 1.5 m; = 0.3 m; = 0.01 m; = 0.01 m; = 7.49·104 MPa; = 0.32. The deformation diagram was assumed as a two-linked broken line defined by the formula:
where
Figure 1. Calculation diagram of an elliptical cylinder.
The calculations were performed in two variants: in the first variant, the above three-field mixed FEM algorithm in a physically nonlinear formulation was implemented; in the second variant, the FEM algorithm in the formulation of the displacement method was used. The results of the variant calculations are presented in tabular form. The tables show the numerical values of normal stresses in the rigid embedment and at the free end of the shell, as well as the values of bending moment (for the first variant of calculation) when varying the degree of refinement of the discretisation grid and the number of loading steps.
Tables 1–3 present the results of the first variant of calculation at successive densification of the mesh of nodes 41×2, 51×2, and 61×2 depending on the number of stages of successive loading.
The selected design scheme allows to calculate the values of physical bending moment in the rigid embedment (Fig. 2):
Figure 2. Section of the shell by the plane ZOY.
It is also obvious that the bending moment at the free end of the elliptical cylinder, as well as the stresses, must be equal to zero.
The rightmost column of Tables 1–3 shows the above-mentioned analytical values of the controlled strength parameters of the shell SSS.
Table 1. Values of controlled SSS parameters at 41×2 node grid.
Point coordinates, y, m; z, m
Stress σ, MPa, moment, M22, N·m
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
347.3
347.5
347.6
347.6
–
-346.7
–347.1
–347.1
–347.1
–
70.74
70.55
70.22
70.16
70.20
1.5;
0.0
–0.192
–0.190
–0.190
–0.191
–
–0.137
–0.136
–0.136
–0.136
–
0.0687
0.0681
0.0681
0.0685
0.000
0.0005
0.0005
0.0005
0.0005
0.000
Table 2. Values of controlled SSS parameters at 51×2 node grid.
Point coordinates, y, m; z, m
Stress σ, MPa, moment, M22, N·m
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
347.4
347.6
347.6
347.6
–
–346.7
–347.1
-347.2
–347.2
–
70.87
70.36
70.32
70.29
70.20
1.5;
0.0
–0.0708
–0.0682
–0.0679
–0.0683
–
–0.0533
–0.0514
–0.0511
–0.0515
–
0.0273
0.0266
0.0265
0.0266
0.000
0.0004
0.0004
0.0004
0.0004
0.000
Table 3. Values of controlled SSS parameters at 61×2 node grid.
Point coordinates, y, m; z, m
Stress σ, MPa, moment, M22, N·m
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
347.5
347.6
347.7
347.7
–
–346.7
–347.1
–347.2
–347.2
–
71.03
70.50
70.37
70.33
70.20
1.5;
0.0
–0.0295
–0.0298
–0.0299
–0.0297
–
–0.0223
–0.0225
–0.0225
–0.0224
–
0.0119
0.0120
0.0120
0.0120
0.000
0.0003
0.0003
0.0003
0.0003
0.000
As follows from the analysis of the data given in Tables 1–3, the three-field version of the mixed FEM in the physically nonlinear formulation demonstrates stable convergence of the computational process, both when the discretisation grid is reduced and when the number of successive loading stages is increased. In addition, the numerical values of the bending moment in the rigid embedment and at the free end practically coincide with their analytical values, which is also a proof of the correctness and high degree of accuracy of finite element solutions obtained by using the developed three-field mixed FEM in the physically nonlinear formulation.
If the classical FEM formulation of the displacement method is applied to the solution of this problem, the results of finite element solutions will be quite different from the above-mentioned ones. For example, Table 4 shows the results of the FEM calculation of an elliptical cylinder in the form of the displacement method with a 61×2 node grid.
Table 4. Stress values of elliptical cylinder FEM in the form of displacement method with 61×2 node grid.
Point coordinates, y, m; z, m
Stress σ22, MPa
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
304.1
274.6
120.3
276.0
–
–303.8
–274.3
–112.7
–275.8
–
1.5;
0.0
–76.44
–66.17
–140.8
–62.41
–
358.4
323.2
330.2
319.4
–
270.7
250.8
246.4
248.9
0.000
As follows from the analysis of Table 4, there is no convergence of the computational process when the number of successive loading steps is increased. In addition, the stresses at the free end of the shell reach unacceptably high values, although they should be equal to zero.
Obviously, when using the FEM in the formulation of the displacement method in physically nonlinear calculations of shells with significant curvature of the medial surface, a much more significant refinement of the discretisation grid is required. Tables 5–7 show the results of elliptical cylinder calculations with 81×2, 101×2, and 121×2 node grids, respectively.
Table 5. Stress values of elliptical cylinder FEM in the form of displacement method with 81×2 node grid.
Point coordinates, y, m; z, m
Stress σ22, MPa
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
343.9
302.9
369.0
264.1
–
–343.4
–302.2
–370.4
–263.0
–
1.5;
0.0
–66.74
–59.54
–32.91
–56.99
–
278.5
259.7
253.6
248.7
–
213.8
178.8
176.7
156.3
0.000
Table 6. Stress values of elliptical cylinder FEM in the form of displacement method with 101×2 node grid.
Point coordinates, y, m; z, m
Stress σ22, MPa
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
345.0
345.2
344.2
341.6
–
–344.5
–344.8
–343.8
–341.2
–
1.5;
0.0
–47.34
–47.49
–47.21
–46.41
–
227.5
227.7
227.2
225.6
–
114.9
115.2
114.6
112.9
0.000
Table 7. Stress values of elliptical cylinder FEM in the form of displacement method with 121×2 node grid.
Point coordinates, y, m; z, m
Stress σ22, MPa
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
346.3
346.6
346.6
346.6
–
–345.8
–346.1
–346.1
–346.1
–
1.5;
0.0
–37.30
–37.41
–37.43
–37.44
–
145.1
145.6
145.7
145.8
–
69.49
69.74
69.80
69.81
0.000
Analysing the tabulated values of normal stresses in Table 5 shows that for 81×2 node grid, there is also no convergence of the computational process as the number of sequential loading steps increases and the stresses at the free end have unacceptably high values.
For 101×2 node grid in rigid embedment (Table 6), satisfactory convergence of the computational process is observed. However, the stresses at the free end have unacceptably high values.
With a 121×2 node grid (Table 7), the values of stresses in the rigid embedment remain almost unchanged and coincide with the values of stresses in the rigid embedment obtained using the three-field version of mixed FEM (Tables 1–3). However, at the free end, the normal stresses still remain very far from zero values, which requires further refinement of the discretisation grid.
Table 8 shows the calculation results of the elliptical cylinder at 201×2 node grid. As can be seen from the analysis of the tabular data, a steady convergence of the computational process is observed in the rigid termination as the number of successive loading steps increases. At the free end of the shell, the stress values have decreased but still remain quite far from zero values.
Table 8. Stress values of elliptical cylinder FEM in the form of displacement method with 201×2 node grid.
Point coordinates, y, m; z, m
Stress σ22, MPa
Number of loading steps
Analytical solution
22
52
82
102
0.0;
0.3
347.2
347.5
347.6
347.6
–
–346.8
–347.0
–347.1
–347.1
–
1.5;
0.0
–17.44
–17.49
–17.50
–17.51
–
40.28
40.41
40.44
40.45
–
16.31
16.37
16.38
16.38
0.000
In addition, it should be noted that when using the traditional formulation of the FEM in the form of the displacement method, the problem of obtaining the numerical value of the bending moment in the rigid embedment arises. To solve this problem, a graphical method can be applied to calculate the bending moment value from the normal stress epiphysis in the rigid embedment, plotted at fixed points of vertical coordinate along the normal to the medial surface of the shell.
Fig. 3 shows the normal stress epiphysis in the rigid embedment with a 201×2 node grid plotted at 9 node points along the surface normal. For convenience of further calculations, the epyure was divided into elementary figures and the centres of gravity in each figure were found. Then, the forces were calculated as the areas of each of the figures. Moments were calculated as products of forces by the corresponding arms
Figure 3. Normal stress diagram in a rigid termination.
The resultant moment can be obtained by summing the moments :
After performing the above calculations, the resultant torque is obtained as follows:
The calculation error was:
Obviously, when building a more detailed stress diagram, for example, by dividing it by height into 17 points, the calculation error can be reduced. However, the error of 0.92 % obtained by dividing the stress diagram into 9 points along the section height is quite acceptable for engineering calculations.
Calculation example 2: A fragment of an elliptical ring with the ratio of ellipse parameters = 6, loaded on the right side with a linear load of intensity = 25 kN/m uniformly distributed along the formations and having a hinged support on the left side, was calculated (Fig. 4). The following initial data were used: = 1.2 m; = 0.2 m; = 0.008 m; = 0.01 m. The physical characteristics of the material and the deformation diagram were taken from the previous calculation example. A reactive force equal to the applied nodal load occurs at the hinge points Thus, the normal stresses at points and must be equal. The equality of the normal stresses at points and can serve as an additional criterion for the correctness of the numerical values of the normal stresses obtained as a result of the solution.
Figure 4. Calculation diagram of an elliptical ring.
The calculations, as in the previous example, were performed in two variants: the first variant used the developed three-field mixed four-node discretisation element; the second variant used a finite element whose stiffness matrix was composed in the formulation of the displacement method. The results of the first variant of the calculation are summarised in Tables 9 and 10. Table 9 shows the values of normal stresses on the inner and outer surfaces of the shell at the points of application of a given load and at the points of hinge support depending on the density of the discretisation grid at a fixed number of loading steps equal to 22.
Table 9. Stress values of the first calculation variant depending on the size of the discretisation grid.
Point coordinates, y, m; t, m
Stress σ22, MPa
Sampling grid