Mathematical Simulations of Deformation for the Rotation Shells with Variable Wall Thickness

Well-posed boundary value problems are constructed for calculating rotation shells of with a stiffness variable along the meridian in two directions, and also with variable bilateral with respect to the reference surface with the shell wall thickness. Algorithms for the numerical integration of systems of differential equations with variable coefficients are discussed.


Introduction
Shell structures are widely used in the creation of structures for modern mechanical engineering, in the oil and gas, chemical and other industries. At the same time, the requirements for ensuring the strength reliability while reducing the weight indicators lead to the need to build more and more reliable models and methods for calculating shell structures with non-homogeneous parameters (in particular, with variable stiffness) [4, 6-8, 10, 13-15, 17].
In the problems of determining the optimal distribution of material [9] or calculating the durability of shells taking into account the degradation of their surface in an aggressive environment [16], the stiffness parameters change at each step of successive approximations. This leads to the necessity of restructuring the grid at each step of the corresponding iterative computational (search) algorithm using the known finite element analysis packages [1, 5].

Basic Relations
The basic equations of the moment theory of shells are obtained under the assumption that shells of revolution (with an arbitrary meridian shape in the general case) and circular (annular) plates are homogeneous, isotropic, thin and elastic. The validity of Kirchhoff's hypotheses is accepted, as well as the smallness of deformations and angles of rotation in comparison with unity. The shell wall thickness is generally considered arbitrary h = h (s, ϕ).
Deformation of the middle surface ε 1 , ε 2 , γ 12 the angles of rotation of the normal to the surface ϑ 1 , ϑ 2 and the parameters of the change in curvature χ 1 , χ 2 , χ 12 are associated with the displacements u, v, w (Fig. 1)  where θ(s) is the angle between the normal to the median surface and the shell rotation axis; r(s) is the radius of the parallel circle. Elastic ratios are taken in the usual form where D = Eh 3 /(12(1 − µ 2 )) is the cylindrical stiffness; R 1 , R 2 are thr radii of curvature of the surface; E, µ are the modulus of elastic and Poisson's ratio, respectively. As a result, the equations of forces and moments can be written as follows: Here, for the force factors, the generally accepted designations are introduced [3]: N 1 , N 2 , S, M 1 , M 2 , M, Q 1 , Q 2 , Р°q 1 , q 2 , q 3 are the meridional, circumferential and normal components of the intensity of the external load, respectively. As the main variables with respect to which the system is written, four quantities u, v, w, ϑ 1 are selected, which characterize the displacements and the four force factors N 1 , S * , Q * 1 , M 1 corresponding to them, where S * = S + 2M R 2 ; Q * 1 = Q 1 + 1 r ∂M ∂ϕ are reduced efforts. After appropriate transformations, the equations of the moment theory for homogeneous isotropic elastic thin-walled shells of variable thickness under asymmetric loading can be, as is known [3,12], reduced to a system of eight partial differential equations, which can be written in the form:

Shells with Stiffness, Variable in two Directions
To solve system (2.4) in the case h = h (s, ϕ), it is proposed to use the method of straight lines [2,4,7,12], the essence of which is to replace the derivatives in the direction of the circumferential coordinate by difference relations, which allows one to obtain a set of one-dimensional boundary value problems along each ith meridian (i = 1, m), which are subsequently solved by the sweep method with orthogonalization according to S.K. Godunov [11] along the nodal points j (j = 1, n).
Assuming that the shell state parameters are sufficiently smooth in the circumferential direction, the partial derivatives with respect to the variable ϕ = ϕ i , (i = 1, 2, ..., m) equations (2.4) are replaced by finite differences of the fourth order of accuracy [2]: where ∆ = ∆ϕ i is the step of the difference grid in the circumferential direction with the approximation error O ∆ 4 . It should be noted that, depending on the method of replacing derivatives with finite differences, different systems of the method of lines are possible, which may differ in the accuracy of the approximating relations for the corresponding derivatives.
After replacing the derivatives with respect to the coordinate ϕ in the system of equations (2.4) and reducing similar terms for the components of the stressstrain state vector Y = (u, v, w, ϑ 1 , N 1 r, S * 1 r, Q * 1 r, M 1 r) T for each i-th strip (i = 1, 2, ..., m), we can obtain a boundary value problem for a system of ordinary differential equations.
Further (in order to avoid cumbersome presentation) such a system for the case of a cylindrical shell is given (sin θ = 1, cos θ = 0, R 1 = ∞, R 2 = R): Here, As for the fulfillment of the boundary conditions, for the case of an open shell at the ends of the variation interval ϕ, one-sided differences are used, where for the straight lines 0 and m, the values of the parameters specified in accordance with the conditions for fixing the contour are taken into account. For a closed cylindrical shell, when the derivatives in the circumferential direction are replaced by their finite-difference expressions, only the central differences are used. The boundary conditions at the meridional edges of the shell (s = 0, s = L where L is the shell length) are taken into account when solving boundary value problems along the meridian s = s i , (i = 1, 2, ..., n), using the sweep method with orthogonalization according to S.K. Godunov [11].

Rotational Shells with a Stiffness Variable Along Meridian
In the case h = h (s) of separating variables, the Fourier method is used. At the same time, the problem of solving the system of partial differential equations (2.4) by expanding the components of the stress-strain state and load into trigonometric series in the circular coordinate [3,10,12,15], reduces, in the general case, to solving systems of t+1 ordinary differential equations for finding harmonics of expansions of the sought functions in Fourier series.
The decomposition of the load, displacements and forces acting in the shell into Fourier series [3,15] along the circumferential coordinate ϕ is carried out in the form where f , in the generally accepted notation [3], stands for the functions u, w, ε 1 , ε 2 , ϑ 1 , χ 1 , χ 2 , N 1 , N 2 , Q 1 , M 1 , M 2 , q 1 , q 3 , whereas ψ can be substituted by functions v, γ 12 , ϑ 2 , χ 12 , S, Q 2 , M, q 2 , f c k , f t k , ψ c k , ψ t k , and f c k , f t k , ψ c k , ψ t k are the coefficients of their expansions in trigonometric series.
With this choice of functions, the expansion coefficients with the superscript t, which correspond to the skew-symmetric deformation of the shell with respect to the zero meridian, are determined by exactly the same system of equations as the coefficients with the index c, which correspond to the symmetric deformation. Therefore, the results of further transformations for these coefficients coincide, which allows them to be carried out only for functions with the index c, omitting this sign.
In this case, the displacements and forces that correspond to the k-th term of the expansion are determined by the formulas [3] u = u k (s) cos kϕ; v = v k (s) sin kϕ; w = w k (s) cos kϕ; ϑ 1 = ϑ 1k (s) cos kϕ; ϑ 2 = ϑ 2k (s) sin kϕ; N 1 = N 1k (s) cos kϕ; The use of the Fourier method (and this is possible only in the case when the shell wall thickness changes only in the meridional direction h = h(s), and remains constant in the circumferential direction) makes it possible to reduce the adopted system of equations of state of the shell in partial derivatives to a system of ordinary differential equations with respect to the coefficients expansions in trigonometric series of the main variables of the stress-strain state, which are the coefficients of the expansion of displacements and force factors. In this case, it is convenient to take for the main unknowns their product by the radius r of the parallel circle N 1k r(s), S * k r(s), Q * 1k r(s), M 1k r(s)) : The expansion coefficients of displacements and forces for each harmonic number k, which are not the main variables, using the relations of the theory of elasticity and the dependences between displacements and deformations, are expressed in terms of the main variables as follows: The inconvenience of the system of equations (4.3) is that the forces and displacements are related to the local coordinate system associated with the normal and tangent to the shell meridian. Therefore, the coefficients of the system have discontinuities when the shell meridian consists of several sections with corner points between them. In this case, it turns out to be necessary to draw up compatibility equations for different areas.
These difficulties, in accordance with [3], can be circumvented by passing to global coordinates. For this, forces and displacements are projected onto the normal to the shell's symmetry axis and onto the axis itself. In this case, instead of displacements u, w, displacements ξ, ζ are introduced, and instead of forces N 1 , Q * 1 , X, Z forces are introduced as follows: The same dependencies are also related to the coefficients of the expansion in the Fourier series of the corresponding functions. Substitution of u k , w k , N 1k , Q * 1k and their derivatives through ξ k , ζ k , X k , Z k into system (4.3) brings it to the form: where the following designations are introduced for the radial and axial loading components q xk = q 1k cos θ + q 3k sin θ; q zk = q 1k sin θ − q 3k cos θ. Since the coefficients of the resulting system of equations do not contain the curvature 1/R 1 of the meridian, they remain continuous even for a shell whose curvature is subject to discontinuity. As a consequence, the main unknowns, referred to a fixed coordinate system, remain continuous for an arbitrary shape of the meridian, including for combined shell, which makes it possible not to compose the docking equations for such cases. As for the force unknowns X k r, Z k r, S * k r, M 1k r, they experience discontinuity of a predetermined magnitude only where concentrated forces are applied to the shells at a specific parallel of the load.
The specified forces and displacements at the ends of the shell are the boundary conditions for the resulting system, which are also decomposed into the corresponding trigonometric series along the circumferential coordinate.

Shells with Two-Sided Relative to the Reference Surface Change in Wall Thickness
Let us consider the case of the shell wall thickness δ(s) variable along the meridian under axisymmetric loading. In this case, both the external and internal components of the shell wall thickness with respect to the reference surface are independent functions of the meridional coordinate s. Let us denote by H(s) and h(s) the distances along the normal direction from the reference surface to the outer and inner surfaces of the shell, respectively, and h * (s) the distance from the middle surface to the reduced surface ( Fig. 5.2), so that It is also assumed that the error associated with the mismatch of the normal to the middle surface and the reference surface can be neglected for thin-walled shells.
Taking into account that the deformations ε 1 , ε 2 and parameters of changing the curvatures χ 1 , χ 2 of the reference surface are expressed through the deformations of the middle surface as follows ε 1 = ε 1cp − h * χ 1 , ε 2 = ε 2cp − h * χ 2 (Fig.  5.2), the relationship between stresses and deformations on the reference surface in accordance with Hooke's law will have the form Note that expressions (5.2) at z = h * coincide with the well-known [3] expressions for calculating the stresses on the middle surface.
Internal forces and moments relative to the reference surface, taking into account the known dependencies and relations (5.2) after introducing the notation will be next: Considering that N 1 sin θ −Q 1 cos θ = F (s); N 1 cos θ +Q 1 sin θ = N , (Fig. 2.1) we represent the efforts N 1 and Q 1 in the form where N is the spacer force; F (s) = P 0 +ˆs n s 0 (q n cos θ − q 1 sin θ) 2πrds is total axial loading; P 0 is axial load; q n , q 1 is distributed normal and meridian loads, respectively. We take as the main variables the radial displacement ξ, the angle of rotation of the normal ϑ, the axial displacement ζ, as well as the spacer force N r and moment M 1 r multiplied by the radius of the parallel circle.
Eliminating from (5.5) taking into account we obtain Similarly, relations (5.6) yield the expression for M 2 From equation (5.5), taking into account (5.7), (5.8), we obtain Substituting further (5.7), (5.8), (5.11) into (5.6), we obtain from which one of the equations of the system of state follows cos θ r N r + (5.12) Substitution of the found expressions (5.5), (5.7)-(5.9) into equation (5.13) gives one more equation of the system The following equation is obtained from (5.14) after substituting into it the values of the variables from (5.7)-(5.9), taking into account that d ds cos θ r ϑ (5.17) +µ cos θ r (N r) + µ sin θ r F (s) 2π − q r r where q r = q 1 cos θ + q n sin θ. Substituting into (5.15) the expression for and from (5.6), (5.7), we obtain d (M 1 r) dr = K 2 1 − µ 2 cos θ r ξ + D 1 − µ 2 cos 2 θ r ϑ + sin θ (N r) + µ cos θ r (M 1 r) − cos θ F (s) 2π . (5.18) To determine the axial displacement, the expression for (5.11) is substituted into the equation of continuity of deformations, which will have the form Thus, the obtained equations (5.14) , (5.12), (5.17) -(5.19) form a system of differential equations with variable coefficients, which describes the stress-strain state of shells of revolution with a two-sided, relative to the reference surface, change along the meridian wall thickness.
For the particular case, when the reference surface coincides with the median, that is h * (s) = 0, it follows that H(s) = h(s). Then (5.4) will have the form K 1 = Eδ/ 1 − µ 2 , K 2 = 0, D = Eδ 3 / 12 1 − µ 2 , where δ (s) = H (s) + h (s) is the shell thickness, and the system of the obtained equations coincides with the well-known system of equations given in [3].
The boundary conditions are the conditions for fixing the ends of the shell. For the numerical solution of the obtained boundary value problem for a system of ordinary differential equations with variable (due to a change in the shell wall thickness) coefficients under the given boundary conditions, a sufficiently effective and repeatedly tested in the problems of mechanics of thin-walled structures [3,12,15] are used the sweep method of S. K. Godunov [11].

Conclusion.
The article presents correct mathematical models describing the state of asymmetrically loaded shells of revolution with variable wall thickness in the meridian and circumferential directions, only along the meridian, as well as with a two-sided change in the wall thickness relative to the reference surface. For all considered cases, boundary value problems for systems of ordinary differential equations with variable coefficients are constructed, the numerical solution of which is carried out by the sweep method.