Skew the Element 45° and the Stiffness Grows 2.5x — Metric Tensors and Constitutive Tensor Transformation

Measure the same deformation twice and the answers split. The element did not change. Neither did the material. Neither did the physical strain. The only swap is the coordinate system the strain was written in. The strain energy still comes out 2.5 times larger. The factor sits in the covariant and contravariant bases and in the metric, and the Python check shows which rotation of the fourth-order tensor pins the energy back at 1.412308.

Drag skew from 0 to 45° or raise |g_2| above 1 and the white arrow stays put while the covariant and contravariant numbers split.

At skew 0 and |g_2| 1.00 the blue parallelogram and the yellow perpendiculars meet, and max | v^i - v_i | turns green. A few degrees of skew pull the two rows apart. The arrow itself does not move. That split is where a skewed element starts reporting the wrong stiffness.

The deformation held still and the energy grew 2.5x#

A shell stiffness matrix is usually assembled like this. At each Gauss point you build the strain-displacement matrix $\mathbf{B}$, multiply by the constitutive matrix $\mathbf{C}$, and integrate $\mathbf{B}^{\mathsf T}\mathbf{C}\mathbf{B}$. The trouble is that these two matrices are born in different coordinate systems.

$\mathbf{C}$ is a material property. It is defined in the frame where the material tests were run: the local orthogonal Cartesian system glued to the shell surface. $\mathbf{B}$, by contrast, comes from differentiating shape functions. Shape functions are written in the element's natural coordinate system $(r^1, r^2, r^3)$, and once the element bends or distorts, that system is neither orthogonal nor normalized.

There is one case where the two frames happen to coincide: a rectangular element on a flat surface. Run your patch test on a planar model only and you validate nothing but that coincidence. Wrap the element onto a curved surface and the coincidence is gone.

A non-orthogonal basis gives every vector two sets of components#

The direction the physical position $\mathbf{x}$ moves when the natural coordinate $r^i$ changes — that is the covariant basis.

Here $\mathbf{x}$ is the Cartesian position vector and $r^i$ are the element's natural coordinates. These three vectors are not mutually orthogonal, and their lengths are not 1.

On a non-orthogonal basis there are two ways to write one vector in components. Decompose it along the base vectors with the parallelogram rule and the coefficients are the contravariant components $v^i$. Drop it perpendicularly onto each base vector and the projections are the covariant components $v_i$.

On an orthonormal basis both constructions land on the same point. A Cartesian-only eye never sees the split. The dial at the top of the note already pulled the two rows apart. The white arrow is the same vector after the coordinates change.

The metric tensor gives back the length that was lost#

What links the two sets of components is the metric tensor.

$g_{ij}$ packs the lengths of the base vectors and the angles between them into one matrix. The diagonal entries $g_{11}, g_{22}, g_{33}$ are the squared lengths; the off-diagonal entries are the cosine of the included angle times those lengths.

To measure a length you must go through this matrix.

The second equality is the useful one. Pair a covariant component with a contravariant one and the metric cancels itself out. Square the contravariant components and add them up ($v^i v^i$) and you do not get a length. The bottom line of the visualization shows that value in red.

This is exactly why continuum mechanics uses contravariant components for stress and covariant components for strain. Virtual work has to be a scalar, and only the pairing of the second Piola-Kirchhoff contravariant components $S^{ij}$ with the Green-Lagrange covariant components $E_{ij}$ makes $S^{ij}E_{ij}$ independent of the frame. Grid metrics do the same job in flow solvers — the Jacobian covered in curvilinear coordinate transformation and grid metrics is, here, literally the collection of base vectors.

The contravariant basis is the rows of the inverse Jacobian#

To convert covariant components back to contravariant ones you need $g^{ij} = [g_{ij}]^{-1}$, and the basis built from that matrix is the contravariant basis.

$\delta^i_j$ is the Kronecker delta. So $\mathbf{g}^1$ is the vector perpendicular to both $\mathbf{g}_2$ and $\mathbf{g}_3$, scaled so that its dot product with $\mathbf{g}_1$ is exactly 1.

In code you never need two matrix inversions. Stack the covariant base vectors as the columns of a matrix $\mathbf{J}$:

$\big(\mathbf{J}^{-1}\big)_{i\bullet}$ denotes the $i$-th row of $\mathbf{J}^{-1}$.

$\mathbf{J}$ is the Jacobian from natural to physical coordinates. It is the same matrix you already compute at every integration point to get $\det\mathbf{J}$. The contravariant basis comes for free.

A fourth-order tensor takes four direction cosines#

Now the main event. Move the constitutive tensor $C^{pqrs}$, defined on the local Cartesian basis $\mathbf{e}_p$, into the natural coordinate system. A second-order tensor picks up two direction cosines; a fourth-order tensor picks up four.

Indices $i,j,k,l$ live in natural coordinates, $p,q,r,s$ in local Cartesian ones. The strain side runs the other way.

The constitutive tensor picks up $\mathbf{J}^{-1}$, the strain picks up $\mathbf{J}$. Contract them and the Jacobians cancel exactly, leaving the energy invariant. Turn that around: transform only the strain and leave the constitutive tensor alone, and four copies of $\mathbf{J}$ survive. Those four copies are the error column in the table. How that column swells with skew and |g_2| shows up as a bar right after the printout.

Energy per skew angle, measured in Python#

I fixed one strain state on an isotropic plane-stress material, then twisted only the coordinate system and computed the strain energy two ways. The material constants match the source notes ($E = 2.1\times10^6$, $\nu = 0.3$).

import numpy as np
 
def natural_basis(skew_deg, stretch=1.0):
    """Jacobian J whose columns are the covariant base vectors g_i = dx/dr^i"""
    a = np.deg2rad(skew_deg)
    g1 = np.array([1.0, 0.0])
    g2 = stretch * np.array([np.sin(a), np.cos(a)])
    return np.column_stack([g1, g2])
 
def plane_stress_tensor(E=2.1e6, nu=0.3):
    """Isotropic plane-stress fourth-order tensor C^{pqrs} in local Cartesian coordinates"""
    lam = E * nu / (1.0 - nu**2)
    mu = E / (2.0 * (1.0 + nu))
    d = np.eye(2)
    return (lam * np.einsum('pq,rs->pqrs', d, d)
            + mu * (np.einsum('pr,qs->pqrs', d, d) + np.einsum('ps,qr->pqrs', d, d)))
 
def rotate_fourth_order(C, Jinv):
    """C^{ijkl} = (g^i.e_p)(g^j.e_q)(g^k.e_r)(g^l.e_s) C^{pqrs}, g^i = i-th row of J^-1"""
    return np.einsum('ip,jq,kr,ls,pqrs->ijkl', Jinv, Jinv, Jinv, Jinv, C)
 
def strain_energy(C, eps):
    return 0.5 * np.einsum('pqrs,pq,rs->', C, eps, eps)
 
def to_voigt2d(C):
    """2D Voigt: (00,11,01) -> 3x3"""
    idx = [(0, 0), (1, 1), (0, 1)]
    return np.array([[C[p, q, r, s] for (r, s) in idx] for (p, q) in idx])
 
# One physical strain state (local Cartesian components). This tensor never changes,
# no matter how the coordinates are chosen.
eps_cart = np.array([[1.0e-3, 4.0e-4],
                     [4.0e-4, -6.0e-4]])
C_cart = plane_stress_tensor()
U_ref = strain_energy(C_cart, eps_cart)
 
print("skew   g11     g12     g22    | U_correct     U_naive      err%")
print("-" * 68)
for skew in [0, 5, 10, 15, 20, 30, 40, 45]:
    J = natural_basis(skew)
    g = J.T @ J                       # metric tensor g_ij
    Jinv = np.linalg.inv(J)           # rows = contravariant basis g^i
    eps_nat = J.T @ eps_cart @ J      # covariant strain components
    C_nat = rotate_fourth_order(C_cart, Jinv)
    U_ok = strain_energy(C_nat, eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)   # code that forgot the transform
    err = 100.0 * (U_bad - U_ok) / U_ok
    print(f"{skew:3d}   {g[0,0]:.3f}  {g[0,1]:+.3f}  {g[1,1]:.3f} |"
          f" {U_ok:.6e}  {U_bad:.6e}  {err:+8.2f}")
 
print()
print("orthogonal (skew=0), only |g2| stretched")
for st in [1.0, 1.5, 2.0]:
    J = natural_basis(0, stretch=st)
    Jinv = np.linalg.inv(J)
    eps_nat = J.T @ eps_cart @ J
    U_ok = strain_energy(rotate_fourth_order(C_cart, Jinv), eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)
    print(f"  stretch={st:.1f}  g22={(J.T@J)[1,1]:.2f}  err% = {100*(U_bad-U_ok)/U_ok:+9.2f}")
 
print()
print(f"Cartesian reference U_ref   = {U_ref:.6e}")
J = natural_basis(30)
Jinv = np.linalg.inv(J)
C_nat = rotate_fourth_order(C_cart, Jinv)
eps_nat = J.T @ eps_cart @ J
print(f"skew=30 after transform     = {strain_energy(C_nat, eps_nat):.6e}  (invariant)")
 
# Does folding into Voigt form give the same number?
Cv = to_voigt2d(C_nat)
ev = np.array([eps_nat[0, 0], eps_nat[1, 1], 2.0 * eps_nat[0, 1]])
print(f"skew=30 via Voigt 3x3       = {0.5 * ev @ Cv @ ev:.6e}")
 
# Strain transformation matrix A in Voigt space: e_v(nat) = A e_v(cart)
def voigt_map(J):
    cols = []
    for e in (np.array([[1, 0], [0, 0]]), np.array([[0, 0], [0, 1]]), np.array([[0, .5], [.5, 0]])):
        n = J.T @ e @ J
        cols.append([n[0, 0], n[1, 1], 2 * n[0, 1]])
    return np.array(cols).T
 
A = voigt_map(J)
Cv_cart = to_voigt2d(C_cart)
Ai = np.linalg.inv(A)
print("Voigt congruence C_nat = A^-T C_cart A^-1  residual =",
      f"{np.max(np.abs(Ai.T @ Cv_cart @ Ai - Cv)):.3e}")

skew   g11     g12     g22    | U_correct     U_naive      err%
--------------------------------------------------------------------
  0   1.000  +0.000  1.000 | 1.412308e+00  1.412308e+00     +0.00
  5   1.000  +0.087  1.000 | 1.412308e+00  1.486003e+00     +5.22
 10   1.000  +0.174  1.000 | 1.412308e+00  1.585621e+00    +12.27
 15   1.000  +0.259  1.000 | 1.412308e+00  1.722495e+00    +21.96
 20   1.000  +0.342  1.000 | 1.412308e+00  1.906550e+00    +35.00
 30   1.000  +0.500  1.000 | 1.412308e+00  2.437219e+00    +72.57
 40   1.000  +0.643  1.000 | 1.412308e+00  3.163176e+00   +123.97
 45   1.000  +0.707  1.000 | 1.412308e+00  3.567692e+00   +152.61
 
orthogonal (skew=0), only |g2| stretched
  stretch=1.0  g22=1.00  err% =     +0.00
  stretch=1.5  g22=2.25  err% =   +105.60
  stretch=2.0  g22=4.00  err% =   +407.84
 
Cartesian reference U_ref   = 1.412308e+00
skew=30 after transform     = 1.412308e+00  (invariant)
skew=30 via Voigt 3x3       = 1.412308e+00
Voigt congruence C_nat = A^-T C_cart A^-1  residual = 9.313e-10

Push skew from 0 to 45° or raise |g_2| above 1 and the element on the left keeps its shape while only the red bar on the right grows. Switch to shear and that bar climbs fastest.

The +152.61% at 45° and the +105.60% at stretch 1.5 are those two bars. stretch lifts the error curve on |g_2| alone, with the angle still at 0.

Every row of the table is the same subtraction. The correct energy stays pinned at 1.412308. Twist the coordinates out to 45° or fold the tensor into Voigt and that number does not move. The untouched side grows. At 5° it is 1.486003, error +5.22%. At 15° it is 1.722495, +21.96%. Divide 1.722495 by 1.412308 and the quotient is about 1.22, which matches the percent in the table. Fifteen degrees is an ordinary angle on a curved mesh. A suite that only runs rectangles sees the 0.00 on the skew=0 row and walks past.

The 45° row is the 2.5 in the title. Divide the wrong energy 3.567692 by 1.412308 and the quotient is about 2.53. The gap alone is 2.155. There $g_{12}$ is +0.707, which is $\sin 45^\circ$. At 30°, $g_{12} = +0.500 = \sin 30^\circ$, the error is +72.57%, and the energy is 2.437219. Zero angle still blows the same column up if the length is wrong. At $|\mathbf{g}_2| = 1.5$, $g_{22} = 2.25$ and the error is +105.60%. At 2.0, $g_{22} = 4.00$ and the error is +407.84%. The cause is not the skewed angle. It is the metric not being the identity. The residual $9.3\times 10^{-10}$ on the last line only says the fourth-order contraction and the Voigt congruence both return 1.412308. It does not rescue the column that skipped the transform.

Folded into Voigt form it becomes a single 6x6 matrix#

Nobody carries a four-index array around in production code. Stress and strain tensors are symmetric, so only six independent components survive and the fourth-order tensor folds into a $6\times6$ matrix (the code above is 2D, so $3\times3$).

The transformation survives the folding. If the Voigt strain vector maps as $\tilde{\boldsymbol\varepsilon}_v = \mathbf{A}\,\boldsymbol\varepsilon_v$, then energy invariance forces a congruence transformation on the constitutive matrix.

The entries of $\mathbf{A}$ are products of direction cosines. The residual $9.3\times10^{-10}$ on the last line of the output confirms that this route agrees with the fourth-order contraction. That is the T matrix you see in real shell codes, and the shear correction factor of 5/6 plus the plane-stress assumption ($\sigma_{33}=0$) go into $\mathbf{C}_v$ first, in the local Cartesian frame. Do not reverse the order — the plane-stress condition only means something in the frame where the thickness direction is defined.

Where the same mistake shows up in finite volume solvers#

This is not a structural-code-only mistake. The same structure shows up when a curvilinear finite volume solver builds the viscous stress. Take the strain rate from derivatives in natural coordinates and multiply Newton's law $\tau = 2\mu \mathbf{D}$ in Cartesian form, and the error column above comes back whole. That is the +152.61% at 45°. Whether the components in hand are physical or covariant, and whether every upper index meets a lower one in the contraction, comes first. A suite that only passes on orthogonal grids never sees that column.

Cell area from the diagonal cross product hands Green–Gauss the same metric. If shear locking and MITC tying is about fixing the element's formulation, the question here is which coordinate system that formulation is read in. Fix the formulation or miss the coordinates, and the planar patch test still passes. One element at 45°, or one with $|\mathbf{g}_2| = 1.5$, is what splits the table.