Numerical Analysis of Deformation Characteristics of Elastic Inhomogeneous Rotational Shells at Arbitrary Displacements and Rotation Angles

: Adequate mathematical models and computational algorithms are developed in this study to investigate speciﬁc features of the deformation processes of elastic rotational shells at large displacements and arbitrary rotation angles of the normal line. A ﬁnite difference method (FDM) is used to discretize the original continuum problem in spatial variables, replacing the differential operators with a second-order ﬁnite difference approximation. The computational algorithm for solving the nonlinear boundary value problem is based on a quasi-dynamic form of the ascertainment method with the construction of an explicit two-layer time-difference scheme of second-order accuracy. The inﬂuence of physical and mechanical characteristics of isotropic and composite materials on the deformation features of elastic spherical shells under the action of surface loading of “tracking” type is investigated. The results of the studies conducted have shown that the physical and mechanical characteristics of isotropic and composite materials signiﬁcantly affect the nature of the deformation of the clamped spherical shell in both the subcritical and post-critical domains. The developed mathematical models and computational algorithms can be applied in the future to study shells of rotation made of hyperelastic (non-linearly elastic) materials and soft shells.


Introduction
There are many papers devoted to the study of shells, but there are still a number of unresolved issues that require attention [1][2][3]. Elastic single and multi-layered shells of rotation under the action of an axisymmetric system of boundary or surface static loads of general or local nature are considered. For a single-layer structure, the mid-surface of the shell is taken as the coordinate surface (datum surface); for a multi-layer structure, the mid-surface of one of the layers or the contact surface of the layers can be taken as the coordinate surface. To describe the deformation processes of elastic rotational shells at large displacements of coordinate surface points and unlimited rotation angles of the normal to it, the changes of Cartesian coordinates x, y are taken as unknowns, which allows for the considered variant of geometrically nonlinear deformation of thin-walled constructions to greatly simplify the structure of initial differential equations in comparison with using traditional components of tangential u and normal w displacements [4,5].
Unlike most traditional mathematical models for the study of shell deformation processes, when the geometric parameters of the shell are identical both in the initial and deformed state (for example, the curvature and the parameters of the Lamet before and after deformation are the same), this work takes into account the difference in the geometry of the shell in the initial and final (after loading-genia) state. This makes it possible to study the features of deformation of elastic multi-layer shells in the region of large displacements and angles of rotation of the normal. The developed mathematical models and computational algorithms can be practically implemented in various fields of science and technology, for example, in aerospace systems (deployable elastic antennas, displacement elements of fuel systems) and in the automotive industry when designing automotive-mobile springs made of multi-layer composites, which undergo significant shape changes during operation, etc. These issues are partially considered in the works [6,7].

Materials and Methods
Before the deformation, an element of an unstrained formative shell of length ds 0 has initial coordinates x 0 , y 0 , curvatures k 10 , k 20, and initial angle θ 0 between the x-axis and the normal to the generatrix (Figure 1). Given the coordinates x 0 , y 0 as linear coordinate functions along the generatrix x 0 = x 0 (s 0 ), y 0 = y 0 (s 0 ), for the initial (undeformed) state, the following relations take place (1): where R 10 is the initial radius of curvature of the surface in the direction of the formant; R 20 is the initial radius of curvature in the circumferential direction ( Figure 1). After deformation, the element will have length ds, coordinates x, y, and curvatures k 1 , k 2 . Considering the coordinates x, y as functions of the linear coordinate s along the deformed generatrix x = x(s), y = y(s), it is possible to write down relations similar to (1) for the deformed state (2): where θ is the rotation angle between the normal to the deformed coordinate surface and the x-axis.
There are no restrictions imposed on rotation angles ∆θ = θ − θ 0 . The deformation components along the generatrix E 11 and in the circumferential direction E 22 are defined as (3): The curvature change parameters along the generatrix K 11 and in the circumferential direction K 22 are determined from the obvious ratios (4): The deformation components at distance z from the coordinate surface, taking into account the hypotheses and assumptions made for thin shells, are distributed according to a linear law (5): We introduce force factors in the cross-section of the shell: longitudinal force T 11 and circumferential force T 22 , transverse force Q 13 , as well as bending moments M 11 and M 22 : M 11 -in longitudinal direction, M 22 -in circumferential direction. A distributed load q = q(s) with components: q 1 = q 1 (s) and q 2 = q 2 (s) acts on the shell locally or over the entire surface. For the conservative load case q 1 = q x (s), q 2 = q y (s), and for the "tracking" load case q 1 = q u (s), q 2 = q w (s). The positive directions for the force factors T 11 , T 22 , Q 13 , M 11 , M 22 and load components q 1 , q 2 are shown in Figure 2. In the case of a single-layer construction of an orthotropic material, the forces and moments in the elastic shell are expressed in terms of the deformed state components according to the following dependencies (6), (7) [2,3]: where where E 1 , E 2 are the Young's moduli in the longitudinal and circumferential directions, respectively; ν 12 , ν 21 are the Poisson's ratios: the first index indicates the direction of deformation and the second index indicates the direction of force.
The following correlation between the elastic characteristics of the material takes place (8): For a multi-layer structure as a whole, the conditions of rigid contact between layers without mutual separation and slippage are assumed, whereby layers of variable thickness h m are generally considered; m is the layer index: 1 ≤ m ≤ M, M is the number of layers ( Figure 3). The force factors in the multi-layer package are expressed in terms of the deformation components of the coordinate surface according to the Formulas (9), (10) [8]: where Formula (10) is obtained by assuming the constancy of the elastic characteristics of the 21 within the m-layer [9][10][11]. The equilibrium equations of the deformed element of a rotational shell have the following form (11) [4,5,8]: The system of differential equations of the developed mathematical model (1)-(11) is supplemented by boundary conditions at the edges s = s 0 and s = s L . The characteristic boundary conditions at the edge s = s 0 can be written in the following way: • rigid fixing (12): • pin-edge fixing (13): The boundary conditions at the edge s = s L are formulated similarly to (12)-(13).

Results and Discussion
When constructing the discrete analog of the original continuum problem (1)-(13), FDM [3,4] is used to discretize the spatial variables. Two grids are introduced in the area of continuous variation of the argument s: a main grid with integer indices i and an auxiliary grid with indices i ± 1/2, the nodes of which are positioned in the middle between the nodes of the main grid i ± 1 (Figure 4). The partial derivatives in differential relations (1)-(13) are approximated by difference operators of second-order accuracy (∆s 2 ). The undeformed state of the shell is given by the grid functions (x 0 ) i , (y 0 ) i . Stress-strain state (SSS) parameters (1)-(10) are related to main and auxiliary grid nodes. Finite-difference approximations of differential relations (1) for the initial state of the rotational shell can be represented in the following way (14), (15): where In the relations (14), (15) the grid functions (x 0 ) i , (y 0 ) i are set from the condition of constancy of the step ∆s 0 = const for the whole computational domain of the shell 0 ≤ s 0 ≤ L 0 , so (16): where N is the number of sampling points: The deformed state parameters of the shell are described by the grid functions of nodal displacements x i and y i with approximation of relations (2) by finite difference operators (17), (18) similar to (14), (15) (Figure 4): where The grid functions of the strain and curvature change components in (3) and (4) are approximated at the corresponding nodes of the main and auxiliary grid as follows (19): As follows from relations (14)- (18), uniform finite-difference approximations are used to describe both the initial state and the deformed state [1,2]. The finite-difference analogs of equilibrium equations (11), discretized with respect to the node point of the main grid i, are written in the following form (20), (21): where (Q 13 ) i−1/2 = 1 where (q 1 ) i , (q 2 ) i are the grid functions of the surface loading components. Since for approximation of VAT parameters in initial and deformed states the same type finite-difference approximations (14)- (19) are used, then at (q 1 ) i = (q 2 ) i = 0 the nondeformed state given by relations (14), (15) is the exact solution of grid equations (20), that confirms the correctness of developed difference approximations (14)-(21). To construct a computational algorithm for solving the system of nonlinear grid equations (20), a quasi-dynamic form of the establishment method is used [3][4][5]9]. To construct the evolutionary problem, the grid analogs of equilibrium equations (20) are replaced by nonstationary equations that have the same form as the shell's equations of motion in a viscous medium [1][2][3]5]. The finite-difference analogs of the equilibrium equations (20) in operator form are represented as (22): where [L ∆s (U k )] i are the corresponding generalised finite-difference operators for the vector U k of the grid displacement functions; (k = 1, 2).
For the case of "tracking" load, to simplify the computational procedure, it is advisable to build an iterative process with respect to the displacements of the local basis u, w, followed by recalculation in the x, y coordinate system (Figure 1). Then the non-stationary setting method equations for the node point of the main grid will be written in the following way (23): where m k = ρh; ρ is the density; The grid functions of displacements x i , y i of the node points in the x, y coordinate system can be recalculated through the grid functions of the displacement components u i , w i in the following way (25), (26) (Figure 1): where Thus, the difference approximation of the non-stationary equations (23) leads to the iterative process (24) of finding the solution to the original stationary problem (22). The computational algorithm for solving the nonlinear problem in the case of conservative loading is constructed similarly to (23), (24). Iteration process parameters-specific viscosities of the medium ε k (i) and step in time ∆t-are determined from the condition of accelerated convergence and stability of the difference scheme [3,4]. The works related to the study of inhomogeneous elastic shells are reviewed in publications . Estimating formulas, taking into account the structure of equations (23), can be written in the form (27) [5,8,9]: where µ 1(k) and µ 2,(k) are the smallest and largest eigenvalues for the corresponding difference operators in equations (22); a ε,(k) and a t,(k) are the correction factors close to unity. The time step ∆t for the whole difference scheme is determined from the condition of the form (28): For nonlinear problems, the exact determination of the boundaries of the difference operator spectra is associated with considerable mathematical difficulties, so µ 1(k) and µ 2(k) are evaluated within the framework of linear relations with appropriate simplifications in the original equations. The evaluation formulas for µ 1(k) and µ 2(k) for the case of a single-layer orthotropic construction can be represented in the following way: • the largest eigenvalues (30): where k 10 , k 20 are the characteristic values of the shell curvature parameters. Based on the developed mathematical models and computational procedures, (1)- (30) there was conducted a study on the influence of physico-mechanical characteristics of isotropic and composite materials on the deformation features of the elastic spherical shell in the subcritical and post-critical regions. A single-layer spherical shell of constant thickness h under the action of a static uniformly distributed external load of the "tracking" type with intensity q w was considered. Geometric parameters of the shell: R 10 /h = 200; x c = 0.707·R 10 ; y c = 0; θ 0 = θ L = 15 0 , where R 10 -the initial radius of curvature ( Figure 5). At the edges of the shell s = s 0 and s = s L the rigid fixing boundary conditions were considered (12). The intensity of the surface-distributed external pressure q w was determined through the value of the load parameter p q in relation to the critical pressure q cr for a closed isotropic spherical shell with radius R = R 10 : The intensity of the surface-distributed external pressure q w was determined through the value of the load parameter p q in relation to the critical pressure q cr for a closed isotropic spherical shell with radius R = R 10 (31): The influence of physical and mechanical characteristics of isotropic and orthotropic materials on peculiarities of nonlinear deformation of the shell in the subcritical and post-critical regions has been investigated. Carbon fiber composite with the following combinations of Young's moduli E 1 and E 2 in formula (6) [3,5,8] was considered as an orthotropic material (32): E 1 > E 2 : E 1 = 11.25·E 2 ; ν 21 = 0.289; E 1 < E 2 : E 1 = 0.083·E 2 ; ν 12 = 0.289.
The value of Young's modulus E of the isotropic material in relation to the moduli E 1 and E 2 of the carbon fiber-reinforced plastic was: E = 0.39·E max , where E max = max(E 1 , E 2 ). Poisson's ratio of an isotropic material: ν = 0.3. The studies were conducted for load parameter values p q = 0.566; 1.18; 2.36. Due to the symmetry of the problem, the calculation domain was 1 2 of a shell structure bounded by the parameters: θ 0 ≤ θ ≤ π/2. In the numerical solution of the problem, the number of sampling points 1 ≤ i ≤ N was assumed to be N = 57. The node point with index i = 1 corresponded to the clamped edge of the shell s 0 , and at the node point i = 57 at θ = π/2 the conditions of solution symmetry were fulfilled. Figure 6 illustrates the change in k 1 /k 10 curvatures at p q = 2.36 for the considered variants of isotropic and orthotropic materials. The horizontal line1 corresponds to the curvature of the spherical shell in the initial (undeformed) state: k 10 = const. Curve 2 describes the change in curvature k 1 /k 10 for an isotropic material (E 1 = E 2 ), and curves 3 (E 1 > E 2 ) and 4 (E 1 < E 2 )-for an orthotropic material. When building the dependencies in Figure 6, Microsoft Office Excel application was used.

Conclusions
The results of the studies conducted have shown that the physical and mechanical characteristics of isotropic and composite materials significantly affect the nature of the deformation of the clamped spherical shell in both the subcritical and post-critical domains. The most optimal variant in terms of increasing the bearing capacity and reducing the deformability of the shell is the carbon fiber reinforced plastic with the Young's modulus ratio E 2 > E 1 . At the same time, maximum displacements u,w in the shell for isotropic material were of order max(u,w) ≈ 5·h at the maximum rotation angles of normal max(∆θ) ≈ 0.15·π, and for orthotropic material maximum values of deformed state parameters were of order max(u,w) ≈ 20·h, max(∆θ) ≈ 0.2·π-for the case E 1 > E 2 , and max(u,w) ≈ 2·h, max(∆θ) ≈ 0.07·π-for E 2 > E 1 . The results obtained can be applied to the use of additive technologies in the aviation and space industries.
The developed mathematical models and computational algorithms can be applied in the future to study shells of rotation made of hyperelastic (non-linearly elastic) materials and soft shells.