import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench import (
build_elasticity_triangle_p1_2d, AffineTriangleMap2D,
ReferenceTriangleP1,
)
from symbolic_fem_workbench.elasticity import (
plane_stress_D, plane_strain_D, B_matrix_triangle_2d,
element_stiffness_BtDB,
)11 2D Linear Elasticity with P1 Triangles
This notebook extends the scalar Poisson problem to vector elasticity. Instead of a single scalar field \(u\), we now solve for two displacement components \((u_x, u_y)\) at each node. The key new concepts are:
- Voigt notation: Compact vector representation of stresses and strains
- Constitutive relations: How stresses relate to strains (linear elastic material)
- Strain-displacement matrix (B): Maps nodal displacements to element strains
- Plane stress vs. plane strain: Different constitutive models for 2D elasticity
We derive the element stiffness matrix \(K_e = t A B^T D B\) where \(D\) is the constitutive matrix and \(B\) is the strain-displacement matrix.
11.1 Imports
11.2 Voigt Notation: From Tensors to Vectors
In 2D elasticity, we work with the symmetric stress and strain tensors:
\[\boldsymbol{\sigma} = \begin{bmatrix} \sigma_{xx} & \tau_{xy} \\ \tau_{xy} & \sigma_{yy} \end{bmatrix}, \quad \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} & \gamma_{xy}/2 \\ \gamma_{xy}/2 & \varepsilon_{yy} \end{bmatrix}\]
Voigt notation compresses these into 3-component vectors (for 2D problems):
\[\mathbf{\sigma} = \begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \tau_{xy} \end{bmatrix}, \quad \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \gamma_{xy} \end{bmatrix}\]
Note: The shear strain \(\gamma_{xy} = 2 \varepsilon_{xy}\) (engineering shear strain).
The advantage: we can write the constitutive relation as a matrix equation: \[\mathbf{\sigma} = \mathbf{D} \boldsymbol{\varepsilon}\]
where \(\mathbf{D}\) is the 3×3 constitutive matrix (also called the elasticity or stiffness matrix).
11.3 Plane Stress vs. Plane Strain
For 2D problems, we make a simplifying assumption about the out-of-plane stress/strain.
E, nu = sp.symbols("E nu", positive=True)
# Plane stress: sigma_zz = 0 (thin plates, no out-of-plane loading)
D_stress = plane_stress_D(E, nu)
print("Plane stress D (thin plates):")
print("\nD =")
display(Math(sp.latex(D_stress)))
# Plane strain: epsilon_zz = 0 (thick bodies, constrained out-of-plane)
D_strain = plane_strain_D(E, nu)
print("\n\nPlane strain D (thick bodies, constrained):")
print("\nD =")
display(Math(sp.latex(D_strain)))Plane stress D (thin plates):
D =\(\displaystyle \left[\begin{matrix}\frac{E}{1 - \nu^{2}} & \frac{E \nu}{1 - \nu^{2}} & 0\\\frac{E \nu}{1 - \nu^{2}} & \frac{E}{1 - \nu^{2}} & 0\\0 & 0 & \frac{E \left(\frac{1}{2} - \frac{\nu}{2}\right)}{1 - \nu^{2}}\end{matrix}\right]\)
Plane strain D (thick bodies, constrained):
D =\(\displaystyle \left[\begin{matrix}\frac{E \left(1 - \nu\right)}{\left(1 - 2 \nu\right) \left(\nu + 1\right)} & \frac{E \nu}{\left(1 - 2 \nu\right) \left(\nu + 1\right)} & 0\\\frac{E \nu}{\left(1 - 2 \nu\right) \left(\nu + 1\right)} & \frac{E \left(1 - \nu\right)}{\left(1 - 2 \nu\right) \left(\nu + 1\right)} & 0\\0 & 0 & \frac{E \left(\frac{1}{2} - \nu\right)}{\left(1 - 2 \nu\right) \left(\nu + 1\right)}\end{matrix}\right]\)
11.4 Strain-Displacement Matrix (B): From Displacements to Strains
For a P1 triangle with 3 nodes and 2 DOF per node (total 6 DOF), the strain-displacement relation is:
\[\boldsymbol{\varepsilon} = \mathbf{B} \mathbf{d}\]
where \(\mathbf{d} = [u_1, v_1, u_2, v_2, u_3, v_3]^T\) are the nodal displacements and \(\mathbf{B}\) is 3×6.
The B matrix is built from the shape function gradients \(\frac{\partial N_i}{\partial x}\) and \(\frac{\partial N_i}{\partial y}\):
\[\mathbf{B} = \begin{bmatrix} \frac{\partial N_1}{\partial x} & 0 & \frac{\partial N_2}{\partial x} & 0 & \frac{\partial N_3}{\partial x} & 0 \\ 0 & \frac{\partial N_1}{\partial y} & 0 & \frac{\partial N_2}{\partial y} & 0 & \frac{\partial N_3}{\partial y} \\ \frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial y} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial y} & \frac{\partial N_3}{\partial x} \end{bmatrix}\]
# Build B for the unit right triangle
xi, eta = sp.symbols("xi eta")
ref = ReferenceTriangleP1(xi, eta)
# Unit right triangle vertices: (0,0), (1,0), (0,1)
x1, y1 = 0, 0
x2, y2 = 1, 0
x3, y3 = 0, 1
geom = AffineTriangleMap2D(xi=xi, eta=eta, x1=x1, y1=y1, x2=x2, y2=y2, x3=x3, y3=y3)
# Physical-space gradients via J^{-T}
JinvT = geom.inverse_transpose_jacobian
dN_dx = []
dN_dy = []
for grad_ref in ref.shape_gradients_reference:
grad_phys = sp.simplify(JinvT * grad_ref)
dN_dx.append(grad_phys[0])
dN_dy.append(grad_phys[1])
B = B_matrix_triangle_2d(dN_dx, dN_dy)
print("B matrix for unit right triangle (3×6):")
print("(Each row represents one strain component: ε_xx, ε_yy, γ_xy)")
print("(Each pair of columns represents one node's DOF: u_i, v_i)\n")
display(Math(r"B = " + sp.latex(B)))B matrix for unit right triangle (3×6):
(Each row represents one strain component: ε_xx, ε_yy, γ_xy)
(Each pair of columns represents one node's DOF: u_i, v_i)
\(\displaystyle B = \left[\begin{matrix}-1 & 0 & 1 & 0 & 0 & 0\\0 & -1 & 0 & 0 & 0 & 1\\-1 & -1 & 0 & 1 & 1 & 0\end{matrix}\right]\)
11.5 Element Stiffness Matrix: K_e = t A B^T D B
For a constant-strain element (P1 triangle), the strain and D are constant over the element. The weak form integral simplifies to:
\[K_e = t \cdot A \cdot B^T D B\]
where: - \(t\) is the plate thickness (for 2D plane stress problems) - \(A\) is the element area - \(B\) is the strain-displacement matrix (3×6) - \(D\) is the constitutive matrix (3×3) - \(K_e\) is the 6×6 element stiffness matrix
# Use the pre-built workflow to get Ke
data = build_elasticity_triangle_p1_2d("plane_stress")
Ke = data["Ke_unit_right_triangle"]
print("6×6 element stiffness matrix for unit right triangle (plane stress):")
print("(with E=1, ν=1/3, t=1 for readability)\n")
# Substitute specific values for display
Ke_specific = Ke.subs({data["E"]: 1, data["nu"]: sp.Rational(1, 3), data["t"]: 1})
display(Math(r"K_e = " + sp.latex(sp.simplify(Ke_specific))))
print("\nProperties:")
print(f" Symmetric: {sp.simplify(Ke - Ke.T) == sp.zeros(6, 6)}")
print(f" Shape: {Ke.shape}")6×6 element stiffness matrix for unit right triangle (plane stress):
(with E=1, ν=1/3, t=1 for readability)
\(\displaystyle K_e = \left[\begin{matrix}\frac{3}{4} & \frac{3}{8} & - \frac{9}{16} & - \frac{3}{16} & - \frac{3}{16} & - \frac{3}{16}\\\frac{3}{8} & \frac{3}{4} & - \frac{3}{16} & - \frac{3}{16} & - \frac{3}{16} & - \frac{9}{16}\\- \frac{9}{16} & - \frac{3}{16} & \frac{9}{16} & 0 & 0 & \frac{3}{16}\\- \frac{3}{16} & - \frac{3}{16} & 0 & \frac{3}{16} & \frac{3}{16} & 0\\- \frac{3}{16} & - \frac{3}{16} & 0 & \frac{3}{16} & \frac{3}{16} & 0\\- \frac{3}{16} & - \frac{9}{16} & \frac{3}{16} & 0 & 0 & \frac{9}{16}\end{matrix}\right]\)
Properties:
Symmetric: True
Shape: (6, 6)11.6 Plane Stress vs. Plane Strain Comparison
Let’s compare the two formulations on the same unit right triangle:
# Build for both formulations
data_ps = build_elasticity_triangle_p1_2d("plane_stress")
data_pe = build_elasticity_triangle_p1_2d("plane_strain")
Ke_ps = data_ps["Ke_unit_right_triangle"]
Ke_pe = data_pe["Ke_unit_right_triangle"]
# Compare with E=1, nu=0.3, t=1
E_val, nu_val, t_val = 1, sp.Rational(3, 10), 1
Ke_ps_num = Ke_ps.subs({data_ps["E"]: E_val, data_ps["nu"]: nu_val, data_ps["t"]: t_val})
Ke_pe_num = Ke_pe.subs({data_pe["E"]: E_val, data_pe["nu"]: nu_val, data_pe["t"]: t_val})
print("Plane stress K_e (E=1, ν=0.3, t=1):")
display(Math(sp.latex(Ke_ps_num)))
print("\nPlane strain K_e (E=1, ν=0.3, t=1):")
display(Math(sp.latex(Ke_pe_num)))
print(f"\nStiffness ratio (plane strain / plane stress):")
ratio = sp.simplify((1 - nu_val) / (1 - 2*nu_val))
print(f" Factor ≈ {float(ratio):.4f}")
print(f" Plane strain is stiffer due to out-of-plane constraint.")Plane stress K_e (E=1, ν=0.3, t=1):\(\displaystyle \left[\begin{matrix}\frac{135}{182} & \frac{5}{14} & - \frac{50}{91} & - \frac{5}{26} & - \frac{5}{26} & - \frac{15}{91}\\\frac{5}{14} & \frac{135}{182} & - \frac{15}{91} & - \frac{5}{26} & - \frac{5}{26} & - \frac{50}{91}\\- \frac{50}{91} & - \frac{15}{91} & \frac{50}{91} & 0 & 0 & \frac{15}{91}\\- \frac{5}{26} & - \frac{5}{26} & 0 & \frac{5}{26} & \frac{5}{26} & 0\\- \frac{5}{26} & - \frac{5}{26} & 0 & \frac{5}{26} & \frac{5}{26} & 0\\- \frac{15}{91} & - \frac{50}{91} & \frac{15}{91} & 0 & 0 & \frac{50}{91}\end{matrix}\right]\)
Plane strain K_e (E=1, ν=0.3, t=1):\(\displaystyle \left[\begin{matrix}\frac{45}{52} & \frac{25}{52} & - \frac{35}{52} & - \frac{5}{26} & - \frac{5}{26} & - \frac{15}{52}\\\frac{25}{52} & \frac{45}{52} & - \frac{15}{52} & - \frac{5}{26} & - \frac{5}{26} & - \frac{35}{52}\\- \frac{35}{52} & - \frac{15}{52} & \frac{35}{52} & 0 & 0 & \frac{15}{52}\\- \frac{5}{26} & - \frac{5}{26} & 0 & \frac{5}{26} & \frac{5}{26} & 0\\- \frac{5}{26} & - \frac{5}{26} & 0 & \frac{5}{26} & \frac{5}{26} & 0\\- \frac{15}{52} & - \frac{35}{52} & \frac{15}{52} & 0 & 0 & \frac{35}{52}\end{matrix}\right]\)
Stiffness ratio (plane strain / plane stress):
Factor ≈ 1.7500
Plane strain is stiffer due to out-of-plane constraint.11.7 Material Parameters and Their Effect
The element stiffness depends on: - E (Young’s modulus): Higher E → stiffer material - ν (Poisson’s ratio): Controls how transverse strains relate to axial strains - t (thickness): For plane stress, stiffer with larger t
# Show how the D matrix changes with Poisson's ratio
D_template = plane_stress_D(1, nu) # E=1 for simplicity
print("Plane stress D for E=1:")
display(Math(sp.latex(D_template)))
print(f"\nKey observations:")
print(f" - Diagonal terms scale with 1/(1-ν²), increasing with ν")
print(f" - Off-diagonal coupling (D[0,1], D[1,0]) = ν·E/(1-ν²)")
print(f" - Shear modulus D[2,2] = E/(2(1+ν)) decreases with ν")
# Numerical examples
print(f"\nFor E=1, plane stress:")
for nu_test in [0.0, 0.2, 0.3, 0.49]:
D_val = D_template.subs(nu, nu_test)
diag_00 = float(D_val[0, 0])
print(f" ν = {nu_test:5.2f}: D[0,0] = {diag_00:8.4f}")Plane stress D for E=1:\(\displaystyle \left[\begin{matrix}\frac{1}{1 - \nu^{2}} & \frac{\nu}{1 - \nu^{2}} & 0\\\frac{\nu}{1 - \nu^{2}} & \frac{1}{1 - \nu^{2}} & 0\\0 & 0 & \frac{\frac{1}{2} - \frac{\nu}{2}}{1 - \nu^{2}}\end{matrix}\right]\)
Key observations:
- Diagonal terms scale with 1/(1-ν²), increasing with ν
- Off-diagonal coupling (D[0,1], D[1,0]) = ν·E/(1-ν²)
- Shear modulus D[2,2] = E/(2(1+ν)) decreases with ν
For E=1, plane stress:
ν = 0.00: D[0,0] = 1.0000
ν = 0.20: D[0,0] = 1.0417
ν = 0.30: D[0,0] = 1.0989
ν = 0.49: D[0,0] = 1.316011.8 Worked Example: Simple 2-Element Cantilever
Consider a small cantilever beam modeled with 2 right triangles. This demonstrates assembly.
# Two-triangle cantilever:
# Element 1: (0,0) - (1,0) - (0,1)
# Element 2: (1,0) - (1,1) - (0,1)
print("Example mesh: Two right triangles forming a square")
print("Nodes: (0,0), (1,0), (0,1), (1,1)")
print("Elem 1: nodes 0-1-2")
print("Elem 2: nodes 1-3-2")
print()
print("With plane stress: E=1000, ν=0.3, t=0.1")
# Get element stiffness for one element
data = build_elasticity_triangle_p1_2d("plane_stress")
Ke_1 = data["Ke_unit_right_triangle"]
Ke_1_num = Ke_1.subs({data["E"]: 1000, data["nu"]: 0.3, data["t"]: 0.1})
# Convert to float for better display
Ke_1_float = sp.Matrix([[float(x) for x in row] for row in Ke_1_num.tolist()])
print("\nElement stiffness Ke (unit right triangle):")
print(Ke_1_float.evalf(3))Example mesh: Two right triangles forming a square
Nodes: (0,0), (1,0), (0,1), (1,1)
Elem 1: nodes 0-1-2
Elem 2: nodes 1-3-2
With plane stress: E=1000, ν=0.3, t=0.1
Element stiffness Ke (unit right triangle):
Matrix([[74.2, 35.7, -54.9, -19.2, -19.2, -16.5], [35.7, 74.2, -16.5, -19.2, -19.2, -54.9], [-54.9, -16.5, 54.9, 0, 0, 16.5], [-19.2, -19.2, 0, 19.2, 19.2, 0], [-19.2, -19.2, 0, 19.2, 19.2, 0], [-16.5, -54.9, 16.5, 0, 0, 54.9]])11.9 Exercise
Verify symmetry: Show that \(K_e = K_e^T\) for both plane stress and plane strain formulations. Why is this property important?
Effect of aspect ratio: Compute the B matrix and \(K_e\) for a tall, thin triangle with vertices at \((0,0)\), \((0.1, 0)\), \((0, 1)\). How do the stiffness values change compared to the unit right triangle?
3D extension: The module also provides
isotropic_3d_D()andB_matrix_tetra_3d(). Inspect these functions and explain how the Voigt notation extends to 6 components for 3D.Poisson effect: Fix E and t, and plot the diagonal entries of \(K_e\) as functions of \(\nu\) for \(\nu \in [0, 0.49]\). Which entries increase or decrease with \(\nu\)?
# Exercise workspace