Elasticity Theory Applications And Numerics

Getting Stiffness Matrices to Play Nice

Most people approach elasticity numerics from the textbook direction. They learn the weak form, they learn FEM, and then they try to assemble a code and immediately hit a wall of singular matrices. I'd rather start from the other end. Let me explain how this actually works when you're trying to get a result before your supervisor starts asking questions.

The fundamental problem in elasticity numerics is that you're taking a continuous partial differential equation—the Navier-Cauchy system for linear elasticity—and converting it into a discrete algebraic system. That conversion process, the weak formulation, is where most mistakes happen. The strong form says div(sigma) + f = 0 in the domain, with boundary conditions on displacement or traction. But you can't just discretize that directly without getting nowhere fast. You multiply by a test function, integrate by parts, and arrive at the weak form where you work in H1 spaces instead of dealing with second derivatives. This isn't optional hand-waving. It's what makes the problem solvable numerically. Once you have the weak form, you pick a discretization. The finite element method dominates this space, and for good reason. You divide your domain into elements—tetrahedra or hexahedra in 3D, triangles or quads in 2D—and approximate the displacement field within each element using polynomial shape functions. Linear triangles are the easiest starting point. The displacement varies linearly across each element, the strain is constant within each element, and the stress is also constant. That constant strain assumption is both the strength and the limitation of linear elements. For most structural problems, quadratic elements give you dramatically better accuracy per degree of freedom. A 10-node tetrahedron or 8-node hexahedron captures curvature in the displacement field that a 4-node tet simply cannot. The tradeoff is computational cost. Quadratic elements roughly double the degrees of freedom and increase the bandwidth of your stiffness matrix. But in my experience, the accuracy improvement usually justifies it unless you're working with massive models where memory is genuinely constrained.

Here's something most introductory treatments don't emphasize enough: the choice of element type matters far more than people realize when you're dealing with incompressible or nearly incompressible materials. When Poisson's ratio approaches 0.5, the material becomes nearly incompressible. A standard displacement-based linear element will lock. You'll get a solution that's orders of magnitude too stiff, and you won't immediately know why because the code runs without errors. I encountered this specifically when modeling rubber-like materials in a hyperelastic setting. The displacements were essentially zero across the entire domain despite a reasonable load. Spent two days debugging before I realized the element formulation was the culprit. Switched to mixed formulation with pressure as an independent degree of freedom, and the problem solved immediately. This is the kind of thing that doesn't show up in a tutorial. Assembling the global stiffness matrix follows a straightforward but meticulous process. For each element, you compute the strain-displacement matrix B, multiply it by the constitutive matrix D—which encodes Young's modulus and Poisson's ratio for isotropic materials, or the full elastic tensor for anisotropic cases—and integrate over the element volume. The result is the element stiffness matrix k_e. Then you map those local degrees of freedom into the global assembly using the element connectivity array. This is the direct sum operation that builds your global K matrix. The assembly step is where most numerical errors sneak in through indexing mistakes. Double-check your connectivity arrays against a small test case where you can compute the answer by hand. The constitutive relationship itself deserves more attention than it gets. For isotropic linear elasticity, D has two independent constants. You can express them as E and nu, or as lambda and mu (the Lamé parameters), or as K and mu. The choice affects numerical behavior. Working in terms of lambda and mu is more stable when you're implementing things from scratch because mu is always the shear modulus and stays well-behaved. E and nu become numerically problematic near nu = 0.5 because lambda goes to infinity. If you're writing your own solver, stick with lambda and mu internally and convert for input and output only.

Boundary conditions are where even experienced practitioners make costly mistakes. Essential (Dirichlet) boundary conditions prescribe displacements. Natural (Neumann) boundary conditions prescribe tractions. The way you apply them changes the solution. For Dirichlet conditions, you can either modify the system to eliminate constrained degrees of freedom or use penalty methods. Elimination is more accurate but requires reordering your matrix. Penalty methods are simpler to implement but introduce a parameter that needs careful tuning. Too small and the constraints aren't enforced. Too large and your matrix becomes ill-conditioned. For production work, I use the elimination approach with a skyline solver or sparse direct solver like MUMPS. Solving the resulting linear system K*u = F is actually the easy part compared to everything else. The matrix is symmetric positive definite for well-posed elasticity problems with proper boundary conditions. That means conjugate gradient is theoretically applicable, but the convergence rate depends entirely on your preconditioner. In practice, a sparse direct solver is often faster for moderate-sized problems because the overhead of iterative methods doesn't pay off until you're dealing with millions of degrees of freedom. For large-scale problems, an algebraic multigrid preconditioner paired with GMRES or CG will handle 10 million DOFs in reasonable time on a decent workstation. Validation is non-negotiable. Before you trust any elasticity simulation, run it against at least one analytical solution. The classic benchmarks are good for a reason: a single edge-cracked plate under tension, a pressurized thick-walled cylinder (Lamé problem), a simply supported beam under uniform load. If your FEM result doesn't match the analytical solution within an acceptable tolerance on a refined mesh, everything else you do with that code is suspect. I've seen teams skip this step and spend weeks chasing bugs that turned out to be fundamental formulation errors.

Mesh convergence is another area where people rush. A mesh that's too coarse will give you optimistic stress values in regions of high gradient. Stress is a derived quantity in FEM—it's computed from the displacement gradient, which is one derivative less smooth than the displacement itself. That means stress convergence is always slower than displacement convergence. If your stress at a point of interest doesn't stabilize as you refine the mesh, you need to keep refining or switch to a formulation that recovers smoother stresses. Superconvergent patch recovery is a standard technique for this. It averages the stress at a node from all surrounding elements, which tends to produce values that converge faster than the raw element stresses. Nonlinear elasticity introduces additional complexity. Material nonlinearity means D is no longer constant—it depends on the current strain state. Geometric nonlinearity means you're working with finite strains and the equilibrium equations must be written in the deformed configuration. Both require an iterative solution procedure. Newton-Raphson is the standard approach. Each iteration involves computing a tangent stiffness matrix, solving for the displacement increment, and updating the configuration. The tangent stiffness must be consistent with the constitutive model. Using an approximate tangent instead of the exact consistent tangent will slow convergence dramatically or cause divergence entirely. I learned this the hard way when switching from a neo-Hookean to a Mooney-Rivlin material model. The code converged in five iterations with the consistent tangent and didn't converge at all with the secant approximation. Contact problems are perhaps the hardest category in practical elasticity numerics. When two bodies touch, the boundary conditions change depending on whether the surfaces are in contact, separating, or sliding. This is a free boundary problem. Penalty-based contact formulations are the most common approach in commercial software. You allow a small amount of penetration and apply a restoring force proportional to the penetration depth. The penalty parameter needs to be large enough to prevent excessive penetration but small enough to avoid ill-conditioning. There's no universal value. You tune it based on your material stiffness and mesh size. Augmented Lagrangian methods are more robust but more complex to implement. For most engineering applications, a carefully tuned penalty formulation gives results that are good enough.

The biggest practical bottleneck in elasticity numerics is usually not the solver—it's the preprocessing. Generating a quality mesh for a complex 3D geometry takes more time than running the simulation. Automatic mesh generators produce reasonable meshes, but they struggle with regions of high stress gradient. You need to identify those regions beforehand and locally refine the mesh. Manual mesh control points are worth the effort. A mesh that's uniformly fine everywhere is almost never the most efficient approach. Targeted refinement around holes, notches, and contact zones typically gives you the same accuracy with a fraction of the degrees of freedom. Post-processing deserves the same care as the solution itself. Principal stresses, von Mises stress, strain energy density—these are all derived quantities that tell you different things about the structural response. Von Mises stress is useful for ductile materials under yielding criteria. Principal stresses matter more for brittle materials where tensile failure is the concern. Neither is universally "correct." Pick the metric that matches your failure criterion and report it consistently. Comparing von Mises stress against a tensile strength limit is a common mistake that gives non-conservative results. There's also the question of whether FEM is the right tool for your problem. Boundary element methods reduce the dimensionality by one because they only discretize the boundary. For infinite domain problems like crack propagation or soil-structure interaction, BEM can be significantly more efficient. The tradeoff is that BEM produces dense matrices instead of sparse ones, and the fundamental solutions required for anisotropic or heterogeneous materials are not always available. For homogeneous isotropic elasticity in unbounded domains, BEM is elegant and accurate. For everything else, FEM is the safer default.

If you're looking for implementation references, the deal.II library is probably the most comprehensive open-source framework for elasticity numerics. It handles adaptive mesh refinement, multiple element types, and various constitutive models out of the box. The documentation is thorough, though the learning curve is steep. For a more minimal starting point, Elmer offers a good balance between accessibility and capability. If you just need to prototype something quickly and don't care about performance, Fenics simplifies the weak form specification considerably, though you lose some control over the linear algebra behind the scenes. The field moves slowly in terms of fundamental methods. The mathematics hasn't changed much in fifty years. What changes is the scale of problems you can solve and the sophistication of the material models you can include. The practical knowledge—the stuff that actually determines whether your simulation produces useful results—is still largely transmitted through osmosis and painful experience rather than documentation.