Stress concentration at a circular hole

A plate in tension, solved in FEniCSx in 2D and 3D and checked against Kirsch and Peterson

FEniCSx

linear elasticity

validation

Linear-elastic plate with a central hole: FEniCSx plane-stress and solid models, validated against Kirsch’s solution and Peterson’s finite-width chart.

Published

October 10, 2026

A plate with a central hole, pulled along its length, is the textbook stress raiser: the hole lifts the stress at its edge to almost four times the remote value. The plate is solved in FEniCSx as a 2D plane-stress model and as a 3D solid, and both are checked against Kirsch’s closed-form solution and Peterson’s finite-width chart.

Figure 1: Normal stress \(\sigma_{xx}/\sigma_\infty\) on the 3D model, mirrored from the one-eighth symmetry model. Drag to rotate, scroll to zoom, hover for values.

\(K_t\) · FE, 2D plane stress 3.78 +1.2% vs Peterson

\(K_t\) · FE, 3D mid-plane 3.81 free surface 3.76

\(K_t\) · Peterson, finite width 3.73 \(d/W\) = 0.4

\(K_t\) · Kirsch, infinite plate 3.00 2D FE is +25.9% above

The model

Plate	400 × 200 × 10 mm, central hole of radius \(a\) = 40 mm (\(d/W\) = 0.4)
Material	linear elastic, \(E\) = 210 GPa, \(\nu\) = 0.3
Load	2 kN on each end, so \(\sigma_\infty = F/(W t)\) = 1 MPa
Symmetry	quarter plate in 2D, one eighth in 3D (planes \(x = 0\), \(y = 0\) and the mid-plane \(z = 0\))
Mesh	linear quadrilaterals (416 in 2D) and hexahedra (1,263 in 3D, 3 layers through the half thickness), 5 mm at the hole growing to 10 mm

Figure 2: The 2D quarter model: its mesh, the two symmetry edges, and the remote stress applied as a traction on the loaded end.

Analytical solution

Kirsch (1898) solved the infinite plate with an Airy stress function. With the remote stress \(\sigma_\infty\) along \(x\) and \(\theta\) measured from the load axis:

\[ \begin{aligned} \sigma_{rr} &= \frac{\sigma_\infty}{2}\left(1 - \frac{a^2}{r^2}\right) + \frac{\sigma_\infty}{2}\left(1 - \frac{4a^2}{r^2} + \frac{3a^4}{r^4}\right)\cos 2\theta \\ \sigma_{\theta\theta} &= \frac{\sigma_\infty}{2}\left(1 + \frac{a^2}{r^2}\right) - \frac{\sigma_\infty}{2}\left(1 + \frac{3a^4}{r^4}\right)\cos 2\theta \\ \sigma_{r\theta} &= -\frac{\sigma_\infty}{2}\left(1 + \frac{2a^2}{r^2} - \frac{3a^4}{r^4}\right)\sin 2\theta \end{aligned} \tag{1}\]

On the hole edge and along the net section \(x = 0\) this reduces to

\[ \sigma_{\theta\theta}(a, \theta) = \sigma_\infty\,(1 - 2\cos 2\theta), \qquad \sigma_{xx}(0, y) = \frac{\sigma_\infty}{2}\left(2 + \frac{a^2}{y^2} + \frac{3a^4}{y^4}\right), \tag{2}\]

so \(K_t = \sigma_{\max}/\sigma_\infty = 3\) at \(\theta = 90°\), whatever the elastic constants. A finite width \(W\) raises it. Peterson’s curve fit for a central hole, referred to the gross stress, is

\[ K_t = \frac{2 + 0.284\,\beta - 0.600\,\beta^2 + 1.32\,\beta^3}{\beta}, \qquad \beta = 1 - \frac{d}{W}, \tag{3}\]

which gives \(K_t\) = 3.73 for \(d/W\) = 0.4 (Pilkey & Pilkey, Peterson’s Stress Concentration Factors, 3rd ed., chart 4.1).

FEniCSx implementation

The plate is a standard linear-elasticity problem: find \(\underline{u} \in V\), with \(\underline{u} \cdot \underline{n} = 0\) on the symmetry edges, such that

\[ \int_\Omega \underline{\underline{\sigma}}(\underline{u}) : \underline{\underline{\varepsilon}}(\underline{v})\,\mathrm d\Omega = \int_{\Gamma_\text{load}} \sigma_\infty\, v_x \,\mathrm d\Gamma \quad \forall\, \underline{v} \in V_0, \qquad \underline{\underline{\sigma}} = \lambda^*\operatorname{tr}(\underline{\underline{\varepsilon}})\,\underline{\underline{I}} + 2\mu\,\underline{\underline{\varepsilon}}, \quad \lambda^* = \frac{E\nu}{1 - \nu^2}. \tag{4}\]

Condensed excerpts from the 2D scripts (... marks omitted lines):

2D/mesh.py

rectangle = gmsh.model.occ.addRectangle(0, 0, 0, half_length, half_width)
hole = gmsh.model.occ.addDisk(0, 0, 0, hole_radius, hole_radius)
1plate, _ = gmsh.model.occ.cut([(2, rectangle)], [(2, hole)])
...
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "CurvesList", hole_curves)
2threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "InField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "SizeMin", hole_element_length)
gmsh.model.mesh.field.setNumber(threshold_field, "SizeMax", element_length)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", refinement_distance)
gmsh.model.mesh.field.setAsBackgroundMesh(threshold_field)
...
gmsh.option.setNumber("Mesh.Algorithm", 8)  # Frontal-Delaunay for quads
3gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.model.mesh.generate(2)

4mesh_data = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

1: Quarter plate with the hole centred at the origin, so \(x = 0\) and \(y = 0\) are symmetry lines.
2: Element size grows from 5 mm at the hole to 10 mm over 5 hole radii.
3: Triangles are recombined into quadrilaterals. In 3D the section is extruded in layers into hexahedra.
4: Physical groups tag the symmetry and loaded edges; they arrive in DOLFINx as facet tags.

2D/simulation.py

1displacement_space = fem.functionspace(domain, ("P", element_degree, (geom_dim,)))
ds_load = ufl.Measure("ds", domain=domain, subdomain_data=tags)

2def symmetry_condition(facet_tag, normal_component):
    component_space = displacement_space.sub(normal_component)
    symmetry_dofs = fem.locate_dofs_topological(component_space, bc_dim, tags.find(facet_tag))
    return fem.dirichletbc(np.float64(0.0), symmetry_dofs, component_space)

boundary_conditions = [
    symmetry_condition(symmetry_x_tag, 0),
    symmetry_condition(symmetry_y_tag, 1),
]

shear_modulus = young_modulus/(2*(1+poisson_ratio))
3lame_plane_stress = young_modulus*poisson_ratio/(1-poisson_ratio**2)

def strain(displacement):
    return ufl.sym(ufl.grad(displacement))

def stress(strain_tensor):
    return lame_plane_stress*ufl.tr(strain_tensor)*ufl.Identity(geom_dim) + 2*shear_modulus*strain_tensor

bilinear_form = ufl.inner(stress(strain(trial)), strain(test))*ufl.dx
4linear_form = remote_stress*test[0]*ds_load(loaded_end_tag)

1: Vector Lagrange space of degree 1.
2: On a symmetry edge only the normal displacement component is fixed.
3: Plane stress enters through \(\lambda^*\). The 3D model uses \(\lambda = E\nu/((1+\nu)(1-2\nu))\) and adds the mid-plane \(z = 0\) as a third symmetry plane.
4: The remote stress is a traction on the loaded end, integrated over its tagged facets.

2D/simulation.py

problem = LinearProblem(
    bilinear_form, linear_form, u=solution, bcs=boundary_conditions,
1    petsc_options=DIRECT_SOLVER,
    petsc_options_prefix="stress_concentration_2d_",
)
problem.solve()

# The stress from the displacement gradient jumps between elements; its L2
# projection onto a continuous space is smooth and more accurate at the hole.
2projection_bilinear_form = ufl.inner(stress_trial, stress_test)*ufl.dx
projection_linear_form = ufl.inner(stress(strain(solution)), stress_test)*ufl.dx

projection_problem = LinearProblem(
    projection_bilinear_form, projection_linear_form, u=stress_function,
3    petsc_options=PROJECTION_SOLVER,
    petsc_options_prefix="stress_concentration_2d_stress_projection_",
)
projection_problem.solve()

1: DIRECT_SOLVER: LU factorisation with MUMPS through PETSc.
2: \(L^2\) projection onto continuous tensors: find \(\underline{\underline{\sigma}}{}_h\) with \(\int \underline{\underline{\sigma}}{}_h : \underline{\underline{\tau}} = \int \underline{\underline{\sigma}}(\underline{u}{}_h) : \underline{\underline{\tau}}\) for all \(\underline{\underline{\tau}}\).
3: PROJECTION_SOLVER: conjugate gradients with a Jacobi preconditioner, since the mass matrix is well conditioned.

Results

The FE stresses are sampled on the hole edge and along the net section, and compared with Equation 2 and Equation 3. The 3D model is sampled at its mid-plane and on its free surface.

Figure 3: Hoop stress around the hole edge, from the load axis (\(\theta = 0°\)) to the net section (\(\theta = 90°\)).

Figure 4: Normal stress across the net section \(x = 0\), from the hole edge (\(y = a\)) to the plate edge (\(y = W/2\)).

Finite width matters. With \(d/W\) = 0.4 the FE peak is +25.9% above Kirsch’s \(K_t = 3\). The curves keep Kirsch’s shape but with a larger amplitude: the compressed flank at \(\theta = 0\) reaches −1.46 \(\sigma_\infty\), against \(-\sigma_\infty\) for Kirsch.
Peterson’s chart captures it. The FE factors are +0.7% to +2.0% from Peterson’s 3.73.
The 3D peak is at the mid-plane. Through-thickness constraint is strongest inside the plate, so \(K_t\) drops from 3.81 at the mid-plane to 3.76 on the free surface, bracketing the plane-stress value 3.78.

Table 1: Stress concentration factors, referred to the gross remote stress.

\(K_t = \sigma_{xx}(0, a)/\sigma_\infty\)	Value	vs Peterson	vs Kirsch
FE · 2D plane stress	3.777	+1.2%	+25.9%
FE · 3D mid-plane	3.809	+2.0%	+27.0%
FE · 3D free surface	3.760	+0.7%	+25.3%
Peterson (finite width)	3.732
Kirsch (infinite plate)	3.000