CompositeShell

Daniel Hoover and Ashok V. Kumar

https://doi.org/10.1016/j.finel.2022.103794

sUMMARY

An immersed boundary 3D shell element was developed that is based on Mindlin-Reissner shell theory assumptions and uses quadratic B-spline approximation for the solution. It enables mesh independent analysis wherein the surface representing the shell geometry is defined independently and not by the mesh. The shell geometry is immersed in a uniform background mesh whose elements are 3D B-spline shell elements. By making the geometric model independent of the mesh, the typical difficulties associated with mesh generation can be avoided. The quadratic 3D B-spline shell elements in the mesh can represent the displacement field as a tangent continuous piece-wise polynomial approximation. The nodes of the element have three translational and no rotational degrees of freedom. The element is formulated under small-deformation and plane stress assumptions. Constructing the element based on shell theory enables the use of effective properties to model thin fiber reinforced laminated composite shell-like structures. The element formulation and the mesh independent approach are validated against a series of common shell element benchmark problems.

Introduction

In traditional finite element method, the mesh approximates the geometry of the structure and that can lead to difficulties in mesh generation when the geometry is complex. Meshless methods [1] seek to avoid constructing a mesh by using only nodes distributed in the domain. But the lack of connectivity between the nodes makes these methods computationally inefficient. In the mesh independent approach, the geometry is represented using equations or tessellated approximations and is immersed in a background mesh. This approach has been referred to by several names such as immersed boundary finite element method (IBFEM) [2], embedded domain method [3], fictitious domain method [4] and finite cell method (FCM) [5]. Hollig [6] developed a mesh independent method using weighted extended B-splines or web-splines. All these methods utilize a background mesh with the geometry represented independently, often using equations, and embedded in this mesh. But they differ in how boundary conditions are applied and how area and volume integrations are computed. The extended finite element method (X-FEM) [7] is also an example of a mesh independent method, applied mainly to fracture mechanics, where the geometry of a crack in the structure is represented as an equation. The main motivation of more recent efforts to developed mesh independent approach is to fully automate the mesh generation process for all geometry. As the background mesh does not have to conform to the geometry, it is easy to generate such a mesh regardless of the complexity of the geometry. If accurate equations of the geometry are utilized, then the inaccuracy associated with approximating the geometry with a mesh can be avoided. The background mesh is used only for approximating the trial and test functions as piece-wise polynomials. Therefore, it can be a uniform structured mesh where all elements are rectangular (2D) or cuboids (3D). As none of the elements are distorted into quadrilaterals (2D) or hexahedrons (3D), the derivatives of the field variables with respect to the global Cartesian coordinates are also piece-wise polynomials and therefore errors associated with quadrature can also be avoided. Many traditional shell elements are based on Mindlin-Reissner shell theory, because formulation that use Kirchhoff Love theory require C¹ continuous interpolation or approximation of the displacement field and because the latter theory is only applicable to very thin shells where the transverse shear strain energy is negligible. Mindlin-Reissner shell theory includes transverse shear strain energy in the formulation and assumes that a through-thickness line normal to the mid-surface remains straight during bending but need not remain normal to the mid surface [8]. The main challenge for Mindlin shell elements is shear locking that occurs due to parasitic shear strain that arises due to the low order interpolations and the kinematic assumptions of the shell theory. Several techniques for overcoming these problems have been developed including selective reduced order integration, assumed natural strains (ANS) or mixed interpolation of tensorial components (MITC) [9] which has led to several triangular [10] and quadrilateral shell elements [11]. More recently, non-traditional approaches including scaled boundary finite element method (SBFEM) [12] and isogeometric method (IGM) [13] have been used to develop effective shell elements. In all these approaches, a mesh is needed that approximates the geometry of the shell and in IGM both the geometry and the solution are constructed using the same basis functions. IGM has also used several types of splines including NURBS [14] and T-splines [15] for geometric representation and for the trial solution. These basis function can provide C¹ or higher degree of continuity and have also been used with Kirchhoff shell theory [16,17]. However, IGM differs from the mesh independent approach because it requires that both the geometry and the displacement field should be approximated using the same basis functions. In contrast, in the mesh independent approach, the geometry can be independently defined as a boundary representation (B-Rep) using parametric curves and surfaces or even as a tessellated approximation. Geometry defined in these formats are readily available from CAD software and it can be directly used for the analysis without further modifications or approximations.

In this article, a mesh independent approach (referred to here as immersed boundary finite element method or IBFEM) is presented for shell-like structures. The term ‘immersed boundary method’ was originally used in computational fluid dynamics (CFD) [18] especially for fluid-structure interaction problems where the interface boundaries are immersed in a Cartesian background mesh. To model shells, the geometry of the shell represented by its mid-surface, is immersed in a uniform background mesh consisting of 3D cuboid elements as shown in Fig. 1. The figure also shows the typical triangulation used within an element for quadrature. This approach has been used in past work to develop immersed boundary shell elements with both a continuum formulation [19] where 3D stress-strain relations are used and with Kirchoff-Love assumptions [20]. Mindlin-Reissner shell theory assumptions are used here to formulate the 3D immersed boundary element, but the transverse shear strains are not assumed to be constant through the thickness. The element uses B-spline approximations for the trial and test functions. By formulating the element using a shell theory, it becomes possible to use effective properties of laminate composites. Much work has been done on modeling laminated composites using layer wise models [21,22], zig-zag models [23], unified framework [24] and variable kinematic models [25]. Other shell element formulations that use B-splines such as the isogeometric approach have also been used for shell formulations based on these theories [26,27]. The feasibility of using shell theory-based lamination models for estimating composite shell behavior using immersed boundary approach is studied here using a simple model wherein the transverse shear strains are assumed to vary through the thickness to the extend allowed by the element shape functions. The feasibility of using more accurate equivalent single layer models and layer-wise models will be studied in the future.

In the immersed boundary approach, as the mesh does not conform to the geometry, the boundaries on which essential boundary conditions must be applied may not have nodes on them. Many methods for imposing essential boundary conditions for mesh independent analysis have been proposed including the penalty method [28], stabilized Lagrange multiplier method [4], and Nitsche’s method [29,30]. In this article, the step boundary method (SBM) [2] is used which is particularly suited for a mesh independent approach because it can impose boundary conditions using the equations of the boundary curves/surfaces. Boundary conditions can be applied even if the boundaries do not have any nodes on them and when the shape functions used for the analysis do not satisfy the Kronecker-delta property. SBM has been shown to be effective for plate elements [31] and for a variety of other applications using both Lagrange and B-spline elements [[32], [33], [34], [35]]. The step boundary method is applied here to handle both isotropic and laminated composite shell models using the immersed boundary shell elements. $E_{L^{2}} = \frac{{‖ u - u_{h} ‖}_{L^{2}}}{{‖ u ‖}_{L^{2}}}, E_{s} = \frac{{‖ u - u_{h} ‖}_{s}^{2}}{{‖ u ‖}_{s}^{2}} .$

Examples

6.1. Square plate

A square plate of dimensions 10×10, thickness t =0.1 inch and subjected to uniform transverse load was analyzed using 3D B-spline shell elements. Fig. 2(a) shows the mesh with a square surface representing the mid-surface of the shell embedded in the mesh. The element size was selected such that it is thinner in the direction normal to the plate. The deflection of the plate was computed with both fully clamped and simply supported boundary conditions on all sides. The magnitude of the transverse pressure load is P = 100 and the material properties are E = 1.12E¹⁰ and Poisson’s ratio is 0.27.

Table 1 shows the maximum displacement at the center of the plate in the z-direction computed using the 3D B-spline shell elements, traditional shell elements and the analytical/series solution [40]. For the comparison with traditional shell, ANSYS SHELL181 element was used with a 100×100 element mesh to obtain a converged solution. The ANSYS element is a four node quadrilateral shell element with six DOF per node that is designed to model thin to moderately thick shells and composite structures [41,42].

Table 1. Square plate: maximum deflection. $w_{\max}$

Element type	Mesh	SS $w_{\max}$ $v_{\max}$	Clamped $w_{\max}$ $v_{\max}$	SS von Mises stress	Clamped von Mises stress
Quadratic B-spline element	20 × 20	4.048E⁻³	1.243 × 10⁻³	3.416 × 10⁵	2.439 × 10⁵
Quadratic B-spline element	40 × 40	4.062E⁻³	1.260 × 10⁻³	3.683 × 10⁵	2.692 × 10⁵
Cubic B-spline	10 × 10	4.060E⁻³	1.261 × 10⁻³	3.310 × 10⁵	2.715 × 10⁵
element	20 × 20	4.075E⁻³	1.262 × 10⁻³	3.547 × 10⁵	2.753 × 10⁵
ANSYS	100 × 100	4.064E⁻³	1.259 × 10⁻³	3.580 × 10⁵	2.643 × 10⁵
Series solution		4.035E⁻³	1.256 × 10⁻³	–	–

Fig. 3 shows a plot of the normalized z-direction displacement versus increasing mesh density obtained using the quadratic element for (a) the simply supported (SS) plate and (b) the clamped plate. The displacements are normalized with respect to the solutions obtained using ANSYS. The rate of convergence of the relative error norms for the simply supported and clamped plates obtained using cubic B-spline elements are shown in Fig. 4(a) and (b) respectively for different thickness. The dotted line in the figure corresponds to a cubic rate of convergence. These plots indicate that the rate of convergence does not deteriorate with decreasing thickness which implies that these elements do not suffer from shear locking even though they are based on Mindlin shell theory. The quadratic B-spline elements do show that the rate of convergence reduces with thickness. However, shown in Table I, they do yield accurate answers for thickness of 0.1 for this example.

6.2. Scordelis-Lo roof

Scordelis-Lo roof problem is a thin curved 3D shell-like structure that is commonly used as a benchmark problem to test the accuracy for shell elements. It has a singly curved geometry as shown in Fig. 5(a) with radius R = 25, length L = 50, a radial angle of α = 80° and thickness t = 0.25. A transverse load of P = 90 is applied downward on the top surface of the structure. The ends of structure on either side of the length are attached to ridged diaphragms and the other two side are free of constraints. The rigid diaphragms prevent radial motion while allowing motion in the axial direction. Due to the symmetry of the geometry, only one-fourth of the structure is modeled. The model consists of the mid-surface of the shell representing its geometry immersed in a 3D background mesh as illustrated in Fig. 9(b). The material used for this problem are: Young’s modulus, E = 4.32E8 and the Poisson’s ratio ν = 0.

Due to symmetry, the maximum displacement occurs at the corner of the symmetric and free edges as shown in Fig. 6. This result obtained using a model containing 2772 nodes compares favorably with the analytical solution of vmin = −0.3086 [43].

Table 2 shows results obtained for the maximum displacement and von Mises stress for quadratic and cubic B-splines for comparison with results from ANSYS for roughly the same number of nodal degrees of freedom. All three results are close and agree reasonably with the analytical solution.

Table 2. Scordelis-Lo Roof results.

Element type	Max Displacement	Max Von Mises stress
Quadratic B-spline	0.3068	3.612 × 10⁵
Cubic B-spline	0.3075	3.740 × 10⁵
ANSYS	0.3067	3.446 × 10⁵

In Fig. 7, the convergence of displacement and relative error norm are shown. The displacement is normalized with respect to the analytical solution. The plot of relative error in s-norm and L² norm versus the element size in a log-log plot shows an approximately quadratic rate of convergence.

6.3. Pinched hemisphere

The pinched hemisphere shown in Fig. 8 is a benchmark problem for shell elements that tests the element’s performance on a double-curved geometry with concentrated applied loads. It is used to test an element’s ability to predict rigid-body rotations in the normal direction to the surface. The geometry is a spherical shell of radius R = 10 and a thickness of t =0.04. Concentrated loads are applied to the surface of the hemisphere with one applied in the positive x-direction at the intersection with the X-axis and the second applied in the negative y-direction at the intersection with the Y-axis, as shown in Fig. 8. Both forces have a magnitude of F = 2. The hemisphere is symmetric such that only one-fourth of the geometry needs to be modeled. To prevent rigid body movement, a point at the intersection of the Z-axis and the surface is fixed. The material properties for this problem are defined as E = 6.825E7 and ν = 0.3. Fig. 12(b) also shows the mesh for the geometry.

The analytical solution available from literature [43] for the pinched hemisphere problem is the maximum displacement which occurs at the point of application of the load and is umax = 0.1850.The results obtained using the 3D quadratic B-spline elements is shown in Fig. 9. Convergence of displacement and relative error norms are plotted in Fig. 10. A quadratic rate of convergence is observed for quadratic B-spline elements by plotting the L² and s-norms relative errors on a log-log plot as shown in Fig. 10. Table 3 shows the results obtained using quadratic and cubic B-spline elements along with results from ANSYS using models that have approximately similar number of degrees of freedom.

Table 3. Pinched hemisphere results.

Element type	Max Displacement	Max Von Mises stress
Quadratic B-spline	0.1802	1.895 × 10⁴
Cubic B-spline	0.1802	1.968 × 10⁴
ANSYS	0.1865	2.302 × 10⁴

6.4. Pinched cylinder

The pinched cylinder problem is another benchmark problem with singly curved geometry, this time with a concentrated load applied. Its deformation is characterized by both bending and membrane forms of deformations. The cylinder is closed at both ends by rigid diaphragms and two opposing concentrated forces are applied inward at the middle of the length on opposite sides of the structure as shown in Fig. 11. The rigid diaphragm prevents displacements in the y-z plane but allows displacement in the x-direction. The length of the cylinder is L = 600, the radius is R = 300 and the thickness is t = 3.0. The applied forces are both F = 1. The material is defined by the Young’s modulus E = 3.0E6 and the Poisson’s ratio υ = 0.3. Due to the symmetry of the geometry along the X-axis, only one eighth of the geometry is modeled geometry is modeled.

Fig. 15 Fig. 11. Pinched cylinder problem (a) Geometry (b) Mesh example.

Fig. 12 shows the results obtained using a background mesh that contains 5856 nodes. Further refinement of the mesh yields slightly higher values and is consistent with ANSYS’s SHELL181 elements which converges towards approximately 1.833×10E−5. The convergence of the max displacement normalized with respect to this value is shown in Fig. 13 (a) and the convergence rate of the relative error norms are compared with the quadratic rate of convergence in Fig. 13 (b). Table 4 summarizes the results obtained using quadratic, cubic B-spline and ANSYS models with similar number of degrees of freedom.

Table 4. Pinched Cylinder results.

Element type	Max Displacement	Max Von Mises stress
Quadratic B-spline	1.825 × 10⁻⁵	0.1881
Cubic B-spline	1.846 × 10⁻⁵	0.1857
ANSYS	1.832 × 10⁻⁵	0.1854

6.5. Composite square plate

A composite square plate subjected to transverse load is used here as a test problem to validate the use of 3D B-spline elements for laminated composites. A square plate of size 10×10 and thickness of 0.01 is subjected to a transverse load of total magnitude q = 10. The plate is constrained in three different ways for two different fiber orientation of the plies of the laminate material. The first boundary condition is fully clamped and the other two are simply supported. Fig. 14 depicts the two simply supported boundary conditions [22].

$Fig. 14 . Composite square plate problem simply supported boundary conditions (a) SS1 (b) SS2.To test the element, two models were created using different stacking sequences: (i) a two-layer cross-ply [0/90] and (ii) a two-layer angled-ply [45/45]. The two-layer cross-ply is non-symmetric and therefore, it is largely characterized by bending and extensional strains, whereas the two-layer angled-ply is antisymmetric and therefore exhibits membrane shear and torsion in addition to bending and extensional strains. Table 5 shows the deflection at the center of the plate computed using IBFEM and compared with the traditional FEM and analytical solutions. The traditional FEM solutions were obtained using a commercial software [44]. Analytical solution is available only for two of the cases and they were obtained using the classical laminate plate theory (CLPT) [22] that is based on Kirchhoff’s plate theory assumption rather than the Mindlin laminate plate theory. The square plate problem was analyzed with clamped boundary conditions as well as both the simply supported boundary condition applied to the two established two-ply laminates.Table 5 . Composite plate maximum displacement.Laminate BCs Analytical Reddy [22] FEM IBFEM {[0 / 90]}_{T} SS1 −1.6955 −1.6839 −1.703 −1.707 {[0 / 90]}_{T} SS2 – −0.8957 −0.9114 −0.9142 {[0 / 90]}_{T} Clamped – −0.3814 −0.3972 −0.3999 {[- 45 / 45]}_{T} SS1 – −0.6773 −0.7033 −0.7072 {[- 45 / 45]}_{T} SS2 −1.0280 −1.0280 −1.055 −1.060 {[- 45 / 45]}_{T} Clamped – −0.3891 −0.4100 −0.4129In Table 5, the column titled “Laminate” describes the stacking sequence and the column titled “BCs” describes the boundary conditions for the problem. Convergence plots for minimum z-displacement and error s-norms were created for each of these problems and are shown in Fig. 15 .Fig. 15 . Composite square plate problem convergence plots (a) Normalized vmin (b) Relative error norms.6.6. Composite Scordelis-Lo roofThe Scordelis-Lo roof problem is once again used here with laminated composite material to test the composite element formulation. The geometry and boundary conditions are the same as that of the isotropic problem shown in Fig. 5 (a). The upward direction was changed to the Z-direction for this composite problem to match the results from Reddy [22]. The dimensions used are radius R = 300 inches, length L = 600 inches, radial angle 80 degrees, and a thickness of t = 6 inches. The transverse load applied has a magnitude of P = 0.625. The material properties of the lamina are assumed to be the same as for the composite square plate. Since the deformation may not be symmetric depending on the stacking sequence, the whole geometry was modeled. A typical mesh and the computed deformation of the shell are shown in Fig. 15 .Fig. 15. Scordelis-Lo roof problem (a) Mesh example (b) Z-direction displacement results.The shell deformation was computed for both the cross-ply and the angle-ply laminates as described before for the composite square plate problem. In addition, ten-layer versions of a cross-ply, or [0/90]10T, and angle-ply, [45/45]10T was also analyzed. The results are compared to FEM solutions found by Reddy [22] and identical problems performed in a commercial FEM software [44] as shown in Table 6 . Unlike for the standard FEM, the results shown for IBFEM are computed using the assumption that transverse shear is not constant, and it varies as allowed by the displacement field assumed within the element.Table 6 . Composite Scordelis-Lo roof maximum z-displacement. Laminate Reddy [22] FEM IBFEM {[0 / 90]}_{T} −0.5082 −0.509 −0.5271 {[- 45 / 45]}_{T} −0.6760 −0.672 −0.6686 {[0 / 90]}_{10T} −0.2940 −0.2943 −0.2963 {[- 45 / 45]}_{10T} −0.4096 −0.4029 −0.4031The convergence of the maximum displacement for each of the laminates as well as the convergence of the relative error s-norms are shown in Fig. 21 (a) and (b) respectively. Fig. 21 . Composite Scordelis-Lo roof problem convergence plots (a) Normalized w_{\min} (b) Relative error norms.University of Florida Facebook X (formerly Twitter) YouTube LinkedIn Instagram Flickr UF Resources UF News UF Calendar myUFL One.UF Directory Campus Weather Campus Map Student Tours Academic Calendar Events Website UF Website Listing Accessibility Text-Only Version Privacy Policy Regulations{"prefetch":[{"source":"document","where":{"and":[{"href_matches":"\/*"},{"not":{"href_matches":["\/wp-*.php","\/wp-admin\/*","\/ashok-kumar\/wp-content\/uploads\/sites\/398\/*","\/ashok-kumar\/wp-content\/*","\/wp-content\/plugins\/*","\/wp-content\/themes\/hwcoe-ufl\/*","\/*\\?(.+)"]}},{"not":{"selector_matches":"a[rel~=\"nofollow\"]"}},{"not":{"selector_matches":".no-prefetch, .no-prefetch a"}}]},"eagerness":"conservative"}]}/* <![CDATA[ */
var hwcoe_ufl_sitedata = {"home_url":"https:\/\/akumar.mae.ufl.edu\/","theme_url":"https:\/\/akumar.mae.ufl.edu\/wp-content\/themes\/hwcoe-ufl","max_main_menu_items":"7","mega_menu":"1"};
/* ]]> */$

Immersed boundary Mindlin-Reissner 3D shell element for modeling isotropic and laminated composite shells