A Nonlinear Shell Finite Element Modeling Method with Vector Rotation Spherical Coordinate Parameterization

Abstract

The present invention discloses a non - linear shell finite - element modeling method with vector - rotation spherical - coordinate parameterization, which relates to the technical fields of non - linear structural dynamics and non - linear finite - element modeling. The method includes the following steps: extracting a two - dimensional mid - surface shell from a thin - walled three - dimensional structure; discretizing the two - dimensional mid - surface shell into regular quadrilateral shell elements, and the quadrilateral shell elements adopt nine - node second - order accurate quadrilateral shell elements; after the element discretization, initializing the values of five parameters in the shape function of the nine - node second - order shell element in the Lagrangian coordinate system and the vector - rotation spherical - coordinate system; establishing a structural motion control equation for the shell element. The present application provides a non - linear shell finite - element modeling method with vector - rotation spherical - coordinate parameterization. The technical solution of the present application will define the mid - surface normal vector and directly derive the large deformation and motion control equations of the shell, and establish a new type of geometrically exact shell element.

Description

Technical Field

[0001] This application relates to the technical fields of non - linear structural dynamics and non - linear finite element modeling, and specifically relates to a parametric finite element shell element based on a vector rotation spherical coordinate system, which can be applied to the large - deformation non - linear modeling of aircraft thin - wall structures. Background Art

[0002] In order to reduce structural weight, thin - wall structures are widely used in civil aircraft design. When using finite element technology to perform stress analysis and strength check on civil aircraft structures, shell elements are very suitable for establishing finite element models of thin - wall structures. Because shell elements have strong flexibility when dealing with complex geometric shapes and boundary conditions. In the past few decades, shell theory research has received great attention. According to the different basic theories adopted, the existing shell element modeling technical solutions can be divided into two categories. One is the dimensionality reduction modeling of solid elements, and the other is direct modeling.

[0003] The first technical solution makes full use of the well - developed three - dimensional solid modeling theory. By reducing the order of the three - dimensional basic equations in continuum mechanics, the two - dimensional shell element control equations are derived, and then the corresponding shell element numerical model is obtained by directly applying the conventional finite element interpolation technology. This theory was first proposed by Ahmad et al. and has since developed rapidly. Representative works include those of Hughes et al., and Bathe and Dvorkin. The linear shell element QUAD4 in the commercial program NASTRAN also belongs to this type of technical solution. In addition, Bathe, Hughes, and Crisfield have also conducted extensive and in - depth discussions on the dimensionality reduction modeling theory of solid elements for shells. The basic theory behind this technical solution is relatively simple and is very suitable for solving linear problems, but the computational efficiency is low. In this solution, based on the parametric theory of finite rotation, two degrees of freedom are used to describe the rotation of the shell element normal vector. For example, Hughes and Liu use two approximate Rodrigues parameters to describe the rotation. Due to the introduction of the small - angle rotation assumption, this shell element cannot handle large rotations.

[0004] Technical Solution 2 is different from Solution 1. Instead of using a reduced three-dimensional equation, it defines a normal vector perpendicular to the neutral plane of the shell element and directly derives the governing equations to describe the deformation and motion of thin-walled structures such as shells, forming a brand-new shell theory. Based on this theory, a more advanced finite element model can be constructed. This type of technical solution was first proposed by Argyris and his collaborators and is called the direct method. Eriksen and Truesdell adopted this method and for the first time simplified the shell into a geometric surface containing two tangential direction vectors, developing the corresponding shell theory. Simo et al. carried out similar work and developed a geometrically exact shell model based on stress-resultants. The concept of a geometrically exact shell model means that when deriving the basic equations of the continuous mechanics of the shell, except for the initially introduced kinematic assumptions, no other assumptions are introduced during the derivation process, and the introduced kinematic assumptions are strictly adhered to during the subsequent derivation process. Therefore, the shell elements established in Technical Solution 2 have very high accuracy and are suitable for dealing with large deformations and large rotations. Simo et al. accurately described the finite rotation using the rotational Rodrigues expression. In addition, Crisfiled and Moita [4] gave the treatment of finite rotation in the co-rotational expression. Different from linear displacements, rotations need to be described by rotation tensors, so the accurate parametric description of rotation vectors is quite complex, consuming a large amount of computer resources and reducing the computational efficiency. Currently, it is still an unsolved problem, seriously restricting the application of Technical Solution 2.

[0005] The existing Technical Solution 1 reduces the three-dimensional basic equations in the continuous medium mechanics to obtain the governing equations of two-dimensional shell elements, and then directly applies the finite element interpolation technique to obtain the corresponding numerical model of the shell elements. The modeling theory is relatively simple, but the computational efficiency is low and the computational accuracy is not high. Usually, two rotational degrees of freedom are used to describe the small-angle rotation of the normal vector of the shell element, which is only suitable for solving linear small-rotation problems and cannot handle large rotations. This technical solution cannot solve the large-deformation problem of thin-walled structures during loading. The structures designed according to this technical solution may be too conservative or may cause catastrophic risks.

[0006] Technical Solution 2 directly derives the governing equations to describe the deformation and motion of thin-walled structures such as shells, which is a brand-new shell theory. However, the improvement in modeling accuracy makes the theoretical derivation quite complex and the computational efficiency extremely low. The parametric description of the rotation of the shell normal vector further complicates the problem. Rotations need to be accurately described by rotation tensors. Even after parameterization, they are still different from displacement vectors and cannot be directly interpolated using finite elements, requiring further special treatment. The low efficiency of Technical Solution 2 limits its application. Summary of the Invention

[0007] To address the above problems, the present application provides a non-linear shell finite element modeling method with vector rotation spherical coordinate parameterization. The technical solution of the present application will define the normal vector of the mid-surface and directly derive the large deformation and motion control equations of the shell, and establish a new geometrically exact shell element. Different from the conventional treatment scheme, instead of parameterizing the rotation of the normal vector, spherical coordinates are introduced to describe this vector, which greatly simplifies the shell theory and improves the calculation efficiency. The spherical coordinate description of the normal vector does not introduce the small angle rotation assumption, and the two spherical coordinate parameters can accurately define the normal vector after arbitrary rotation. The arbitrary large rotation of the shell structure is described by the relative amount obtained by subtracting the initial value from the spherical coordinate parameters of the normal vector in the deformed configuration, solving the problem that the first technical solution cannot handle large rotations. The combination of the three translational displacement vector components and the two spherical coordinate components of the normal vector accurately describes the displacement and rotation of the shell element. According to the definition of the Green-Lagrange strain tensor, a new expression between the new strain and the translational displacement and rotation spherical coordinates is strictly derived, retaining the second-order term of the deflection strain. This new strain-displacement and spherical coordinate expression is a strictly geometrically non-linear description and can handle large deformation problems. In addition, the introduction of the spherical coordinates of the normal vector eliminates the need for further rotation parameterization, and it can strictly describe the rotation of the shell structure itself, saving the calculation cost and improving the calculation efficiency. Therefore, the present technical solution solves the problems of complex rotation parameters and low calculation efficiency in the second technical solution.

[0008] The present invention is implemented through the following technical solutions:

[0009] According to the first aspect of the technical solution of the present invention, a non-linear shell finite element modeling method with vector rotation spherical coordinate parameterization is provided, and the method includes the following steps:

[0010] (1) Extract a two-dimensional mid-surface shell from the thin-walled three-dimensional structure;

[0011] (2) Discretize the two-dimensional mid-surface shell into regular quadrilateral shell elements, and the quadrilateral shell elements adopt nine-node second-order accuracy quadrilateral shell elements;

[0012] (3) After the unit discretization of the five parameters in the shape function of the nine-node second-order shell element, initialize the values in the Lagrangian coordinate system and the vector rotation spherical coordinate system;

[0013] (4) Establish a structural operation control equation for the shell element.

[0014] Further, in step (2), the shape function of the nine-node second-order accuracy quadrilateral shell element is expressed as

[0015]

[0016] where H i is the unit interpolation function, ui For the coordinate degrees of freedom u of each node i = [u1 u2 u3 θ φ] T , where u1, u2, and u3 are the three components of the linear displacement, and φ and θ are the two components of the rotational displacement described in spherical coordinates.

[0017] Furthermore, step (3) specifically includes:

[0018] (31) The interpolation function H of the nine-node second-order shell element i Adopts the Lagrangian quadratic element to express the nodal displacement motion, and describes the linear displacement of each node in the Lagrangian coordinate system;

[0019] (32) The interpolation function H of the nine-node second-order shell element i Adopts the vector rotation spherical coordinate system to express the large rotation deformation of the nodal angular displacement, and introduces the spherical coordinate definition of the normal vector.

[0020] Furthermore, in step (31), the thickness h of the shell and the mid-surface area Ω are given in the initial reference coordinate system B. In the reference configuration, the position coordinates x0 of any material point on the mid-surface of the shell and the material coordinate ζ along the mid-surface normal vector are given, then the position coordinates of any material point in the shell are

[0021]

[0022] where α1 and α2 are the material point coordinates of the mid-surface of the shell, and the coordinates α1, α2, and ζ are called curvilinear coordinates. Then the position vector of the material point in the deformed configuration is

[0023]

[0024] In the formula, x0 is the position vector of the material point on the mid-surface of the shell, u(α1,α2) is the mid-surface displacement vector of the shell, is the mid-surface normal vector after the shell is deformed.

[0025] Furthermore, in step (32), the normal vector is defined as

[0026]

[0027] As Figure 3 shown, after the normal vector is orthogonally decomposed in the spherical coordinate system, it can be described by the direction sine function and direction cosine function of the spherical coordinate variables θ and φ.

[0028] Furthermore, in step (4), the motion control equation of the shell element is

[0029]

[0030]

[0031] Among them, the symbol represents the first-order derivative with respect to time, (·) ,1 and (·) ,2 represent the spatial derivatives with respect to the coordinates α1 and α2, (·) T denotes matrix transpose; and represent the inertial force contribution terms, N 1, N 2 and N 3 are the elastic forces caused by deformation, M 1 and M 2 are the elastic moments caused by deformation; the right side of the equation is the external load contribution terms, f is the external force, m is the external moment, and the symbol (·) * denotes the external moment m described in the local coordinate system; E f matrix is the spherical coordinate transformation matrix, E B is the reduced matrix after truncating the column corresponding to the third component of the rotation vector of the E f matrix.

[0032] In the shell theory, along the normal direction of the middle surface, that is, along the thickness direction of the shell, there is no deformation. Therefore, N 3 has no derivative and only appears in the moment balance equation. In addition, for the corresponding moment M 3 = 0, which means that only 5 degrees of freedom are retained after ignoring the in-plane rotation degrees of freedom, including 3 linear displacements and 2 angular displacements. The subscripts 1, 2, and 3 in the above equations represent the three principal axis directions of the local coordinate system. The subscripts 1 and 2 are in the neutral layer, and the subscript 3 represents the normal direction.

[0033] Ignoring the inertial force in the above equation, the static deformation control equation is obtained.

[0034] Furthermore, in step (4), the following strain-displacement relationship is given

[0035]

[0036]

[0037]

[0038]

[0039]

[0040] The above strain-displacement relation is applicable to the case of arbitrary displacements, large rotations, and large strains of elastic bodies with thin-walled structures. Here, ε ij , where i, j = 1, 2, 3 are the non-zero components of the Green-Lagrange strain tensor, and are dimensionless orthogonal tangential vectors on the mid-surface of the plate and shell, and R1 and R2 are the radii of curvature of the thin-shell element. Additionally, ζ 2 corresponds to the second-order term of the deflection strain.

[0041] Furthermore, in step (4), the structural operation control equation of the shell element is derived from Hamilton's principle. The specific derivation process is as follows:

[0042] Taking the derivative of the position vector of an arbitrary material point in the time domain, the velocity vector

[0043]

[0044] is obtained. Here, the vector X is the position vector of an arbitrary material along the normal direction of the mid-surface of the thin shell after deformation in the inertial system. The vector u is the displacement of the reference point on the mid-surface of the thin shell. Thus, the kinetic energy K of the shell can be written as

[0045]

[0046] where ρ is the material density. The variation of the kinetic energy can be written as

[0047]

[0048] In the formula, h and represent the linear momentum vector and the angular momentum vector respectively, and their specific expressions can be obtained from the variation process of the kinetic energy.

[0049] The variation of the strain energy of the shell is

[0050] (9)

[0051] where τ *ij (i, j = 1, 2, 3) is the Cauchy stress tensor described in the convected coordinate system; substituting the strain components into the expression of the variation of the strain energy and integrating along the thickness, the equivalent form of the variation of the strain energy can be obtained

[0052]

[0053] where, N 1, N 2, and N 3 are the elastic forces described in the coordinate system of the mid-surface of the thin shell, M 1, and M2 is the elastic moment described in the middle surface coordinate system of the thin shell. The virtual work of external forces is expressed as

[0054] δW = ∫ Ω (δ u T f + δψ T m )dΩ (11)

[0055] where δψ represents the virtual angular displacement, f and m represent the forces and moments acting on the unit area of the middle surface of the shell. Substituting the variation of kinetic energy, the variation of strain energy, and the virtual work of external forces into Hamilton's principle, the motion control equations of the shell are finally obtained

[0056]

[0057]

[0058] Ignoring the inertial forces in the above equations, the static deformation control equations can be obtained.

[0059] Furthermore, the motion control equations of the shell structure in step (4) are divided into two cases. One is the ordinary differential equation of shell dynamics, and the other is the nonlinear geometric equation of shell structure statics. Both are solved by the linearized equation iteration method, predicting the deformation increment each time. When the deformation does not meet the convergence requirement, it is necessary to update the stiffness, deformation, and elastic load, repeating the above steps with a maximum of 10 iterations; when the deformation iteration result meets the convergence requirement, that is, the linear iteration converges, the response of any large deformation of the structure is the required result, and the above solution process ends.

[0060] According to the second aspect of the technical solution of the present invention, there is provided a use of the nonlinear shell finite element modeling method according to any one of the above aspects in modeling to construct a parametric finite element shell element based on a vector rotation spherical coordinate system.

[0061] According to the third aspect of the technical solution of the present invention, there is provided a parametric finite element shell element based on a vector rotation spherical coordinate system, characterized in that the parametric finite element shell element based on a vector rotation spherical coordinate system is modeled and constructed by using the nonlinear shell finite element modeling method according to any one of the above aspects.

[0062] The beneficial effects of the present invention are as follows:

[0063] 1. Except for introducing the initial small deflection assumption, no other assumptions are introduced in this element, so large deformation problems can be directly handled.

[0064] 2. Instead of parameterizing the rotation of the normal vector, the vector is defined in the spherical coordinate system, and two spherical coordinate parameters are used to describe any rotation of the normal vector.

[0065] 3. Compared with the tensor description of rotation, the spherical coordinate parameterization of rotation is extremely simplified and efficient.

[0066] 4. The technical solution proposed in this application takes both efficiency and accuracy into account. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 Shows the flowchart of the technical solution proposed in this application;

[0068] Figure 2 Shows a 9-node 2nd-order quadrilateral shell element;

[0069] Figure 3 Shows the spherical coordinate description of the normal vector;

[0070] Figure 4 Shows a twisted flat plate configuration;

[0071] Figure 5 Shows the deformation of the twisted flat plate;

[0072] Figure 6 Shows the time history of the deformation of the free end of the twisted flat plate;

[0073] Figure 7 Shows a semi-cylindrical thin-walled shell and the finite element mesh;

[0074] Figure 8 Shows the time history of the displacement of the loading point of the semi-cylindrical thin-walled shell;

[0075] Figure 9 Shows the time history of the centroid curvature strain of the semi-cylindrical thin-walled shell.

[0076] The following further describes the embodiments of the present invention in conjunction with the accompanying drawings. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0077] In addition, the embodiments of the present invention described below in conjunction with the accompanying drawings are exemplary and are only used to explain the embodiments of the present invention, rather than limiting the present invention.

[0078] The technical solution of this application will introduce in detail the theoretical details of the new shell element and the large deformation and nonlinear vibration analysis of thin-walled structures completed by using this technical solution. Among them, the modeling method of this application provides a simple and efficient finite element shell element for solving large geometric deformations while ensuring accuracy. In the spherical coordinate system, three displacement components and two spherical coordinate parameters are used to describe the large deformation of the shell element, solving the problem that the traditional shell element cannot handle large rotations after reduction in Technical Solution 1, and solving the problem of low calculation efficiency caused by large-scale parametric description of the rotational deformation of the element in Technical Solution 2. This element is a new type of high-precision and high-efficiency shell element, and its modeling technical process is as Figure 1 shown.

[0079] The present invention is realized through the following technical solutions: a nonlinear shell finite element modeling method with vector rotation spherical coordinate parameterization, and the method includes the following steps:

[0080] (1) Extract the two-dimensional mid-surface shell from the thin-walled three-dimensional structure;

[0081] (2) Discretize the mid-surface shell into regular quadrilateral shell elements. The number of elements can be determined based on the in-plane size of the thin-walled structure. For improving the analysis accuracy, the conventional nine-node second-order accuracy quadrilateral shell element is adopted, such as Figure 2 shown. The element shape function can be expressed as

[0082]

[0083] where H i is the element interpolation function, and u i is the coordinate degree of freedom of each node

[0084] (3) After the five parameters in the nine-node second-order shell element shape function are discretized, initialize their values in the corresponding coordinate system:

[0085] (31) The nine-node second-order shell element interpolation function H i uses the Lagrangian quadratic element to express the nodal displacement motion: the linear displacement of each node is described in the Lagrangian coordinate system.

[0086] In the initial reference coordinate system B, the thickness h and the mid-surface area Ω of the shell are given. In the reference configuration, the position coordinate x0 of any material point on the shell mid-surface and the material coordinate ζ along the mid-surface normal vector are given. The position coordinate of any material point in the shell can be written as

[0087]

[0088] where α1 and α2 are the material point coordinates of the shell mid-surface. The coordinates α1, α2 and ζ are called curvilinear coordinates. Correspondingly, the position vector of the material point in the deformed configuration can be written as

[0089]

[0090] In the formula, \(x_0\) is the position vector of the material point on the middle surface of the shell, and \(u(\alpha_1,\alpha_2)\) is the displacement vector of the middle surface of the shell. is the normal vector of the middle surface after the shell is deformed. Equation (2) shows that the deformation and displacement of any material point can be decomposed into the displacement \(u(\alpha_1,\alpha_2)\) of the material point on the middle surface and the rotation of the corresponding normal vector

[0091] (32) The interpolation function \(H\) of the nine-node second-order shell element i It is reconstructed by expressing the large rotation deformation of the nodal angular displacement in vector rotation spherical coordinates.

[0092] The technical solution proposed in this application introduces the spherical coordinate definition of the normal vector and directly solves it as an unknown quantity. It avoids the parameterization of the rotation tensor, greatly simplifies the problem, and can directly apply finite element interpolation at the same time. The normal vector can be defined as

[0093]

[0094] For the corresponding geometric description, see Figure 3 .

[0095] The spherical coordinates \(\varphi\) and \(\theta\) are solved as unknowns through the control equations of the shell. In the reference configuration, the initial values of these two coordinates are \(\varphi_0\) and \(\theta_0\), which define the initial normal vector

[0096] Existing technical solutions all use the parameterized expressions of the rotation tensor to describe the rotation of the normal vector. The present invention solves the problems of the existing technical solution 1 that adopts the linear small rotation assumption, introduces two rotation parameters and cannot handle large rotations; and the existing technical solution 2 that uses the rotation parameterization forms such as the Rodrigues expression to describe finite rotations, which further complicates the rotation problem and cannot directly apply finite element interpolation.

[0097] (4) Complete the definition of the values of the shape functions in each coordinate system, and establish the structural operation control equations for the shell element. The shell element is a geometrically nonlinear structural finite element; based on the above-described spherical coordinate description of rotation for large rotation deformation of the structure, which has strong nonlinear characteristics, the basic equations of the shell element are re-derived, and finally a new non-linear geometrically accurate shell element is developed.

[0098] In the strain definition of the shell element, the Mindlin / Reissner hypothesis is introduced: the material line perpendicular to the middle surface of the shell remains straight after deformation and does not undergo tensile deformation. In this hypothesis, it is not restricted that the material line remains perpendicular to the middle surface after the shell deforms, so shear deformation will occur. This hypothesis is applicable to medium-thickness shells. According to the definition of the Green-Lagrange strain tensor in the theory of elasticity, the technical solution proposed in this application gives the following strain-displacement relationship

[0099]

[0100]

[0101]

[0102]

[0103]

[0104] The above strain-displacement relationship is a completely new expression because it introduces spherical coordinate vectors and is suitable for arbitrary displacements, large rotations, and large strains. Among them, ζ 2 corresponds to the second-order term of the deflection strain. The existing technical solution two considers this term as a high-order small quantity and ignores it. The obtained strain-displacement relationship is only suitable for small strain cases. The technical solution proposed in this application uses the fact that during the numerical integration process from the initial time t i to the current time t f , the incremental change of the strain is very small, and approximates the second-order term of the deflection strain corresponding to ζ 2 by its initial value at time t i and updates it step by step in each integration. Therefore, the new shell element can solve the large deformation problem.

[0105] The motion control equation of the shell element is derived from Hamilton's principle. By differentiating the position vector of any material point in the time domain, Equation (2), the velocity vector can be obtained

[0106]

[0107] Thus, the kinetic energy of the shell is

[0108]

[0109] The variation of the kinetic energy can be written as

[0110]

[0111] In the formula h and represent the linear momentum moment vector and the angular momentum moment vector respectively.

[0112] The variational of the strain energy of the shell is

[0113] δV = ∫ Ω ∫ h (τ * 11δε 11 + τ *22 δε 22 + τ *12 δε 12 + τ *13 δε 13 + τ *23 δε 23 )dζdΩ (9)

[0114] where τ *ij (i, j = 1, 2, 3) is the Cauchy stress tensor described in the convected coordinate system; Substituting the strain components into the expression of the variational of the strain energy and integrating along the thickness, the equivalent form of the variational of the strain energy can be obtained

[0115]

[0116] The virtual work of the external force can be expressed as

[0117] δW = ∫ Ω (δ u T f + δψ T m )dΩ (11)

[0118] where f and m represent the force and moment acting on the unit area of the middle surface of the shell. Substituting the variational of the kinetic energy, the variational of the strain energy and the virtual work of the external force into Hamilton's principle, the motion control equation of the shell is finally obtained

[0119]

[0120]

[0121] Neglecting the inertial force in the above equation, the static deformation control equation can be obtained.

[0122] Based on the above structural element control node deformation equation, following the classical finite element analysis process, including assigning the mass property, the unit normal vector, the material vector and the unit stiffness matrix, the load vector and the boundary constraints, and the damping characteristics of the finite element; the deformation calculation based on the Newton-Raphson iterative algorithm and the post-processing of the calculation results to obtain the structural deformation, strain and stress distributions.

[0123] The motion control equation adopts the linearized equation iteration method. When the deformation does not meet the convergence requirement, after calculating the deformation increment of the current step, the stiffness, deformation, and elastic force vector are updated, and further iteration is performed. The above steps are repeated up to 10 times in total; if the deformation iteration meets the convergence condition, that is, the iteration converges, the large deformation response of the structure for any request is the result. End the above solution process.

[0124] In steps (1) to (3), based on the above-derived formulas and combined with the finite element theory, the technical solution proposed in this application has completed the development of a nine-node quadratic precision quadrilateral shell element. This element not only has the function of static analysis, but also has the functions of modal decomposition and dynamic analysis.

[0125] Embodiment

[0126] According to the specific implementation manner of a non-linear finite element modeling method based on vector rotation spherical coordinate parameterization in this embodiment, a new type of 5-degree-of-freedom shell element with three displacement components and two spherical coordinate parameters can directly perform finite element interpolation, and it is a new type of high-precision and high-efficiency shell element. The technical flow chart of the element modeling is as Figure 1 shown.

[0127] Large deformation bending analysis of a twisted flat plate

[0128] A cantilever flat plate with a rectangular cross-section is twisted 90° along the span direction and a concentrated load is applied at the free end, as Figure 4 shown. The length of the twisted flat plate is L = 0.3048 m, the width b = 0.0279 m, and the thickness is 0.00813 m, which is discretized into a 2×8 element mesh (2 shell elements in the width direction and 8 shell elements in the span direction). The material properties of the rectangular plate are set as: Young's modulus E = 200×10 9 Pascals, and the Poisson's ratio ν = 0.22.

[0129] To verify the accuracy of the developed shell element and its ability to handle large deformations, a dynamic analysis was performed on the element. The predicted results were compared and verified with the calculated results of the geometrically exact beam element [5] model. Under the clamped boundary condition, the deformations of the twisted rectangular plate were predicted using two different elements: the currently developed geometrically exact shell element and the publicly published geometrically exact beam element. At the free end, the twisted flat plate deformed under the action of a concentrated shear load of P = 2.224×10 4 Newtons. The deformed configuration is as Figure 5 shown. The shear load linearly increased from zero to 2.224×10 4 Newtons within 5 seconds and then remained at that value.

[0130] The results show that the developed geometrically exact shell element has the ability to predict large deformations. For the comparison and verification with the calculation results of the geometrically exact beam element model, refer to Figure 6 . The time history of the deformation at the midpoint of the free end is described in the figure. The solid line represents the result of the beam model, and the dashed line represents the result of the shell model. The calculation results of the two are basically consistent.

[0131] Stress calculation of a semi-cylindrical thin-walled shell

[0132] One end of the semi-cylindrical thin-walled shell is fixed and the other end is free. The radius is 1.2 m and the length is 2 m. Its loading conditions are as Figure 7 shown. The magnitudes of the concentrated loads are F1 = F2 = F3 = 1.0×10 2 Newtons. The thin-walled structure is discretized into a 4×8 finite element mesh (4 shell elements in the circumferential direction and 8 shell elements in the arc direction). Material properties: Young's modulus E = 73×10 9 Pascals, Poisson's ratio ν = 0.3, and the wall thickness is 0.006 m.

[0133] The dynamic response of the thin-walled structure under the concentrated load is predicted using the new shell element.

[0134] The time history of the displacement and rotation at the loading point can be referred to Figure 8 , and the curvature strain of the geometric centroid of the thin-walled shell can be referred to Figure 9 . Figure 8 And 9 In, the calculation results of the technical solution proposed in this application are compared with the calculation results of the finite element program DYMORE [6] . Figure 8 And 9 In, the solid line is the calculation result of DYMORE, and the dashed line is the calculation result of the new geometrically nonlinear shell element. The calculation results of the two are consistent, indicating that the new geometrically nonlinear shell element has good accuracy.

[0135] The embodiments of the present invention have been described above in conjunction with the accompanying drawings. However, the present invention is not limited to the above specific embodiments. The above specific embodiments are merely illustrative and not restrictive. Under the inspiration of the present invention, those of ordinary skill in the art can also make many forms without departing from the spirit and scope protected by the claims of the present invention. These all fall within the protection scope of the present invention.

Claims

1. A non - linear shell finite - element modeling method with vector rotation spherical - coordinate parameterization, characterized in that, The method is applied to the large deformation non-linear modeling of the thin-walled structure of an aircraft, and includes the following steps: (1) Extract the two-dimensional mid-surface shell from the three-dimensional thin-walled structure; (2) Discretize the two-dimensional mid-surface shell into regular quadrilateral shell elements, and the quadrilateral shell elements adopt nine-node second-order accurate quadrilateral shell elements; Among them, the shape function of the nine-node second-order accurate quadrilateral shell element is expressed as Among them, H i is the unit interpolation function, and u i is the coordinate degree of freedom of each node, u i = [u1 u2 u3 θ φ] T , where u1, u2, and u3 are the three components of the linear displacement, and φ and θ are the two components of the rotational displacement described in spherical coordinates; (3) After the discretization of the five parameters in the shape function of the nine-node second-order shell element, initialize the values in the Lagrangian coordinate system and the vector rotation spherical coordinate system; Among them, step (3) specifically includes: (31) The interpolation function H of the nine-node second-order shell element i The Lagrangian quadratic element is used to represent the nodal displacement motion, and the linear displacements of each node are described in the Lagrangian coordinate system; (32) The interpolation function H of the nine-node second-order shell element i The large rotation deformation of the nodal angular displacement is expressed by using the vector rotation spherical coordinate system, and the normal vector is defined by introducing the spherical coordinate. (4) Establish the motion control equation of the shell structure for the shell element; Among them, the motion control equation of the shell element is where the symbol represents the first-order derivative with respect to time, (·) ,1 and (·) ,2 represent the spatial derivatives with respect to coordinates α1 and α2, (·) T denotes matrix transpose; and represent the inertial force contribution terms, N 1, N 2 and N 3 are the elastic forces caused by deformation, M 1 and M 2 are the elastic moments caused by deformation; the right side of the equation is the external load contribution term, f is the external force, m is the external moment, and the symbol (·) * denotes the external moment m described in the local coordinate system; E f matrix is the spherical coordinate transformation matrix, and E B is the reduced matrix after truncating the column corresponding to the third component of the rotation vector of the E f matrix.

2. The non - linear shell finite - element modeling method according to claim 1, characterized in that, In step (31), the thickness h and the mid-surface area Ω of the shell are given in the initial reference coordinate system B. In the reference configuration, the position coordinates x0 of any material point on the mid-surface of the shell and the material coordinate ζ along the mid-surface normal vector are given. Then, the position coordinates of any material point in the shell are Among them, α1 and α2 are the coordinates of the material points on the mid-surface of the shell. The coordinates α1, α2 and ζ are called curvilinear coordinates. Then in the deformed configuration The position vector of the material point is where \(x_0\) is the position vector of the material point on the shell middle surface, \(u(\alpha_1,\alpha_2)\) is the displacement vector of the shell middle surface, is the normal vector of the middle surface after the shell deforms.

3. The non - linear shell finite - element modeling method according to claim 2, characterized in that, The normal vector in step (32) is defined as 4. The non - linear shell finite - element modeling method according to claim 3, characterized in that, In step (4), it includes the strain-displacement relationship: The above strain-displacement relation is applicable to the case of arbitrary displacements, large rotations, and large strains in elastic bodies of thin-walled structures, where ε ij , i, j = 1, 2, 3 are the non-zero components of the Green-Lagrange strain tensor, and are dimensionless orthogonal tangential vectors of the mid-surface of the plate and shell, R1 and R2 are the radii of curvature of the thin shell element, and ζ 2 corresponds to the second-order term of the deflection strain.

5. The non - linear shell finite - element modeling method according to claim 4, characterized in that, In step (4), the structural operation control equation of the shell element is derived from Hamilton's principle. The specific derivation process is as follows: Derive the velocity vector by taking the derivative of the position vector of any material point in the time domain, where the vector X is the position vector defined in the inertial system for any material point on the normal of the middle surface of the thin shell after deformation, and the vector u is the displacement of the reference point on the middle surface of the shell, and the kinetic energy K of the shell is written as Among them, ρ is the material density, and the variation of kinetic energy is written as In the formula, h and represent the linear momentum moment vector and the angular momentum moment vector respectively, and the specific expressions are obtained from the variational process of kinetic energy; the variational of the strain energy of the shell is δV = ∫ Ω ∫ h (τ *11 δε 11 + τ *22 δε 22 + τ *12 δε 12 + τ *13 δε 13 + τ *23 δε 23 ) dζ dΩ where τ *ij (i, j = 1, 2, 3) is the Cauchy stress tensor described in the body coordinate system; Substituting the strain components into the strain energy variational expression and integrating along the thickness gives an equivalent form of the strain energy variation Among them, N 1、 N 2 and N 3 are elastic forces described in the middle surface coordinate system of the thin shell, M 1 and M 2 are elastic moments described in the middle surface coordinate system of the thin shell, and the virtual work of the external force is expressed as δW = ∫ a (δ u T f + δψ T m ) dΩ where f and m represent the forces and moments acting on the unit area of the middle surface of the shell. Substituting the variation of kinetic energy, the variation of strain energy, and the virtual work of external forces into Hamilton's principle, the motion control equations of the shell are finally obtained Neglect the inertial force in the above equation, and the static deformation control equation is obtained.

6. The non-linear shell finite element modeling method according to claim 1, wherein, In step (4), the dynamic equilibrium equation of the shell structure is an ordinary differential equation, and the statics of the shell structure is a non-linear geometric equation. Both are solved by the linearized equation iteration method, and the deformation increment is predicted according to the iteration result; When the deformation does not meet the convergence requirement, update the stiffness, deformation and elastic load, and repeat the above steps. The maximum number of total iterations is 10 times; when the deformation iteration result meets the convergence requirement, that is, the linear iteration converges, the arbitrary large deformation response of the structure is the result to be obtained, and the above solution process is ended.

Images