QF Solver
Get started
White box / Physical formulations / Shell elements

MITC4 · four-node quadrilateral shell

Physical formulation · assumptions, essential equations, full derivation, algorithm, code and verification references.

← All formulations and methodsPhysical formulation
Browse another sheet Physical formulations / Shell elements
01 Family
02 Chapter
03 Sheet
Physical formulationsShell elementsSheet 1 / 2
Physical formulationROUTE-DEPENDENT · separate modal, transient, and harmonic evidence

Bilinear Q4 midsurface interpolation, six nodal DOFs, membrane, bending, transverse shear, and a weak drilling stabilizer.

01 Model and scope

Reissner–Mindlin kinematics retain transverse shear by allowing rotations independent of the midsurface slope. MITC4 samples covariant shear at tying points and interpolates it across the element; membrane, bending, shear and drilling terms remain distinct contributions.

Route assumptions and calculation checks

  1. MITC tying points interpolate covariant transverse shear to reduce shear locking.
  2. The drilling term is a weak penalty and requires sensitivity checks.
  3. Maturity remains specific to analysis, mass, damping, laminate, and benchmark route.

Useful foundations and prerequisites

Interpolation, mapping and assemblyStrain, stress and rotation measures

02 Essential equations and symbols

u(x,y,z)=u0−zry,v(x,y,z)=v0+zrx,w(x,y,z)=w0u(x,y,z)=u_0-zr_y,\quad v(x,y,z)=v_0+zr_x,\quad w(x,y,z)=w_0ε(z)=εm+zκ\boldsymbol\varepsilon(z)=\boldsymbol\varepsilon_m+z\boldsymbol\kappaKℓ=Km+Kb+Ks+Kd\mathbf K_\ell=\mathbf K_m+\mathbf K_b+\mathbf K_s+\mathbf K_dγξ(η)=12[(1−η)gA+(1+η)gC],γη(ξ)=12[(1−ξ)gD+(1+ξ)gB]\gamma_\xi(\eta)=\tfrac12[(1-\eta)g_A+(1+\eta)g_C],\quad \gamma_\eta(\xi)=\tfrac12[(1-\xi)g_D+(1+\xi)g_B]

Definitions and conventions

u0,v0,w0u_{0},v_{0},w_{0}
Midsurface translations.
rₓ,rᵧ
Section rotations about the local axes.
εm,κ,γs\varepsilon{}m, \kappa{}, \gamma{}s
Membrane strain, curvature and transverse shear.
A, Db, S
Section matrices for membrane, bending and shear response.
1234ACBDShear is interpolated from tying points24 DOFs: 6 at each of 4 nodesMembrane + bending + shear + drilling4 Gauss points (the 2 × 2 rule)
MITC4 tying points and local shell termsThe quadrilateral uses four transverse-shear tying locations (A, C, D, B). They build assumed covariant shear fields; they are distinct from the four stiffness-integration points.
5 detailed steps

03 Step by step derivation

Follow the assumptions, equations and discrete operations in order. The full explanation stays visible.

1 · Flat facet, local frame, and six DOFs per node

MITC4 is a four-node Reissner–Mindlin shell facet. Each corner carries local midsurface translations (u₀,v₀,w₀) and rotations (rₓ,rᵧ,r_z); the local basis is built from the facet diagonals and then nodal DOFs are transformed between local and global frames. Four nodes give 24 element DOFs.

Ni(ξ,η)=14(1+ξiξ)(1+ηiη)N_i(\xi,\eta)=\tfrac14(1+\xi_i\xi)(1+\eta_i\eta)u(x,y,z)=u0−zry,v(x,y,z)=v0+zrx,w(x,y,z)=w0u(x,y,z)=u_0-zr_y,\qquad v(x,y,z)=v_0+zr_x,\qquad w(x,y,z)=w_0qe=[u1,v1,w1,rx1,ry1,rz1,…]T∈R24\mathbf q_e=[u_1,v_1,w_1,r_{x1},r_{y1},r_{z1},\ldots]^T\in\mathbb R^{24}

Rotations remain independent of the slope of w. That permits transverse shear, unlike Kirchhoff–Love kinematics.

2 · Separate membrane, curvature, and shear energy

The in-plane strain changes linearly through thickness; integrating the isotropic plane-stress law gives extensional, bending, and shear resultants. With t the thickness, E and ν material constants, and κₛ a shear correction, membrane stiffness scales with t, bending stiffness with t³, and shear stiffness with t.

ε(z)=εm+zκ\boldsymbol\varepsilon(z)=\boldsymbol\varepsilon_m+z\boldsymbol\kappaDp=E1−ν2[1ν0ν1000(1−ν)/2]\mathbf D_p=\frac{E}{1-\nu^2}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&(1-\nu)/2\end{bmatrix}A=tDp,Db=t312Dp,S=κsGtI\mathbf A=t\mathbf D_p,\qquad \mathbf D_b=\frac{t^3}{12}\mathbf D_p,\qquad \mathbf S=\kappa_sGt\mathbf I

The sign of curvature and the top/bottom stress labels follow the oriented local normal; reversing the facet orientation changes those labels.

3 · MITC tying field and shear locking

Direct bilinear interpolation of transverse shear can over-constrain a thin shell as t/L becomes small. MITC4 samples covariant shear at four tying points: A and C for the ξ component, D and B for the η component. Linear interpolation across the opposite parent coordinate constructs the assumed shear field, which is then mapped to physical components by the surface Jacobian.

gξ(η)=12[(1−η)gA+(1+η)gC]g_{\xi}(\eta)=\tfrac12[(1-\eta)g_A+(1+\eta)g_C]gη(ξ)=12[(1−ξ)gD+(1+ξ)gB]g_{\eta}(\xi)=\tfrac12[(1-\xi)g_D+(1+\xi)g_B]γs=Bsqe\boldsymbol\gamma_s=\mathbf B_s\mathbf q_e

The tying operator is the central MITC step. It is not a general cure for every distortion or shell boundary-layer error; patches, thin-shell benchmarks, and mesh refinement remain necessary.

4 · Element stiffness, drilling stabilization, and integration

The local 24×24 stiffness is accumulated as separate membrane, bending, shear, and drilling terms at four 2×2 Gauss points. r_z is a drilling rotation rather than a physical thin-shell kinematic; a weak penalty stabilizes it and must not dominate the physical energy. Local stiffness is transformed back to global DOFs.

Ke=Km+Kb+Ks+Kd\mathbf K_e=\mathbf K_m+\mathbf K_b+\mathbf K_s+\mathbf K_dKm=∫ABmTABm dA,Kb=∫ABbTDbBb dA\mathbf K_m=\int_A\mathbf B_m^T\mathbf A\mathbf B_m\,dA,\quad \mathbf K_b=\int_A\mathbf B_b^T\mathbf D_b\mathbf B_b\,dAγd=rz−12(v,x−u,y),Kd=∫AkdBdTBd dA\gamma_d=r_z-\tfrac12(v_{,x}-u_{,y}),\qquad \mathbf K_d=\int_Ak_d\mathbf B_d^T\mathbf B_d\,dAKeglobal=TdofTKelocalTdof\mathbf K_e^{global}=\mathbf T_{dof}^T\mathbf K_e^{local}\mathbf T_{dof}

The implementation exposes stiffness components separately, making drilling-energy sensitivity inspectable. Four-point integration applies to this documented route.

5 · Mass, evidence, and limits

The consistent shell mass carries surface inertia ρt for translations and rotary inertia ρt³/12 for rₓ,rᵧ. Drilling rotation has zero physical inertia and is deliberately massless; a reducer handles those DOFs where needed. Verification separates membrane, bending and shear patches, thin-shell locking studies, Cook, Scordelis–Lo, pinched cylinder, modal and dynamic cases. The card's maturity remains tied to each named analysis/material/benchmark record.

Me=∫ANTdiag⁡(ρt,ρt,ρt,ρt3/12,ρt3/12,0)N dA\mathbf M_e=\int_A\mathbf N^T\operatorname{diag}(\rho t,\rho t,\rho t,\rho t^3/12,\rho t^3/12,0)\mathbf N\,dA∫Aρt dA=ρtA,∫Aρt312 dA=ρt3A12\int_A\rho t\,dA=\rho tA,\qquad \int_A\rho\frac{t^3}{12}\,dA=\rho\frac{t^3A}{12}

Known limits include flat-facet kinematics, large rotations, warping, follower pressure, and dependence on drilling penalty. No lumped MITC4 mass route is claimed by the cited formulation.

04 Algorithm and implementation

Source excerpt · v0.2.11

QF Solver keeps membrane, bending, transverse shear, and drilling contributions separate. The loop evaluates each at the four shell Gauss points; a laminate coupling contribution is included when supplied.

Code
for xi, eta in self.gauss_pts:
    mats = self.strain_matrices_local(coords_2d, xi, eta, tying)
    Km += mats.Bm.T @ active_material.membrane_matrix @ mats.Bm * mats.detJ
    Kb += mats.Bb.T @ active_material.bending_matrix @ mats.Bb * mats.detJ
    coupling = getattr(active_material, "coupling_matrix", None)
    if coupling is not None:
        Kmb += (mats.Bm.T @ coupling @ mats.Bb + mats.Bb.T @ coupling.T @ mats.Bm) * mats.detJ
    Ks += mats.Bs.T @ active_material.shear_matrix @ mats.Bs * mats.detJ
    if active_material.drilling_stiffness > 0.0:
        Kd += mats.Bd.T @ mats.Bd * active_material.drilling_stiffness * mats.detJ
Read the complete versioned implementation →
Implementation

05 Verification and limitations

ROUTE-DEPENDENT · separate modal, transient, and harmonic evidence

Cook, Scordelis–Lo, pinched-cylinder, modal-plate, Newmark, harmonic, and laminate cases have separate records.

Each reference supports its own model, observable and analysis route. Read the original contract and limitations before transferring a conclusion.

Shell benchmark reference setConsolidated qualification registryv0.2.11 verification index
Capability matrix →Recorded cases and original sources →Run a maintained example →
Recorded repository calculations

See this method in a calculation.

Each case keeps its mesh, measured quantities, source record and limitations attached to the result.

QF Solver and Code_Aster displacement UZ at point E across five NAFEMS R0031 meshes.VNV-COMP-NAFEMS-R0031-CODEASTER-004NAFEMS R0031 laminated shell correlationQF Solver UZ at point E: −1.064851 mmRead model, result and limits →
Browse every calculation case →

06 Sources and related sheets

Full derivationElement matrices and loadsOriginal four-node plate formulationGeneral shell formulation referenceShell benchmark reference setImplementationStatic inputModal inputConsolidated qualification registryv0.2.11 verification index

Continue the reasoning

MITC3+ · three-node triangular shell
Back to sheet contents ↑