# A Geometrically Exact Isogeometric Blended Shell: Formulation, Benchmarking, and Automotive Application.

INTRODUCTIONWe present a geometrically exact isogeometric blended shell formulation. Our formulation is inspired by the blended shell formulation presented in [1]. However, we avoid geometric approximations leading to far more accurate solutions in the presence of curved geometry. We restrict our blending to include both five and six degree-of-freedo m shear deformable Reissner-Mindlin shell kinematics as well as a three degree-of-freedom Kirchhoff-Love shell kinematic. However, our framework is general, and additional shell kinematics can be included (for example the shear deformable five degree-of-freedom rotation-free formulations in [2, 3]) with minimal effort. The blending simplifies structural modeling with smooth, higher-order shell elements in that kinks, complex shell intersections, and boundary conditions can be handled in a simple manner. Additionally, proper blending can dramatically reduce the number of degrees-of-freedom required to solve a problem. For example, five or six degree-of-freedom Reissner-Mindlin control points can be used to capture kinks and intersections while three degree-of-freedom control points can be used elsewhere since the basis is smooth enough to handle a higher-order Kirchhoff-Love kinematic description. This is critical to employing isogeometric shell elements for industrial sized problems. Eliminating rotational degrees-of-freedom, by increasing the smoothness and order of the basis, also leads to elements that are shear locking free, an important consequence of employing the isogeometric paradigm. We benchmark the blended shell against the most common standard commercial shell finite elements. We make these comparisons on several carefully selected benchmark problems in the linear elastic regime.

Kirchhoff-Love (KL) shell theory assumes zero transverse shear strains and is appropriate for thin shells. A finite element implementation of KL theory requires higher-order [C.sup.1]-continuous basis functions. Traditionally, building finite element meshes for [C.sup.1]-continuous basis functions was prohibitively difficult. This resulted in the widespread adoption of thick-shell or shear deformable theories like the Reissner-Mndlin (RM) shell theory as the basis for shell finite elements. RM theory only require [C.sup.[degrees]] basis funtions. However, RM shells introduce additional rotational degrees-offreedom which complicate the implementation. With the advent of isogeometric analysis (IGA) and the adoption of smooth Computer-Aided Design (CAD) basis functions in engineering analysis the use of KL shell theory in the finite element method has greatly increased [4, 5, 6, 7] and resulted in new possibilities for the modeling of complex shell configurations.

B-Spline and Nurbs Fundamentals

A univariate B-spline basis is defined by a knot vector [XI] = {[[xi].sub.1], [[xi].sub.2], ..., [[xi].sub.n+p+1]}, which consists of a non-decreasing sequence of real numbers, [[xi].sub.i] [less than or equal to] [[xi].sub.i+1], i= I, ..., n+ p+ 1 where p is the degree of the B-spline basis functions and n is the number of basis functions. The Ath B-spline basis function of degree p, denoted by [??], can be recursively defined as

[mathematical expression not reproducible] (1)

[mathematical expression not reproducible] (2)

A B-spline curve can be written as

[mathematical expression not reproducible] (3)

where [P.sub.A] is called a control point. To exactly capture conic sections a rational variant of the B-spline basis function is often used. A Non-Uniform Rational B-spline (NURBS) curve can be written as

[mathematical expression not reproducible] (4)

where

[mathematical expression not reproducible] (5)

is a rational basis function and [w.sub.A] is the weight corresponding to the A[t.sup.h] basis function. For notational clarity we often suppress the degree superscript on the basis functions and use [N.sub.A] and [R.sub.A] interchangeably.

THE PRINCIPLE OF VIRTUAL WORK

In this paper the geometrically exact blended shell formulation is restricted to linear elastic problems. The starting point is the principle of virtual work in the spatial configuration,

[mathematical expression not reproducible] (6)

where [sigma] is the Cauchy stress, [epsilon] is the small strain tensor, t is the traction, and b is the body force.

SHELL GEOMETRY AND KINEMATICS

The initial shell body [[omega].sub.0] is described by the reference surface X and an inextensible director vector field |D| = 1 which defines the shell thickness. X is parameterized by the convective coordinates [s.sub.1] and [s.sub.2] while the thickness of the shell is parameterized by [s.sub.3], which ranges from -1 to 1. The reference configuration X in the global coordinates can be written as

[mathematical expression not reproducible] (7)

where h is the shell thickness and the established convention for Latin and Greek indices (i.e., i= 1, 2, 3 and [alpha] = 1, 2) is adopted. The reference director vector D([s.sub.a]) is chosen to be the unit normal to the reference surface

[mathematical expression not reproducible] (8)

where (.), a denotes [partial derivative](.)[partial derivative]d[s.sub.a]. Similarly, the current configuration x can be written as

[mathematical expression not reproducible] (9)

where x([s.sub.a] ) is the reference surface in the current configuration and d is the current director vector. Note that for the KL shell theory the current director vector d is set to the unit normal vector n with respect to the current reference surface, i.e.,

[mathematical expression not reproducible] (10)

where

[mathematical expression not reproducible] (11)

The RM shell theory takes the shear deformation into account so that d[not equal to]n.

The displacement of the shell is written as

[mathematical expression not reproducible] (12)

[mathematical expression not reproducible] (13)

[mathematical expression not reproducible] (14)

where u is the displacement of the reference surface and [delta]d is the increment of the director.

For a blended shell the increment of the current director, [delta]d, is split into contributions from several different kinematic descriptions. The parts of the domain affected by each kinematic description are determined by placing surface control points into sets. The set containing all of the control points which define the surface is denoted by S. We denote the set of control points corresponding to KL shell theory by KL3, the set corresponding a standard RM 5-parameter shell theory by RM5, and the set corresponding to a standard RM 6-parameter shell theory by RM6. Note that

S = K L3 [union] RM5 [union] RM6 (15)

and

K L3 [intersection] RM5 [intersection] RM6 = [empty set]. (16)

Then the displacement u can be further specialized to

[mathematical expression not reproducible] (17)

Note that control points in KL3 are assumed to be associated with basis functions defined over smooth portions of the domain. The control points in RM5 and RM6 are associated with basis functions which are smooth or non-smooth.

THE GRADIENT EXTRACTION OPERATOR

The evaluation of the current director, Ad, for different shell formulations requires that certain quantities be lifted to control points. The lifting is achieved through a simple application of Bezier projection [8]. We compute a gradient extraction operator G = [[G.sub.BA]] such that

[mathematical expression not reproducible] (18)

Using G the gradient of any spline field, [psi], with respect to the convective coordinates can be computed as

[mathematical expression not reproducible] (19)

[mathematical expression not reproducible] (20)

[mathematical expression not reproducible] (21)

where

[mathematical expression not reproducible] (22)

For blended shells the lifting operation is restricted to element domains.

THE KL SHELL DIRECTOR, [delta][d.sup.KL3]

The KL shell director presented in this paper makes no geometric approximations. This leads to very accurate shell formulations for highly curved geometry regardless of mesh density. We have that

p = [X,.sub.1] x [X,.sub.2] (23)

p,[alpha] = [X,[X,.sub.1[alpha]] x [x[X,.sub.2] + [X,.sub.1], x [x[X,.sub.1],2[alpha]] (24)

[mathematical expression not reproducible] (25)

and

[mathematical expression not reproducible] (26)

[mathematical expression not reproducible] (27)

The linearized increment [delta][d.sup.KL3] can be written as

[mathematical expression not reproducible] (28)

[mathematical expression not reproducible] (29)

[mathematical expression not reproducible] (30)

[mathematical expression not reproducible] (31)

[mathematical expression not reproducible] (32)

[mathematical expression not reproducible] (33)

and its derivative as

[mathematical expression not reproducible] (34)

where the superscript on E denotes the node set which is used in the underlying sum. In this case the node set is KL3.

GEOMETRIC APPROMIMATIONS TO [delta][d.sup.KL3]

Several approximations to (28) suggest themselves. These approximations can be achieved by lifting various geometric quantities to control points. The approximation employed in [1] lifts P, [x.sub.,1], and [x.sub.,2] in (28) to the control points. To accomplish this, for each control point we compute

[mathematical expression not reproducible] (35)

[mathematical expression not reproducible] (36)

[mathematical expression not reproducible] (37)

[mathematical expression not reproducible] (38)

and set

[mathematical expression not reproducible] (39)

We then approximate [delta][d.sup.KL3] as

[mathematical expression not reproducible] (40)

[mathematical expression not reproducible] (41)

[mathematical expression not reproducible] (42)

[mathematical expression not reproducible] (43)

[mathematical expression not reproducible] (44)

[mathematical expression not reproducible] (45)

We demonstrate the impact of these geometric approximations on solution accuracy in Figure 4.

RM 5 DOF SHELL DIRECTOR, [delta][d.sup.RM5]

For a Reissner-Mindlin shell the current director d may not be normal to the reference surface due to transverse shear strains. A standard Reissner-Mindlin shear deformable shell theory with two rotations can be formulated as

[mathematical expression not reproducible] (46)

[mathematical expression not reproducible] (47)

[mathematical expression not reproducible] (48)

where

[mathematical expression not reproducible] (49)

and [[beta].sub.A] = [[[beta].sub.A1] [[beta].sub.A2]].sup.T]. Note that [??] is the fiber basis constructed at node A using standard techniques [9]. Note that smooth geometry but only C[degrees] basis functions are required by this theory.

RM 6 DOF SHELL DIRECTOR, [delta][d.sup.RM6]

A standard Reissner-Mindlin shear deformable shell theory with three rotations can be formulated as

[mathematical expression not reproducible] (50)

[mathematical expression not reproducible] (51)

[mathematical expression not reproducible] (52)

where [[omega].sub.A] = [[[[omega].sub.A1] [[omega].sub.A2] [[omega].sub.A3]].sup.T]. Note that only C[degrees] basis functions are required by this theory. Since three global rotations are used this kinematic description can be used in the presence of kinks and intersections with other structural members.

FINAL BLENDING

Plugging the definition of each kinematic description into (17) yields

[mathematical expression not reproducible] (53)

which can be further simplified to

[mathematical expression not reproducible] (54)

Standard approaches can then be used to discretize the virtual work equations [9] given the blended form of the displacement.

NUMERICAL RESULTS

We present four computational examples to explore the behavior of our shell element: a straight cantilever beam, a curved cantilever beam, the Scordelis-Lo roof, and a complex automotive hood inner. Comparisons with standard commercial shell elements in Abaqus and Nastran are included where appropriate. The automotive hood inner is a very complex geometry modeled as a Boundary Representation (BREP) and T-spline. The T-spline is analyzed directly without geometry cleanup or mesh generation, demonstrating the utility of our approach as the basis for integrated isogeometric design/analysis procedures.

In all cases we leverage RM5 and RM6 nodes for the imposition of boundary conditions. Interior nodes are set to KL3, RM5, or RM6 to demonstrate the behavior of the blending. We denote an example where boundary nodes are set to RM6 and interior nodes are set to KL3 by RM6/KL3. Other blendings follow the same notational pattern. Different quadrature schemes are used to alleviate shear and membrane locking. We denote a Gauss quadrature scheme with p + 1 points in each direction over an element by QP1, a reduced Gauss quadrature scheme with p points in each direction over an element by QPO, and the non-uniform, reduced Gauss quadrature scheme described in [10, 11, 12] by QNU. Note that QNU is essentially a one-point quadrature scheme for all orders away from the shell boundaries so it is the most efficient quadrature scheme of the three.

STRAIGHT CANTILEVER BEAM

The straight cantilever beam has a length L = 6 in, a width b = 0.2 in, and a thickness t = 0.1 in, as shown in Figure 1. A traction [q.sub.z] = 5 lbf/in is applied to the right end. The material has Young's modulus E = 1.0 x [10.sup.7] psi, and a Poisson's ratio v = 0.3. This example is frequently used to test the element capabilities to capture simple deformation modes, including linearly varying strains and curvatures. Theoretically, an analytical solution can be usually obtained by simplifying a cantilever beam or plate under end load as a Bernoulli beam model due to the consistent underlying mechanical behaviors. Thus, the displacement at the free end is monitored and compared against the exact displacement computed from Bernoulli beam theory. In this case, the exact maximum displacement at the free end is [u.sub.z] = 0.432 in [13]. We use both quadratic and cubic B-splines with maximal smoothness to model the beam and compare the results against Abaqus and Nastran in terms of accuracy per degree-of-freedom. We also employ all three quadrature schemes QP1, QPO, QNU.

The results are shown in Figure 2. In all cases, the blended RM5/KL3 elements converge fastest. This is because the rotational degrees-of-freedom are only used to impose the clamped boundary conditions on the left end. Everywhere else only three degrees-of-freedom per node are used, essentially cutting the number of degrees-of-freedom in half. This gain in efficiency is a direct consequence of smoothness enabling the use of the KL theory. It is also evident that higher-order smoothness is beneficial. For example, for p = 3 and the QNU quadrature scheme, the IGA results are far superior to the Abaqus and Nastran shell elements even in the context of this simple linear elastic problem.

CURVED CANTILEVER BEAM

Figure 3 shows a schematic for a curved cantilever beam. We will use this problem to illustrate the importance of using exact geometric quantities in the isogeometric blended shell formulation. We will compare the solution computed using the exact formulation described in Section to one computed using the geometric approximations described in Section (denoted Approx KL3 in the plots). Note that these geometric approximations are the same as those described in the first blended isogeometric shell paper [1]. The beam has a radius R= 10 in, a width b = 1 in, and the thickness t varies. Young's modulus and Poisson's ratio are E = 1.0 x [10.sup.3] psi and v = 0, respectively. The curved beam is clamped at one end and subjected to a traction, [q.sub.x] = O.1[t.sup.3] lbf/i[n.sup.4], at the other end. An analytical solution [2] based on Bernoulli beam theory gives a value of approximately 0.942 in for the radial displacement at the free end. The beam is modeled with quadratic B-splines of maximal smoothness. We also employ all three quadrature schemes QP1, QPO, QNU.

The results are shown in Figure 4. This problem suffers from both shear and membrane locking. RM5/KL3 is essentially shear locking free but still suffers from membrane locking as shown in Figure 4a. However, moving to the left, the reduced quadrature schemes alleviate locking producing locking free solutions for RM5/KL3 and RM5/RM5 in Figure 4c. Notice, however, that RM5/Approx KL3 remains inaccurate for all quadrature schemes. This is due to the geometric error in the underlying approximations which are employed in the shell formulation. These cannot be overcome through reduced quadrature schemes.

SCORDELIS-LO ROOF

The Scordelis-Lo roof problem, shown in Figure 5, tests a shell element's ability to handle combined membrane and bending. An 80[degrees] arc of a cylinder with radius R = 300 in, length L = 600 in, and thickness t = 3 in is supported on each end by a rigid diaphragm. It is then loaded with its own weight. The material has Young's Modulus E= 3.0 x [10.sup.6] psi, Poisson's ratio v = 0.3, and mass density [rho] = 5.3971 x [10.sup.-4] lbf. [s.sup.2]/[in.sup.4]. The acceleration of gravity is g = 386 in/[s.sup.2]. Only one quarter of the geometry is modeled due to the symmetry. We use both quadratic and cubic NURBS with maximal smoothness to model the shell and compare the results against Abaqus and Nastran in terms of accuracy per degree-of-freedom. We also employ all three quadrature schemes QP1, QP0, QNU.

The displacement on the free edge at L/2 is monitored and compared against a reference displacement of [u.sub.ref],= 3.59 in, based on Reissner-Mindlin shear deformable shell theory [14]. The results are shown in Figure 6. As before, the superiority of the blended isogeometric shells is evident. In this case, locking effects are less pronounced and accuracy improvements due to QNU are not as obvious. However, it should be noted that QNU is essentially a one-point quadrature scheme away from the shell boundaries, so the efficiency gains are significant.

CAR HOOD INNER

In this example, we demonstrate the potential for integrated design/analysis coupled with the accuracy gains provided by the smooth isogeometric basis. As shown in Figure 7, both a BREP and T-spline [15, 16, 17] model of a car hood inner were created. The BREP model was used for FEA where the geometry was cleaned up and meshed before a vibration calculation was performed with Nastran. The BREP models consists of 2090 trimmed NURBS patches. The T-spline, on the other hand, was directly analyzed without any geometry cleanup or mesh generation. In this example, Young's modulus E = 3.05 x [10.sup.7] psi, Poisson's ratio v = 0.3, mass density [rho] = 7.327 x [10.sup.-4] lbf. [s.sup.2]/[in.sup.4], and thickness t = 2.362 x [10.sup.-2] in. Bicubic isogeometric RM6 shells were used in the isogeometric vibration analysis. Frequencies up to 100 Hz were calculated, as shown in Figure 8. For all modes, the T-spline results match the NASTRAN results.

CONCLUSION

In this contribution, a goemetrically exact isogeometric blended shell formulation is presented. Several benchmarks are computed which include comparisons with standard FEA shell elements. The superiority of IGA is demonstrated in all cases. A complex car hood inner was then modeled both for FEA and IGA. In the case of IGA, no geometry cleanup or mesh generation is required. It is shown that the vibratory modes computed by FEA and IGA match.

REFERENCES

[1.] Benson D., Hartmann S., Bazilevs Y., Hsu M.-C, Hughes T., Blended isogeometric shells, Computer Methods in Applied Mechanics and Engineering 255 (0) (2013) 133 - 146.

[2.] Echter R., Oesterle B., Bischoff M., A hierarchic family of isogeometric shell finite elements, Computer Methods in Applied Mechanics and Engineering 254 (2013) 170 - 180.

[3.] Oesterle B., Ramm E., Bischoff M., A shear deformable, rotation-free isogeometric shell formulation, Computer Methods in Applied Mechanics and Engineering 307 (2016) 235 - 255.

[4.] Kiendl J., Bletzinger K.-U., Linhard J., Wchner R., Isogeometric shell analysis with kirchhofflove elements, Computer Methods in Applied Mechanics and Engineering 198 (4952) (2009) 3902 - 3914.

[5.] Kiendl J., Hsu M.-C., Wu M. C., Reali A., Isogeometric kirchhoff-love shell formulations for general hyperelastic materials, Computer Methods in Applied Mechanics and Engineering 291 (2015) 280-303.

[6.] Kiendl J., Bazilevs Y., Hsu M.-C., Wchner R., Bletzinger K.-U., The bending strip method for isogeometric analysis of kirchhofflove shell structures comprised of multiple patches, Computer Methods in Applied Mechanics and Engineering 199 (3740) (2010) 2403 - 2416.

[7.] Benson D., Bazilevs Y., Hsu M.-C., Hughes T., A large deformation, rotation-free, isogeometric shell, Computer Methods in Applied Mechanics and Engineering 200 (1316) (2011) 1367 - 1378.

[8.] Thomas D. C., Scott M. A., Evans J. A., Tew K, Evans E. J., Bezier projection: a unified approach for local projection and quadrature-free refinement and coarsening of NURBS and T-splines with particular application to isogeometric design and analysis, Computer Methods in Applied Mechanics and Engineering 284 (2015) 55-105.

[9.] Hughes T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola, NY, 2000.

[10.] Adam C., Bouabdallah S., Zarroug M., Maitournam H., Improved numerical integration for locking treatment in isogeometric structural elements, part i: Beams, Computer Methods in Applied Mechanics and Engineering 279 (2014) 1-28.

[11.] Adam C., Hughes T., Bouabdallah S., Zarroug M., Maitournam H, Selective and reduced numerical integrations for nurbs-based isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 284 (2015) 732-761.

[12.] Adam C., Bouabdallah S., Zarroug M., Maitournam H., Improved numerical integration for locking treatment in isogeometric structural elements. part ii: Plates and shells, Computer Methods in Applied Mechanics and Engineering 284 (2015) 106-137.

[13.] Macneal R. H., Harder R. L., A proposed standard set of problems to test finite element accuracy, Finite elements in Analysis and Design 1 (1) (1985)3-20.

[14.] Gallagher R. H., Finite elements for thin shells and curved members, John Wiley & Sons, 1976.

[15.] Bazilevs Y., Calo V. M., Cottrell J. A., Evans J. A., Hughes T. J. R., Lipton S., Scott M. A., Sederberg T. W., Isogeometric analysis using T-splines, Computer Methods in Applied Mechanics and Engineering 199 (5-8) (2010) 229-263.

[16.] Sederberg T. W., Zheng J., Bakenov A., Nasri A., T-splines and T-NURCCs, ACM Trans. Graph. 22 (2003) 477-484.

[17.] Scott M. A., Li X., Sederberg T. W., Hughes T. J. R., Local refinement of analysis-suitable T-splines, Computer Methods in Applied Mechanics and Engineering 213 (2012) 206 - 222.

Zhihui Zou, David Willoughby, and Michael A. Scott

Brigham Young University

Mohamed El-Essawi, Raghuraman Baskaran, and James Alanoly

Ford Motor Company

Printer friendly Cite/link Email Feedback | |

Author: | Zou, Zhihui; Willoughby, David; Scott, Michael A.; El-Essawi, Mohamed; Baskaran, Raghuraman; Alanoly |
---|---|

Publication: | SAE International Journal of Passenger Cars - Mechanical Systems |

Article Type: | Report |

Date: | Jul 1, 2017 |

Words: | 3561 |

Previous Article: | Development of Prediction Method for Engine Compartment Water Level by Using Coupled Multibody and Fluid Dynamics. |

Next Article: | Predicted Machining Dynamics for Powertrain Machining. |

Topics: |