Uniformly Distributed Load on a Simply Supported Beam — FEM Verification Benchmark

It's the most fundamental verification problem for FEM. For a simply supported beam under a uniformly distributed load $w$, an exact deflection solution is obtained from Euler-Bernoulli beam theory. The maximum deflection at the center is

The purpose of this benchmark is to compare this theoretical solution with the numerical solution to confirm the accuracy of the element type. Particularly important is that when using first-order elements (linear elements), shear locking occurs, causing the beam to behave stiffer than it actually is, and the deflection is underestimated. Quantitatively understanding this phenomenon and verifying the improvement effect of reduced integration or the B-bar method are practical points.

Yes. Precisely because it's simple, it's so useful. It's included in the NAFEMS benchmark collection and is also recommended as the first step for Code Verification (verification of mathematical correctness of the code) in ASME V&V 10-2006. When introducing a new solver, confirming agreement with the theoretical solution within 0.1% using this problem first is the industry's standard procedure. For example, even in the structural analysis departments of automotive manufacturers, this benchmark is run as a regression test every time the solver is updated.

Let's start from the governing equation of the Euler-Bernoulli beam. The deflection curve for a simply supported beam of length $L$ under a uniformly distributed load $w$ (N/m) is

where $x$ is the distance from the left end, $E$ is Young's modulus, and $I$ is the second moment of area. Substituting the center $x = L/2$ gives the maximum deflection.

Next, the bending moment distribution.

It is maximum at the center: $M_{\max} = wL^2/8$. From this, the maximum bending stress is

$Z = I/c$ is the section modulus, and $c$ is the distance from the neutral axis to the outermost fiber. For a rectangular cross-section $b \times h$, $I = bh^3/12$, $c = h/2$, $Z = bh^2/6$.

The shear force is

Maximum at the supports $V_{\max} = wL/2$, zero at the center. Since Euler-Bernoulli theory assumes that cross-sections remain plane and normal (ignoring shear deformation), its accuracy is high for thin beams with a large span/depth ratio $L/h$.

Timoshenko beam theory additionally considers shear deformation. The center deflection is

The second term is the deflection increment due to shear deformation. $\kappa$ is the shear correction factor (5/6 for rectangular cross-section), $G$ is the shear modulus $G = E/\{2(1+\nu)\}$.

Euler-Bernoulli is sufficient for $L/h \geq 20$. For deep beams with $L/h < 10$, Timoshenko theory or 3D solid elements are needed. In benchmarks, typically $L/h = 50$ or so is used to verify under conditions where Euler-Bernoulli theory holds exactly.

The standard benchmark settings are as follows.

$L/h = 50$, so Euler-Bernoulli theory holds sufficiently. Second moment of area $I = bh^3/12 = 100 \times 20^3/12 = 66{,}667$ mm$^4$.

Here $Z = bh^2/6 = 100 \times 20^2/6 = 6{,}667$ mm$^3$.

Comparison results with the above parameters.

HEX8 with full integration shows about 17% error, and TET4 even shows 30% displacement error. That's the power of shear locking. On the other hand, HEX8 with reduced integration or HEX20 have errors within 0.3%. Clearly showing this difference is the significance of the simply supported beam benchmark.

To understand the essence, consider the limitations of the shape functions of first-order elements.

The displacement field of a first-order hexahedral element (HEX8) is linear in each direction. In pure bending, the cross-section rotates, causing horizontal displacements in opposite directions on the top and bottom surfaces. Theoretically, the shear strain $\gamma_{xy} = 0$ should hold at this time, but a linear displacement field cannot satisfy this condition.

In other words, parasitic shear strain that shouldn't exist is generated. The shear stress corresponding to this strain acts as extra stiffness, making the beam appear stiffer than it actually is, i.e., underestimating the deflection. This is shear locking.

Exactly. Refining the mesh improves the aspect ratio of individual elements and alleviates it, but convergence is very slow. Two elements in the thickness direction are far from enough; about 8 to 16 elements are needed. Since the computational cost becomes enormous, using locking countermeasure techniques is far more efficient in practice.

There are three main approaches.

1. Reduced Integration
Use reduced integration (1×1×1=1 point) instead of full integration (2×2×2=8 points for HEX8). Reducing integration points can eliminate the energy contribution of parasitic shear strain. However, the trade-off is the generation of hourglass mode, a zero-energy deformation mode. Even if it deforms like an hourglass, the energy is evaluated as zero.

2. B-bar Method (Selective Reduced Integration)
Only the volumetric change (dilation) component is calculated with reduced integration, while the deviatoric component maintains full integration. Effective for both volumetric locking and shear locking. Abaqus's C3D8R + hourglass control and Nastran's CHEXA reduced are based on this idea.

3. EAS Method (Enhanced Assumed Strain)
In addition to the usual displacement field, additional strain modes are assumed inside the element. Additional strain parameters are introduced as internal degrees of freedom and made invisible externally through static condensation. Ansys's SOLID185 (keyopt(2)=2) corresponds to this.

Artificial stiffness (hourglass stiffness) is added to suppress hourglass mode. If hourglass energy exceeds 5-10% of the total internal energy, it's a danger sign. In Abaqus, *SECTION CONTROLS, HOURGLASS=ENHANCED; in LS-DYNA, IHQ=6 (Belytschko-Bindeman) are proven settings. You can check the hourglass energy ratio with this benchmark problem and judge it as acceptable if it's below 1%.

According to FEM a priori error estimation (Cea's lemma), for $p$-th order elements, the energy norm error decreases as $O(h^p)$, and the $L^2$ norm error for displacement decreases as $O(h^{p+1})$.

That is, for linear elements ($p=1$), halving the element size $h$ reduces the displacement error to about $1/4$; for quadratic elements ($p=2$), to about $1/8$. Confirming whether this theoretical convergence rate is reproduced in actual mesh convergence tests is the core of Code Verification.

Specifically, analyze with three or more different mesh densities and check on a log-log plot whether the slope is close to the theoretical value $p+1$. If the slope deviates significantly from the theoretical value, suspect bugs, incorrect boundary conditions, or the influence of singular points.

The basic flow is as follows.

Step 1: Create Geometry
A rectangular solid of $L \times b \times h = 1000 \times 100 \times 20$ mm. You can also use a half-model ($L/2$) utilizing symmetry.

Step 2: Generate Mesh
Start with a coarse mesh (10 divisions in span direction, 2 divisions in thickness direction), then progressively refine. Recommended mesh convergence test: 4 stages: span direction 10→20→40→80.

Step 3: Define Material
Linear elastic. $E = 200{,}000$ MPa, $\nu = 0.3$.

Step 4: Set Boundary Conditions
Left end: $u_y = 0$ (vertical constraint) + $u_z = 0$ (out-of-plane displacement constraint)
Right end: $u_y = 0$ + $u_z = 0$
At one end: $u_x = 0$ (constrain axial rigid body motion)

Step 5: Apply Load
Surface pressure $p = w/b = 1.0/100 = 0.01$ MPa on the top surface. Or line load $w = 1.0$ N/mm.

Step 6: Solve & Post-process
Compare the maximum deflection and maximum stress at the center cross-section with the theoretical solution.

There are three common pitfalls.

1. Over-constraint
Setting $u_x = 0$ at both ends even though it's "simply supported." Axial elongation like thermal expansion is constrained, generating parasitic axial force. Set $u_x = 0$ at only one end.

2. Kink Generation Due to Point Constraint
If a support point is constrained with just a single point (single node) in solid elements, stress concentration occurs there. The correct approach is to constrain nodes along the support edge or face, or define a support line using MPC (Multi-Point Constraint).

3. Forgetting Out-of-Plane Constraint
In 3D models, $u_z = 0$ must be set in the out-of-plane direction to prevent rigid body rotation about the z-axis. Forgetting this leads to a singular stiffness matrix and the solver throws an error.

Here is a checklist for passing the benchmark.