Finite Element Method for Shells

Thin-walled structures, or shells, are ubiquitous in nature and engineering, from bird wings and eggshells to aircraft fuselages and skyscrapers. Their remarkable combination of lightweight form and structural strength comes from their curvature, which allows them to carry loads through a complex interplay of in-plane stretching and out-of-plane bending. Capturing this behavior computationally presents a significant challenge: how can we create a digital model that is both accurate and efficient? The Finite Element Method (FEM) provides a powerful answer, offering a framework to translate the physics of shells into a language that computers can solve.

This article delves into the core of modern shell analysis, focusing on one of the most successful and widely used approaches. The primary knowledge gap it addresses is the gap between simply using shell elements in software and truly understanding how they work, including their inherent limitations and the ingenious solutions developed to overcome them. The reader will gain a robust conceptual understanding of the method, from its fundamental building blocks to its application in solving sophisticated engineering problems.

First, we will explore the ​​Principles and Mechanisms​​ behind the popular "degenerated solid" shell element. This includes its kinematic foundation, the critical problem of numerical locking that can plague thin shell analysis, and the clever techniques devised to cure it. We will then examine ​​Applications and Interdisciplinary Connections​​, demonstrating how these verified and robust elements are used to tackle real-world challenges like structural buckling, contact mechanics, and the design of advanced composite materials.

In our journey to understand the world through computation, we often seek a kind of elegant simplicity. We want to capture the essence of a physical object without getting lost in every minute detail. Imagine trying to describe the majestic curve of a bird's wing or the taut surface of a drum. These are shells—structures defined by their thinness, where the behavior is a beautiful and complex dance between stretching and bending. How can we teach a computer to understand this dance? The answer, born of remarkable ingenuity, is one of the most powerful ideas in computational mechanics.

One way to describe a shell would be to start from scratch, writing down complicated two-dimensional equations full of arcane terms for curvature. This is the classical path, but it is a difficult one, paved with mathematical brambles. A more intuitive and, as it turns out, more versatile approach is to start with something we already understand very well: a simple, three-dimensional solid brick. Then, we "degenerate" it.

This is the ​​degenerated solid approach​​, an idea pioneered by engineers like Ahmad, Irons, and Zienkiewicz. Imagine your solid brick is a reference block, like a Rubik's cube, defined by three coordinates $(\xi, \eta, \zeta)$ that each run from -1 to 1. We'll decide that the coordinates $\xi$ and $\eta$ will trace the surface of our shell, and $\zeta$ will run through its thickness. Now, we invent a mapping that takes any point in this simple reference cube and places it in 3D space to form our shell element. The position of any point $\mathbf{x}$ in the shell is given by a wonderfully simple rule:

Let's break this down. $\mathbf{x}_0(\xi, \eta)$ is the position of the shell’s ​​midsurface​​—think of it as the middle page of a book. The term $t$ is the shell's thickness. The coordinate $\zeta$ is our dimensionless thickness parameter; $\zeta=-1$ is the bottom surface, $\zeta=0$ is the midsurface, and $\zeta=+1$ is the top surface.

The "magic" is in the vector $\mathbf{d}(\xi, \eta)$. This is called the ​​director vector​​. It's a vector defined at every point on the midsurface that initially points straight through the thickness, normal to the surface. It "directs" the line of points that makes up the shell's thickness. The beauty of this formulation is that we can describe the entire shell's geometry and deformation just by knowing what its midsurface is doing (how $\mathbf{x}_0$ moves) and how its director is behaving (how $\mathbf{d}$ rotates). Instead of tracking every point in a 3D solid, we only need to track a 2D surface and a field of vectors attached to it.

Using the isoparametric concept, we interpolate these quantities from values at the element's nodes (its corners). For a four-node element, we use four shape functions $N_i(\xi, \eta)$ to blend the nodal positions $\mathbf{x}_i$ and nodal directors $\mathbf{d}_i$:

This equation is the heart of the degenerated shell element. It tells us that what seems like a complex 3D object can be built from simple 2D building blocks and a director. By allowing the director $\mathbf{d}$ to rotate independently of the midsurface, this formulation naturally embodies the ​​Mindlin-Reissner​​ kinematic theory, which accounts for the sliding of layers relative to one another—the transverse shear deformation—that is important in moderately thick shells. If we were to force the director to always stay perfectly perpendicular to the deformed midsurface, we would recover the older, more restrictive ​​Kirchhoff-Love​​ theory.

So, what does a node of our new shell element need to "know" to describe its motion? A node in a standard 3D solid element only knows how to move; it has three ​​degrees of freedom (DOFs)​​: translation in the $x$, $y$, and $z$ directions. Our shell element node is more sophisticated. It can translate (3 DOFs), but it must also describe the orientation of its director vector.

In classical shell theory, the director's orientation is described by its rotation about two axes in the plane of the shell, say $\theta_x$ and $\theta_y$. These two rotations describe bending and shearing beautifully. So, a classical shell node has $3$ translations + $2$ rotations = $5$ DOFs.

But what about the third possible rotation, the "spinning" of the director about its own axis? This is called the ​​drilling rotation​​, $\theta_z$. For a single, flat plate, this rotation does absolutely nothing. It produces no strain and therefore no strain energy. It is a ghost degree of freedom. So, why would we even consider it?

The reason is a practical one. Imagine you are building a structure where two shell elements meet at a sharp angle. If each node only carries 5 DOFs, how do you correctly transfer moments between the elements? The problem becomes much simpler if every node in our system speaks the same language: a full set of 3 translations and 3 rotations. Therefore, most modern degenerated shell elements are implemented with ​​6 DOFs per node​​ ($u_x, u_y, u_z, \theta_x, \theta_y, \theta_z$).

But we just said the drilling rotation has no physical stiffness! By adding it, we have created a "zero-energy mode"—a way for the model to deform without any resistance, leading to a singular stiffness matrix and a failed analysis. To solve this, a tiny amount of artificial stiffness is associated with the drilling rotation, just enough to stabilize the element and prevent this "pinwheeling" behavior at element junctions, but not enough to affect the real physical response. This is a perfect example of the pragmatic compromises made in computational engineering, balancing theoretical purity with practical robustness.

We've built what seems to be a simple, powerful tool. But this simplicity hides a dangerous flaw. Under certain conditions—specifically when the shell becomes very thin—these elements can become pathologically, absurdly stiff. They refuse to bend, as if "locked" in place. This is not a real physical effect; a thin sheet of metal or paper bends very easily. ​​Locking​​ is a pure pathology of the discretization—a failure of our simple polynomial shape functions to capture the subtle physics of thin structures.

Shear Locking: When Your Deck of Cards Turns into a Brick

Imagine a thick deck of cards. When you bend it, the cards slide past each other. This is ​​transverse shear​​. Our Mindlin-Reissner element is designed to capture this. Now, imagine a single, very thin sheet of paper. When you bend it, it does so without any significant shearing; the fibers through its thickness remain perpendicular to the bent surface (the Kirchhoff-Love limit).

Here's the problem: as our shell element becomes very thin, its shear stiffness (proportional to thickness $t$) becomes enormous compared to its bending stiffness (proportional to $t^3$). The physics demands that the shear strain must go to zero to keep the total energy finite. But our simple bilinear element is not "smart" enough to represent a state of pure bending without simultaneously producing small, spurious shear strains throughout its volume. The element tries to bend, but this inadvertently creates parasitic shear. The huge shear stiffness penalizes this parasitic shear so severely that it effectively prevents the element from bending at all. It "locks." The deck of cards has effectively turned back into a solid, unbendable brick.

Membrane Locking: The Curse of the Curve

Another, more subtle form of locking appears when we model ​​curved​​ shells. Consider a shallow arch. The physics of curved structures involves an intricate coupling between out-of-plane bending and in-plane "membrane" stretching. The shell's curvature means that a simple transverse displacement $w$ will induce in-plane membrane strains, according to relations like $\varepsilon \sim u_{,\alpha} + b w$, where $b$ is the curvature.

An elegant deformation mode for a curved shell is "inextensional bending," where it bends without stretching its midsurface. To do this, the in-plane displacements $u$ must adjust in a very specific way to perfectly cancel out the membrane strain caused by the transverse displacement $w$. Once again, our simple bilinear shape functions are often too crude to manage this delicate choreography. In trying to bend, a locked element inevitably generates spurious membrane strains.

Just like with shear, the membrane stiffness (scaling with $t$) is vastly greater than the bending stiffness (scaling with $t^3$) for a thin shell. This spurious stretching is met with immense resistance, and the element wrongly reports a huge stiffness, locking up the bending motion. This has disastrous consequences in, for example, buckling analysis, where a locked element will drastically overpredict the load a shell can carry before it snaps through.

This story of locking sounds like a tragedy of good intentions. But the history of the finite element method is one of clever solutions triumphing over numerical gremlins. Engineers have developed brilliant strategies to cure locking, which can be thought of as forms of "cheating smartly."

The Art of Not Looking Too Closely: Reduced Integration

One of the earliest and simplest cures for shear locking is ​​selective reduced integration​​. The reasoning goes like this: if the element is producing spurious shear strains that cause trouble, maybe we should just be less picky about measuring them. Instead of calculating the shear energy at multiple locations (e.g., four Gauss points) inside the element, which is called "full integration," we calculate it at only a single point, usually the element's center.

It turns out that for many simple element shapes, the spurious shear strains happen to be zero at the element's center. By sampling the shear energy only at this "sweet spot," we effectively make the element blind to its own parasitic shear, freeing it to bend correctly.

This is a powerful trick, but it's a bit of a devil's bargain. By sampling the energy at only one point, the element also becomes blind to certain real deformation modes. These are non-physical, wobbly motions called ​​hourglass modes​​, which have zero energy and can destroy a simulation. Thus, elements using reduced integration must often be paired with "hourglass stabilization" schemes, which add back just enough artificial stiffness to control these wobbly modes. It's a testament to the fact that in engineering, there is rarely a free lunch.

Assuming a Better Reality: Enhanced and Assumed Strains

A more elegant and robust solution is to attack the problem at its source. If the strains derived from the displacement interpolation are flawed, why not just... assume a better strain field? This is the core idea behind ​​assumed strain methods​​, such as the Assumed Natural Strain (ANS) and Mixed Interpolation of Tensorial Components (MITC) families of elements.

In these methods, we abandon the strain field calculated directly from the displacements. Instead, we construct a new, independent strain field inside the element. This assumed field is carefully designed to be simple enough to avoid the spurious components that cause locking, but rich enough to represent the true physical states, like constant strain and pure bending.

For instance, in an ​​Assumed Natural Strain (ANS)​​ formulation for shear locking, one might sample the shear strains at a few key locations (like the midpoints of the element's edges) and then use a simple interpolation to define the "assumed" shear strain field everywhere else. This constructed field is guaranteed to be well-behaved and can be shown to correctly pass fundamental consistency checks, like the "patch test," which ensures the element can exactly represent a state of constant strain when it should. These methods represent a deeper understanding of the problem, moving from a simple numerical trick to a theoretically sound reformulation of the element itself. They are the reason why modern finite element software can reliably and accurately predict the behavior of the most complex shell structures, from aircraft fuselages to civil engineering marvels.

We have spent some time learning the principles behind the finite element method for shells—the "grammar," if you will, of how we describe these elegant structures in a language a computer can understand. But knowing grammar is one thing; writing poetry is another. Now, we are going to see the poetry. We will explore how these abstract ideas breathe life into solutions for real-world problems, transforming our computer from a glorified calculator into a veritable crystal ball for engineers and scientists.

The journey from a set of mathematical equations to a reliable prediction of a skyscraper's response to an earthquake or a car's behavior in a crash is not one of blind faith. It is a craft, a science, and an art. It is in the application of these principles that we find their true power and beauty.

Before we simulate a jumbo jet, we must be certain that our digital building blocks aren't made of sand. How do we know the computer isn't telling us a convenient, but dangerously wrong, story? We can't just ask it. Instead, we must be clever. We must design a "digital obstacle course"—a series of rigorous verification tests grounded in the fundamental laws of physics.

First, we demand that our shell element pass the ​​patch test​​. Imagine a small patch of elements. If we apply boundary conditions that should produce a simple, constant state of stretching or a pure, constant curvature, does the element assembly reproduce this state exactly? If it can't handle these simplest of all possible worlds, it has no hope of being correct in a complex one. This is the computational equivalent of checking if a new calculator knows that $1+1=2$.

Next, we check for ​​rigid-body invariance​​. If we take an object and simply move it or rotate it without deforming it, it should experience no internal stresses or strains. Our element must agree. We can check this by examining the stiffness matrix of an unconstrained element; it must have exactly six zero-energy modes in three-dimensional space, corresponding to three translations and three rotations. If it has fewer, it means the element locks up and resists rigid motion; if it has more, it means the element has extra, unphysical floppy modes.

Finally, we test for ​​convergence​​. As we refine our mesh, using more and smaller elements to represent our structure, does our computed answer get closer to the true, exact answer? A well-behaved element does, and it does so at a predictable rate that depends on the complexity of its underlying mathematical basis. By performing these checks, we build confidence that our numerical microscope is not flawed, and that the images it shows us bear a genuine resemblance to reality.

Even with an element that passes all our verification tests, strange gremlins can appear. These are numerical pathologies, artifacts of discretization that can corrupt our solution. The most notorious of these is "locking."

Consider the paradox of a thin shell. It derives its strength from its curvature, but its utility from its flexibility in bending. Now, what happens when we try to model this with simple, low-order elements? As a shell gets very, very thin, a phenomenon called ​​transverse shear locking​​ can occur. The physical bending energy of a thin shell scales with the cube of its thickness, $t^3$, while its transverse shear energy scales linearly with thickness, $t$. A poorly designed element, when forced to bend, can generate a large amount of spurious, non-physical shear energy. In the thin limit, this parasitic shear energy completely dominates the true bending energy, making the element artificially stiff—it "locks." This isn't just a numerical curiosity; it can lead to catastrophic underestimation of deflections and a dangerous overestimation of a structure's buckling load. A bridge that the computer declares safe might, in reality, be poised for collapse.

This pathology is beautifully illustrated in the "pinched cylinder" benchmark, a famous torture test for shell elements. Here, both shear locking and its cousin, ​​membrane locking​​ (which pollutes bending with spurious membrane energy), can conspire to give wildly incorrect results.

How do we exorcise these gremlins? The solutions are a testament to the ingenuity of computational scientists. One common technique is ​​reduced integration​​, where we cleverly compute the element's stiffness by sampling the internal strains at fewer, more strategic points. By avoiding the points that are most sensitive to generating spurious energy, we can often alleviate locking. However, this fix can introduce a new gremlin: ​​hourglassing​​. The element can become too flexible, deforming in floppy, zero-energy patterns that resemble an hourglass, contaminating the solution with checkerboard-like modes. The cure? We add a tiny amount of artificial "hourglass control" stiffness, just enough to suppress the unphysical modes without affecting the real physics. It is a beautiful and delicate dance of trade-offs, a true art form within the science of simulation.

With our digital tools verified and our numerical gremlins tamed, we can now turn our attention to problems that were, for most of history, utterly intractable.

The Grace and Violence of Buckling

Think of the sudden snap of a flexible ruler compressed between your hands, or the crinkling of a soda can as you step on it. This is buckling—a sudden, often catastrophic loss of stability. Shells are particularly susceptible to this behavior. Predicting it is one of the most critical and challenging tasks in structural engineering.

Here, we discover another profound lesson: not only does the element formulation matter, but the accuracy of the geometry itself is paramount. Consider the "snap-through" of a shallow spherical cap under pressure. If we model this curved shell with a coarse mesh of flat facets, we are fundamentally misrepresenting its geometry. One might intuitively think that making the shell "flatter" would make it seem weaker. The reality is more subtle. The flatter, discretized shell must develop higher compressive membrane stresses to resist the same external pressure. This higher pre-stress amplifies the geometric stiffness contribution that leads to instability, causing the model to predict buckling at a lower load than is correct. A small error in geometry can lead to a dangerously non-conservative prediction of failure. The lesson is clear: "garbage in, garbage out" applies not just to data, but to geometry itself. The solution involves both improving the model—using higher-order elements or isogeometric analysis that can perfectly represent the curvature—and using powerful numerical algorithms that can trace the complex, snapping equilibrium path beyond the stability limit.

When Worlds Collide: The Challenge of Contact

So much of the world is governed by things touching each other. A car crash, a metal stamping press forming a fender, a surgeon inserting a prosthetic joint—these are all contact problems. The governing physics is deceptively simple to state: two objects cannot occupy the same space at the same time. But how do we teach a computer this fundamental rule?

The finite element method provides an elegant framework. We designate one surface as the "master" and another as the "slave." For any point on the slave surface, we find the closest point on the master surface. The signed distance between them, measured along the normal to the master surface, is defined as the "gap," $g_n$. The non-penetration condition is then simply the inequality $g_n \ge 0$. Computational contact algorithms are sophisticated methods for finding the forces required to enforce this simple geometric constraint at millions of points simultaneously. This seemingly simple idea unlocks the ability to simulate and design for an enormous class of complex, real-world phenomena.

Few real-world structures are made of a single, monolithic shell. They are almost always complex assemblies of different components, and they are increasingly built from advanced, multi-layered materials.

Building with Digital LEGOs

Think of an airplane. The skin is a shell, but the wings are reinforced with internal spars and ribs (beams), and the landing gear is a complex assembly of solid components. To analyze the whole system, we need to connect different types of elements together. How do you weld a 2D shell element to a 1D beam element in a simulation?

You must ensure kinematic compatibility: their translations must match, and their rotations must be consistent. But here lies another subtlety. A standard shell element has five degrees of freedom per node: three translations and two rotations that bend the shell. It has no physical stiffness associated with a "drilling" rotation about its own normal. A beam element, however, does have a meaningful torsional stiffness. If we blindly force all three rotational degrees of freedom to be identical at the connection, we would be artificially constraining the beam's torsion, leading to an incorrect result. The art lies in understanding the physics of each idealization. The correct approach is to tie the common degrees of freedom—the translations and the two bending rotations—while leaving the unshared "drilling/torsion" degree of freedom free to rotate. This allows us to build vast, complex models from a library of specialized "digital LEGOs."

The Fabric of Modern Engineering: Composite Materials

From the chassis of a Formula 1 car to the wings of a Boeing 787, advanced composite materials have revolutionized engineering. These materials, made of layers of fibers embedded in a matrix, are incredibly strong and lightweight. Their properties, however, are directional—a carbon fiber laminate is far stiffer along the fiber direction than across it.

The shell element provides the perfect framework for modeling these anisotropic materials. At each integration point within our shell, we can define a multi-layered stack. For each layer (or "lamina"), we provide its thickness, its material properties, and, crucially, its fiber orientation angle. The finite element software then performs the necessary tensor transformations to compute the layer's stiffness in the element's local coordinate system. By integrating these properties through the thickness, the shell "homogenizes" the complex layered behavior into an equivalent, but highly anisotropic, shell stiffness. This allows us to design and analyze structures with properties tailored for specific performance goals. Once again, we find that the "drilling" degree of freedom, while not strictly necessary from a continuum mechanics standpoint, can be added with a numerically small, vanishing penalty to improve element performance, a further example of the artful compromises made in computational mechanics.

Shell models are, by nature, approximations that reduce a 3D reality to a 2D surface. But can we use this simplified model to peer back into the full three-dimensional world, especially in critical regions where the 2D approximation might break down?

The Hidden Stresses that Kill Composites

The Achilles' heel of composite laminates is often not failure of the fibers themselves, but ​​delamination​​—the peeling apart of adjacent layers. This failure mode is driven by "interlaminar stresses," the out-of-plane shear ($\tau_{xz}, \tau_{yz}$) and normal ($\sigma_{zz}$) stresses that exist between the layers. A 2D shell model does not compute these directly.

However, we can "recover" them in a post-processing step. By taking the in-plane stresses calculated by the shell element and integrating the fundamental 3D equations of equilibrium through the laminate's thickness, we can solve for the distribution of these hidden interlaminar stresses. This is an incredibly powerful technique. To trust this recovered 3D stress field, we must subject it to its own verification checklist: it must satisfy 3D equilibrium, it must meet the traction-free boundary conditions on all free surfaces, and its integral must be consistent with the shear force resultants from the parent shell model. This allows engineers to use computationally efficient shell models to pinpoint regions at high risk of delamination, a critical step in the safety certification of composite aircraft and vehicles.

Sharpening the Picture: Stress Recovery and Adaptive Meshing

A raw finite element solution provides a stress field that is typically discontinuous and "jagged" across element boundaries. While often highly accurate at specific integration points (the Gauss points) inside the element, this raw output can be difficult to interpret. The Zienkiewicz-Zhu (ZZ) recovery technique provides a way to create a smooth, continuous, and even more accurate stress field from this discrete data. The idea is to fit a smooth polynomial patch over the highly accurate Gauss point values in the vicinity of each node. This process, often called Superconvergent Patch Recovery (SPR), gives us a much cleaner and more reliable picture of the stress distribution.

This isn't just about creating aesthetically pleasing contour plots. The difference between the raw, jagged stress field and the smooth, recovered one serves as a brilliant a posteriori error estimator. In regions where this difference is large, our solution is likely inaccurate. We can use this information to drive an adaptive meshing process, where the computer automatically refines the mesh in areas of high estimated error and re-runs the analysis, repeating the cycle until a desired level of accuracy is achieved everywhere.

We have seen that the finite element method for shells is far more than a black box. It is a rich and fascinating field, a delicate interplay between physics, mathematics, and computer science. From the foundational checks that build our trust in the method, to the clever tricks used to tame its pathologies, to the powerful ways it connects to different engineering disciplines and materials, we find a consistent theme: simple, powerful ideas can be composed and extended to solve problems of incredible complexity. The journey from a simple concept like a thin sheet of paper to the full virtual prototype of a next-generation aircraft is a long one, but it is paved with the beautiful and insightful principles we have explored.