01
Model and constitutive relations A double-clamped beam is loaded at its midpoint by the force F F F [2] .
Material Linear elastic
Strain Small
Displacement Moderately large
Figure 1. Double-clamped beam, geometry, load, and displacement notation.Points A A A B B B
{ u A = 0 v A = 0 φ A = 0 \left\{\begin{array}{c}u_A=0\\[4pt]v_A=0\\[4pt]\varphi_A=0\end{array}\right. ⎩ ⎨ ⎧ u A = 0 v A = 0 φ A = 0 { u B = 0 v B = 0 φ B = 0 \left\{\begin{array}{c}u_B=0\\[4pt]v_B=0\\[4pt]\varphi_B=0\end{array}\right. ⎩ ⎨ ⎧ u B = 0 v B = 0 φ B = 0 The displacement-dependent axial strain makes the problem nonlinear. For the Euler–Bernoulli beam model, the generalized strain and work-conjugate generalized stress vectors are:
ε = { ε 0 κ } , σ = { N M } \htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\left\{\begin{array}{c}\varepsilon_0\\[4pt]\kappa\end{array}\right\},\qquad \htmlClass{fea-math-bold-italic}{\boldsymbol{\sigma}}=\left\{\begin{array}{c}N\\[4pt]M\end{array}\right\} ε = { ε 0 κ } , σ = { N M } In this planar analysis, T T T T y T_y T y M M M M z M_z M z
Hooke's law is:
σ = D ε , D = [ E A 0 0 E I ] \htmlClass{fea-math-bold-italic}{\boldsymbol{\sigma}}=\htmlClass{fea-math-bold-italic}{\boldsymbol D}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}},\qquad \htmlClass{fea-math-bold-italic}{\boldsymbol D}=\left[\begin{array}{cc}EA&0\\[4pt]0&EI\end{array}\right] σ = D ε , D = [ E A 0 0 E I ] or, component by component,
N = E A ε 0 , M = E I κ N=EA\varepsilon_0,\qquad M=EI\kappa N = E A ε 0 , M = E I κ ε 0 \varepsilon_0 ε 0 axial strain of the neutral axis κ \kappa κ curvature of the neutral axis A , I A,\ I A , I cross-sectional area and second moment of area 02
Cubic beam finite element Figure 2. Two-node cubic beam finite element and its six nodal degrees of freedom.With the normalized coordinate ξ = x / L \xi=x/L ξ = x / L Section 5.1 [1] .
u ( x ) = [ 1 − ξ , 0 , 0 , ξ , 0 , 0 ] u e l = N u u e l u(x)=\left[1-\xi,\ 0,\ 0,\ \xi,\ 0,\ 0\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_u\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} u ( x ) = [ 1 − ξ , 0 , 0 , ξ , 0 , 0 ] u e l = N u u e l v ( x ) = [ 0 , 2 ξ 3 − 3 ξ 2 + 1 , ( ξ 3 − 2 ξ 2 + ξ ) L , 0 , − 2 ξ 3 + 3 ξ 2 , ( ξ 3 − ξ 2 ) L ] u e l = N v u e l v(x)=\left[0,\ 2\xi^3-3\xi^2+1,\ (\xi^3-2\xi^2+\xi)L,\ 0,\ -2\xi^3+3\xi^2,\ (\xi^3-\xi^2)L\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_v\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} v ( x ) = [ 0 , 2 ξ 3 − 3 ξ 2 + 1 , ( ξ 3 − 2 ξ 2 + ξ ) L , 0 , − 2 ξ 3 + 3 ξ 2 , ( ξ 3 − ξ 2 ) L ] u e l = N v u e l φ ( x ) = d v d x = [ 0 , 6 ξ 2 − 6 ξ L , 3 ξ 2 − 4 ξ + 1 , 0 , − 6 ξ 2 + 6 ξ L , 3 ξ 2 − 2 ξ ] u e l = N v , x u e l \varphi(x)=\frac{dv}{dx}=\left[0,\ \frac{6\xi^2-6\xi}{L},\ 3\xi^2-4\xi+1,\ 0,\ \frac{-6\xi^2+6\xi}{L},\ 3\xi^2-2\xi\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} φ ( x ) = d x d v = [ 0 , L 6 ξ 2 − 6 ξ , 3 ξ 2 − 4 ξ + 1 , 0 , L − 6 ξ 2 + 6 ξ , 3 ξ 2 − 2 ξ ] u e l = N v , x u e l u e l = { u 1 v 1 φ 1 u 2 v 2 φ 2 } \htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\left\{\begin{array}{c}u_1\\[3pt]v_1\\[3pt]\varphi_1\\[3pt]u_2\\[3pt]v_2\\[3pt]\varphi_2\end{array}\right\} u e l = ⎩ ⎨ ⎧ u 1 v 1 φ 1 u 2 v 2 φ 2 ⎭ ⎬ ⎫ Here L L L u e l \htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} u e l
03
Moderately large-displacement strain Figure 3. Differential geometry used to derive the axial strain.The exact axial strain of the neutral axis obtained from Figure 3 is:
ε 0 = ( d x + d u ) 2 + d v 2 − d x d x = ( 1 + d u d x ) 2 + ( d v d x ) 2 − 1 \varepsilon_0=\frac{\sqrt{(dx+du)^2+dv^2}-dx}{dx}=\sqrt{\left(1+\frac{du}{dx}\right)^2+\left(\frac{dv}{dx}\right)^2}-1 ε 0 = d x ( d x + d u ) 2 + d v 2 − d x = ( 1 + d x d u ) 2 + ( d x d v ) 2 − 1 ε 0 = 1 + 2 d u d x + ( d u d x ) 2 + ( d v d x ) 2 − 1 \varepsilon_0=\sqrt{1+2\frac{du}{dx}+\left(\frac{du}{dx}\right)^2+\left(\frac{dv}{dx}\right)^2}-1 ε 0 = 1 + 2 d x d u + ( d x d u ) 2 + ( d x d v ) 2 − 1 The square root is expanded as:
1 + a ≈ 1 + 1 2 a − 1 8 a 2 , a = 2 d u d x + ( d u d x ) 2 + ( d v d x ) 2 \sqrt{1+a}\approx1+\frac12a-\frac18a^2,\qquad a=2\frac{du}{dx}+\left(\frac{du}{dx}\right)^2+\left(\frac{dv}{dx}\right)^2 1 + a ≈ 1 + 2 1 a − 8 1 a 2 , a = 2 d x d u + ( d x d u ) 2 + ( d x d v ) 2 After terms of third or higher degree are neglected, the moderately large-displacement approximation is:
ε 0 ≈ d u d x + 1 2 ( d v d x ) 2 \varepsilon_0\approx\frac{du}{dx}+\frac12\left(\frac{dv}{dx}\right)^2 ε 0 ≈ d x d u + 2 1 ( d x d v ) 2 The generalized strain vector is therefore:
ε = { ε 0 κ } = { d u d x + 1 2 ( d v d x ) 2 d 2 v d x 2 } \htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}=\left\{\begin{array}{c}\varepsilon_0\\[5pt]\kappa\end{array}\right\}=\left\{\begin{array}{c}\dfrac{du}{dx}+\dfrac12\left(\dfrac{dv}{dx}\right)^2\\[10pt]\dfrac{d^2v}{dx^2}\end{array}\right\} ε = { ε 0 κ } = ⎩ ⎨ ⎧ d x d u + 2 1 ( d x d v ) 2 d x 2 d 2 v ⎭ ⎬ ⎫ The derivatives of the interpolation functions are:
d u d x = 1 L [ − 1 , 0 , 0 , 1 , 0 , 0 ] u e l = N u , x u e l \frac{du}{dx}=\frac1L\left[-1,\ 0,\ 0,\ 1,\ 0,\ 0\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{u,x}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} d x d u = L 1 [ − 1 , 0 , 0 , 1 , 0 , 0 ] u e l = N u , x u e l d v d x = [ 0 , 6 ξ 2 − 6 ξ L , 3 ξ 2 − 4 ξ + 1 , 0 , − 6 ξ 2 + 6 ξ L , 3 ξ 2 − 2 ξ ] u e l = N v , x u e l \frac{dv}{dx}=\left[0,\ \frac{6\xi^2-6\xi}{L},\ 3\xi^2-4\xi+1,\ 0,\ \frac{-6\xi^2+6\xi}{L},\ 3\xi^2-2\xi\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} d x d v = [ 0 , L 6 ξ 2 − 6 ξ , 3 ξ 2 − 4 ξ + 1 , 0 , L − 6 ξ 2 + 6 ξ , 3 ξ 2 − 2 ξ ] u e l = N v , x u e l κ = d 2 v d x 2 = [ 0 , 12 ξ − 6 L 2 , 6 ξ − 4 L , 0 , − 12 ξ + 6 L 2 , 6 ξ − 2 L ] u e l = N κ u e l \kappa=\frac{d^2v}{dx^2}=\left[0,\ \frac{12\xi-6}{L^2},\ \frac{6\xi-4}{L},\ 0,\ \frac{-12\xi+6}{L^2},\ \frac{6\xi-2}{L}\right]\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{\kappa}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} κ = d x 2 d 2 v = [ 0 , L 2 12 ξ − 6 , L 6 ξ − 4 , 0 , L 2 − 12 ξ + 6 , L 6 ξ − 2 ] u e l = N κ u e l Hence,
ε 0 = N u , x u e l + 1 2 ( N v , x u e l ) 2 , κ = N κ u e l \varepsilon_0=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{u,x}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}+\frac12\left(\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}\right)^2,\qquad \kappa=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{\kappa}\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el} ε 0 = N u , x u e l + 2 1 ( N v , x u e l ) 2 , κ = N κ u e l 04
Potential energy and nonlinear equilibrium The total potential energy, summed over all beam elements, is:
Π = 1 2 ∑ e l ∫ 0 L ε T D ε d x ⏟ U − u T F ⏟ V \mathit{\Pi}=\underbrace{\frac12\sum_{el}\int_0^L\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol D}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}\,dx}_{U}\;\underbrace{-\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol F}}_{V} Π = U 2 1 e l ∑ ∫ 0 L ε T D ε d x V − u T F Π = 1 2 ∑ e l ∫ 0 L ( E A ε 0 2 + E I κ 2 ) d x ⏟ U − u T F ⏟ V \mathit{\Pi}=\underbrace{\frac12\sum_{el}\int_0^L\left(EA\varepsilon_0^2+EI\kappa^2\right)dx}_{U}\;\underbrace{-\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol F}}_{V} Π = U 2 1 e l ∑ ∫ 0 L ( E A ε 0 2 + E I κ 2 ) d x V − u T F U U U V V V F \htmlClass{fea-math-bold-italic}{\boldsymbol F} F u \htmlClass{fea-math-bold-italic}{\boldsymbol u} u
The minimum-total-potential-energy condition is:
Ψ = ∂ Π ∂ u = 0 \htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}=\frac{\partial\mathit{\Pi}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}}=\htmlClass{fea-math-bold-italic}{\boldsymbol0} Ψ = ∂ u ∂ Π = 0 This condition produces the nonlinear algebraic system Ψ ( u ) = 0 \htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}(\htmlClass{fea-math-bold-italic}{\boldsymbol u})=\htmlClass{fea-math-bold-italic}{\boldsymbol0} Ψ ( u ) = 0 i i i
Δ u ( i ) = − [ K t ( u ( i ) ) ] − 1 Ψ ( u ( i ) ) \mathit{\Delta}\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i)}=-\left[\htmlClass{fea-math-bold-italic}{\boldsymbol K}_t\!\left(\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i)}\right)\right]^{-1}\htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}\!\left(\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i)}\right) Δ u ( i ) = − [ K t ( u ( i ) ) ] − 1 Ψ ( u ( i ) ) u ( i + 1 ) = u ( i ) + Δ u ( i ) \htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i+1)}=\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i)}+\mathit{\Delta}\htmlClass{fea-math-bold-italic}{\boldsymbol u}^{(i)} u ( i + 1 ) = u ( i ) + Δ u ( i ) The residual and tangent stiffness are:
Ψ = ∂ U ∂ u − F = ∑ e l f e l − F \htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}=\frac{\partial U}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}}-\htmlClass{fea-math-bold-italic}{\boldsymbol F}=\sum_{el}\htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}-\htmlClass{fea-math-bold-italic}{\boldsymbol F} Ψ = ∂ u ∂ U − F = e l ∑ f e l − F K t ( u ) = ∂ 2 U ∂ u 2 = ∑ e l ∂ 2 U e l ∂ u e l 2 = ∑ e l k e l \htmlClass{fea-math-bold-italic}{\boldsymbol K}_t(\htmlClass{fea-math-bold-italic}{\boldsymbol u})=\frac{\partial^2U}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}^2}=\sum_{el}\frac{\partial^2U_{el}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2}=\sum_{el}\htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el} K t ( u ) = ∂ u 2 ∂ 2 U = e l ∑ ∂ u e l 2 ∂ 2 U e l = e l ∑ k e l 05
Element residual and tangent stiffness The deformation energy of one finite element is:
U e l = 1 2 ∫ 0 L ε T D ε d x = 1 2 ∫ 0 L ( E A ε 0 2 + E I κ 2 ) d x U_{el}=\frac12\int_0^L\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol D}\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}\,dx=\frac12\int_0^L\left(EA\varepsilon_0^2+EI\kappa^2\right)dx U e l = 2 1 ∫ 0 L ε T D ε d x = 2 1 ∫ 0 L ( E A ε 0 2 + E I κ 2 ) d x Using Gaussian quadrature,
U e l = 1 2 ∑ G = 1 n G η G ( E A ε 0 G 2 + E I κ G 2 ) L 2 , ε 0 G = ε 0 ( ξ G ) , κ G = κ ( ξ G ) U_{el}=\frac12\sum_{G=1}^{n_G}\eta_G\left(EA\varepsilon_{0G}^{2}+EI\kappa_G^{2}\right)\frac L2,\qquad \varepsilon_{0G}=\varepsilon_0(\xi_G),\quad \kappa_G=\kappa(\xi_G) U e l = 2 1 G = 1 ∑ n G η G ( E A ε 0 G 2 + E I κ G 2 ) 2 L , ε 0 G = ε 0 ( ξ G ) , κ G = κ ( ξ G ) Two Gauss points are sufficient, n G = 2 n_G=2 n G = 2 0 ≤ ξ ≤ 1 0\leq\xi\leq1 0 ≤ ξ ≤ 1
ξ G = 1 2 ( 1 ± 3 3 ) , η G = 1 \xi_G=\frac12\left(1\pm\frac{\sqrt3}{3}\right),\qquad \eta_G=1 ξ G = 2 1 ( 1 ± 3 3 ) , η G = 1 The element unbalanced-force vector and tangent stiffness matrix follow by differentiation:
f e l = ∂ U e l ∂ u e l , k e l = ∂ f e l ∂ u e l = ∂ 2 U e l ∂ u e l 2 \htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}=\frac{\partial U_{el}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}},\qquad \htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el}=\frac{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}=\frac{\partial^2U_{el}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2} f e l = ∂ u e l ∂ U e l , k e l = ∂ u e l ∂ f e l = ∂ u e l 2 ∂ 2 U e l f e l = ∑ G = 1 2 ( E A ε 0 ∂ ε 0 ∂ u e l + E I κ ∂ κ ∂ u e l ) G L 2 \htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}=\sum_{G=1}^{2}\left(EA\varepsilon_0\frac{\partial\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}+EI\kappa\frac{\partial\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}\right)_G\frac L2 f e l = G = 1 ∑ 2 ( E A ε 0 ∂ u e l ∂ ε 0 + E I κ ∂ u e l ∂ κ ) G 2 L k e l = ∑ G = 1 2 [ E A ε 0 ∂ 2 ε 0 ∂ u e l 2 + E A ∂ ε 0 ∂ u e l ( ∂ ε 0 ∂ u e l ) T + E I κ ∂ 2 κ ∂ u e l 2 + E I ∂ κ ∂ u e l ( ∂ κ ∂ u e l ) T ] G L 2 \htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el}=\sum_{G=1}^{2}\left[EA\varepsilon_0\frac{\partial^2\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2}+EA\frac{\partial\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}\left(\frac{\partial\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}\right)^T+EI\kappa\frac{\partial^2\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2}+EI\frac{\partial\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}\left(\frac{\partial\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}\right)^T\right]_G\frac L2 k e l = G = 1 ∑ 2 [ E A ε 0 ∂ u e l 2 ∂ 2 ε 0 + E A ∂ u e l ∂ ε 0 ( ∂ u e l ∂ ε 0 ) T + E I κ ∂ u e l 2 ∂ 2 κ + E I ∂ u e l ∂ κ ( ∂ u e l ∂ κ ) T ] G 2 L The required first and second derivatives are:
∂ ε 0 ∂ u e l = N u , x T + ( N v , x T N v , x ) u e l , ∂ κ ∂ u e l = N κ T \frac{\partial\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{u,x}^{T}+\left(\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}\right)\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el},\qquad \frac{\partial\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{\kappa}^{T} ∂ u e l ∂ ε 0 = N u , x T + ( N v , x T N v , x ) u e l , ∂ u e l ∂ κ = N κ T ∂ 2 ε 0 ∂ u e l 2 = N v , x T N v , x , ∂ 2 κ ∂ u e l 2 = 0 6 × 6 \frac{\partial^2\varepsilon_0}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2}=\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x}^{T}\htmlClass{fea-math-bold-italic}{\boldsymbol N}_{v,x},\qquad \frac{\partial^2\kappa}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}^2}=\htmlClass{fea-math-bold-italic}{\boldsymbol0}_{6\times6} ∂ u e l 2 ∂ 2 ε 0 = N v , x T N v , x , ∂ u e l 2 ∂ 2 κ = 0 6 × 6 Element vectors and matrices are assembled in the standard way:
Ψ = ∑ e l f e l − F , K t = ∑ e l k e l = ∑ e l ∂ f e l ∂ u e l \htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}=\sum_{el}\htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}-\htmlClass{fea-math-bold-italic}{\boldsymbol F},\qquad \htmlClass{fea-math-bold-italic}{\boldsymbol K}_t=\sum_{el}\htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el}=\sum_{el}\frac{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}}{\partial\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el}} Ψ = e l ∑ f e l − F , K t = e l ∑ k e l = e l ∑ ∂ u e l ∂ f e l The iterative process stops when the displacement-correction norm is smaller than the prescribed tolerance.
06
Example 1 — double-clamped beam The beam has a rectangular cross-section w × t = 5 m m × 1 m m w\times t=5\ \mathit{mm}\times1\ \mathit{mm} w × t = 5 mm × 1 mm d = 100 m m d=100\ \mathit{mm} d = 100 mm E = 2 × 10 5 M P a E=2\times10^5\ \mathit{MPa} E = 2 × 1 0 5 MPa F = 100 N F=100\ \mathit N F = 100 N
v C = − 1.86555 m m , ∣ v C ∣ d = 0.0186555 , N F = 8.01210 v_C=-1.86555\ \mathit{mm},\qquad \frac{|v_C|}{d}=0.0186555,\qquad \frac NF=8.01210 v C = − 1.86555 mm , d ∣ v C ∣ = 0.0186555 , F N = 8.01210 Although ∣ v C ∣ / d < 0.02 |v_C|/d<0.02 ∣ v C ∣/ d < 0.02
Figure 6. Shear-force diagram for F = 100 N F=100\ \mathit N F = 100 N Figure 7. Bending-moment diagram for F = 100 N F=100\ \mathit N F = 100 N Figure 8. Constant axial-force ratio along the beam: N / F = 8.012 N/F=8.012 N / F = 8.012 07
Example 2 — cantilever Consider a cantilever with w × t = 5 m m × 1 m m w\times t=5\ \mathit{mm}\times1\ \mathit{mm} w × t = 5 mm × 1 mm L = 100 m m L=100\ \mathit{mm} L = 100 mm E = 2 × 10 5 M P a E=2\times10^5\ \mathit{MPa} E = 2 × 1 0 5 MPa F = 4 N F=4\ \mathit N F = 4 N
l i n e a r m o d e l : v m a x L = F L 2 3 E I = 0.160 p r e s e n t e d m e t h o d , 40 e l e m e n t s : v m a x L = 0.1424 e x a c t e l a s t i c a : v m a x L = 0.1183 \begin{aligned}\mathit{linear model:}\quad &\frac{v_{max}}L=\frac{FL^2}{3EI}=0.160\\[5pt]\mathit{presented method, 40 elements:}\quad &\frac{v_{max}}L=0.1424\\[5pt]\mathit{exact elastica:}\quad &\frac{v_{max}}L=0.1183\end{aligned} linearmodel : presentedmethod , 40elements : exactelastica : L v ma x = 3 E I F L 2 = 0.160 L v ma x = 0.1424 L v ma x = 0.1183 The error of the moderately large-displacement approximation relative to the exact elastica result is about 20%. The approximation is suitable when the slope φ \varphi φ 10 ∘ ≈ 0.175 r a d 10^\circ\approx0.175\ \mathit{rad} 1 0 ∘ ≈ 0.175 rad
φ 2 = ( d v d x ) 2 ≈ 0.03 \varphi^2=\left(\frac{dv}{dx}\right)^2\approx0.03 φ 2 = ( d x d v ) 2 ≈ 0.03 This term is then small compared with unity [2] . In the present example, φ m a x ≈ 0.21 r a d ≈ 12 ∘ \varphi_{max}\approx0.21\ \mathit{rad}\approx12^\circ φ ma x ≈ 0.21 rad ≈ 1 2 ∘
08
MATLAB and ANSYS files The MATLAB package contains the complete nonlinear analysis, result recovery, diagram generation, and a README. The code follows the element equations above: deriv.m evaluates ε 0 \varepsilon_0 ε 0 κ \kappa κ stiff.m integrates and assembles f e l \htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el} f e l k e l \htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el} k e l main.m performs the load-step and Newton–Raphson loops.
deriv.m Strain and curvature
ep = Nux*uel + 1/2*(Nvx*uel)^2;
ep1 = Nux' + (Nvx'*Nvx)*uel;
ep2 = Nvx'*Nvx;
ka = Nka*uel;
ka1 = Nka';MATLAB compatibility
Tested in MATLAB R2024a. The code uses standard MATLAB syntax and is expected to be compatible with newer MATLAB versions. No additional toolbox is required.
The ANSYS macro uses symmetry: it models half of the beam, imposes the midpoint symmetry conditions, and applies F / 2 = 50 N F/2=50\ \mathit N F /2 = 50 N 100 N 100\ \mathit N 100 N
PREVIOUS SECTION ← 5.3. Isoparametric Two-Node Beam Element — Timoshenko Beam Model NEXT SECTION 7.1. Cantilever Beam → 
