FEA TipsFinite Element Theory & MATLAB Programs
Home/Chapter 24/Addendum 24a

Addendum 24a

Corrected curvature expressions

Chapter 24 uses the practical simplified curvature expressions. This addendum presents the complete expressions including second-order terms, together with their derivatives and the deriv_corrected.m routine.

01

Corrected curvature expressions

The exact expressions of the curvatures, including the second-order terms, are:

 κ={κx‾κy‾κz‾}={dψx‾ds+12(dφz‾dsφy‾−dφy‾dsφz‾)dψy‾ds+12(dφx‾dsφz‾−dφz‾dsφx‾)dψz‾ds+12(dφy‾dsφx‾−dφx‾dsφy‾)}={dψx‾dsdψy‾dsdψz‾ds}+{12(dφz‾dsφy‾−dφy‾dsφz‾)12(dφx‾dsφz‾−dφz‾dsφx‾)12(dφy‾dsφx‾−dφx‾dsφy‾)}⏟Secondorderterms\, \htmlClass{fea-math-bold-italic}{\boldsymbol{\kappa =}}\left\{ \begin{matrix}{\kappa }_{\overline{x}}\\{\kappa }_{\overline{y}}\\{\kappa }_{\overline{z}}\end{matrix}\right\} \htmlClass{fea-math-bold-italic}{\boldsymbol{=}}\left\{ \begin{matrix}\frac{d{\psi }_{\overline{x}}}{ds}+\frac{1}{2}\left( \frac{d{\varphi }_{\overline{z}}}{ds}{\varphi }_{\overline{y}}-\frac{d{\varphi }_{\overline{y}}}{ds}{\varphi }_{\overline{z}}\right) \\\frac{d{\psi }_{\overline{y}}}{ds}+\frac{1}{2}\left( \frac{d{\varphi }_{\overline{x}}}{ds}{\varphi }_{\overline{z}}-\frac{d{\varphi }_{\overline{z}}}{ds}{\varphi }_{\overline{x}}\right) \\\frac{d{\psi }_{\overline{z}}}{ds}+\frac{1}{2}\left( \frac{d{\varphi }_{\overline{y}}}{ds}{\varphi }_{\overline{x}}-\frac{d{\varphi }_{\overline{x}}}{ds}{\varphi }_{\overline{y}}\right) \end{matrix}\right\} \htmlClass{fea-math-bold-italic}{\boldsymbol{=}}\left\{ \begin{matrix}\frac{d{\psi }_{\overline{x}}}{ds}\\\frac{d{\psi }_{\overline{y}}}{ds}\\\frac{d{\psi }_{\overline{z}}}{ds}\end{matrix}\right\} \htmlClass{fea-math-bold-italic}{\boldsymbol{+}}{\underbrace{\left\{ \begin{matrix}\frac{1}{2}\left( \frac{d{\varphi }_{\overline{z}}}{ds}{\varphi }_{\overline{y}}-\frac{d{\varphi }_{\overline{y}}}{ds}{\varphi }_{\overline{z}}\right) \\\frac{1}{2}\left( \frac{d{\varphi }_{\overline{x}}}{ds}{\varphi }_{\overline{z}}-\frac{d{\varphi }_{\overline{z}}}{ds}{\varphi }_{\overline{x}}\right) \\\frac{1}{2}\left( \frac{d{\varphi }_{\overline{y}}}{ds}{\varphi }_{\overline{x}}-\frac{d{\varphi }_{\overline{x}}}{ds}{\varphi }_{\overline{y}}\right) \end{matrix}\right\} }}_{\mathit{Second order terms}}

For the isoparametric two-node beam element, for which the rotations ψxˉ\psi_{\bar x}, ψyˉ\psi_{\bar y}, and ψzˉ\psi_{\bar z} have a linear variation along the element, we can write numerically:

 κ={κx‾κy‾κz‾}≅{ψx‾2−ψx‾1L+12(ψz‾2−ψz‾1L ψy‾2+ψy‾12−ψy‾2−ψy‾1L ψz‾2+ψz‾12)ψy‾2−ψy‾1L+12(ψx‾2−ψx‾1L ψz‾2+ψz‾12−ψz‾2−ψz‾1L ψx‾2+ψx‾12)ψz‾2−ψz‾1L+12(ψy‾2−ψy‾1L ψx‾2+ψx‾12−ψx‾2−ψx‾1L ψy‾2+ψy‾12)}\, \htmlClass{fea-math-bold-italic}{\boldsymbol{\kappa =}}\left\{ \begin{matrix}{\kappa }_{\overline{x}}\\{\kappa }_{\overline{y}}\\{\kappa }_{\overline{z}}\end{matrix}\right\} \htmlClass{fea-math-bold-italic}{\boldsymbol{\cong }}\left\{ \begin{matrix}\frac{{{\psi }_{\overline{x}}}_{2}-{{\psi }_{\overline{x}}}_{1}}{L}+\frac{1}{2}\left( \frac{{{\psi }_{\overline{z}}}_{2}-{{\psi }_{\overline{z}}}_{1}}{L}\, \frac{{{\psi }_{\overline{y}}}_{2}+{{\psi }_{\overline{y}}}_{1}}{2}-\frac{{{\psi }_{\overline{y}}}_{2}-{{\psi }_{\overline{y}}}_{1}}{L}\, \frac{{{\psi }_{\overline{z}}}_{2}+{{\psi }_{\overline{z}}}_{1}}{2}\right) \\\frac{{{\psi }_{\overline{y}}}_{2}-{{\psi }_{\overline{y}}}_{1}}{L}+\frac{1}{2}\left( \frac{{{\psi }_{\overline{x}}}_{2}-{{\psi }_{\overline{x}}}_{1}}{L}\, \frac{{{\psi }_{\overline{z}}}_{2}+{{\psi }_{\overline{z}}}_{1}}{2}-\frac{{{\psi }_{\overline{z}}}_{2}-{{\psi }_{\overline{z}}}_{1}}{L}\, \frac{{{\psi }_{\overline{x}}}_{2}+{{\psi }_{\overline{x}}}_{1}}{2}\right) \\\frac{{{\psi }_{\overline{z}}}_{2}-{{\psi }_{\overline{z}}}_{1}}{L}+\frac{1}{2}\left( \frac{{{\psi }_{\overline{y}}}_{2}-{{\psi }_{\overline{y}}}_{1}}{L}\, \frac{{{\psi }_{\overline{x}}}_{2}+{{\psi }_{\overline{x}}}_{1}}{2}-\frac{{{\psi }_{\overline{x}}}_{2}-{{\psi }_{\overline{x}}}_{1}}{L}\, \frac{{{\psi }_{\overline{y}}}_{2}+{{\psi }_{\overline{y}}}_{1}}{2}\right) \end{matrix}\right\}

Finally, we get:

{κx‾=1L(ψx‾2−ψx‾1)+12L(ψy‾1ψz‾2−ψy‾2ψz‾1)κy‾=1L(ψy‾2−ψy‾1)+12L(ψz‾1ψx‾2−ψz‾2ψx‾1)κz‾=1L(ψz‾2−ψz‾1)+12L(ψx‾1ψy‾2−ψx‾2ψy‾1)\left\{ \begin{matrix}{\kappa }_{\overline{x}}=\frac{1}{L}\left( {{\psi }_{\overline{x}}}_{2}-{{\psi }_{\overline{x}}}_{1}\right) +\frac{1}{2L}\left( {{\psi }_{\overline{y}}}_{1}{{\psi }_{\overline{z}}}_{2}-{{\psi }_{\overline{y}}}_{2}{{\psi }_{\overline{z}}}_{1}\right) \\{\kappa }_{\overline{y}}=\frac{1}{L}\left( {{\psi }_{\overline{y}}}_{2}-{{\psi }_{\overline{y}}}_{1}\right) +\frac{1}{2L}\left( {{\psi }_{\overline{z}}}_{1}{{\psi }_{\overline{x}}}_{2}-{{\psi }_{\overline{z}}}_{2}{{\psi }_{\overline{x}}}_{1}\right) \\{\kappa }_{\overline{z}}=\frac{1}{L}\left( {{\psi }_{\overline{z}}}_{2}-{{\psi }_{\overline{z}}}_{1}\right) +\frac{1}{2L}\left( {{\psi }_{\overline{x}}}_{1}{{\psi }_{\overline{y}}}_{2}-{{\psi }_{\overline{x}}}_{2}{{\psi }_{\overline{y}}}_{1}\right) \end{matrix}\right.

Refer to Sections 34.1 and 34.3 to understand how these equations were derived.

02

First and second derivatives

Accordingly, the first two derivatives of κxˉ\kappa_{\bar x}, κyˉ\kappa_{\bar y}, and κzˉ\kappa_{\bar z} become:

dκx‾ du‾el=12L  {000−2ψz‾2−ψy‾20002−ψz‾1ψy‾1}\frac{d{\kappa }_{\overline{x}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}}=\frac{1}{{2L}^{\, }}\, \left\{ \begin{matrix}0\\0\\0\\-2\\{{\psi }_{\overline{z}}}_{2}\\{{-\psi }_{\overline{y}}}_{2}\\0\\0\\0\\2\\-{{\psi }_{\overline{z}}}_{1}\\{{\psi }_{\overline{y}}}_{1}\end{matrix}\right\}
dκy‾ du‾el=1L  {000−ψz‾2−2ψx‾2000ψz‾12−ψx‾1}\frac{d{\kappa }_{\overline{y}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}}=\frac{1}{{L}^{\, }}\, \left\{ \begin{matrix}0\\0\\0\\{{-\psi }_{\overline{z}}}_{2}\\-2\\{{\psi }_{\overline{x}}}_{2}\\0\\0\\0\\{{\psi }_{\overline{z}}}_{1}\\2\\{{-\psi }_{\overline{x}}}_{1}\end{matrix}\right\}
dκz‾ du‾el=1L  {000ψy‾2−ψx‾2−2000−ψy‾1ψx‾12}\frac{d{\kappa }_{\overline{z}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}}=\frac{1}{{L}^{\, }}\, \left\{ \begin{matrix}0\\0\\0\\{{\psi }_{\overline{y}}}_{2}\\{{-\psi }_{\overline{x}}}_{2}\\-2\\0\\0\\0\\{{-\psi }_{\overline{y}}}_{1}\\{{\psi }_{\overline{x}}}_{1}\\2\end{matrix}\right\}
d2κx‾ du‾el2=12L  [000000000000000010−100000−101]\frac{{d}^{2}{\kappa }_{\overline{x}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}^{2}}=\frac{1}{{2L}^{\, \, }}\left[ \begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & 1\\0 & & & & & & & & & & -1 & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & -1 & & & & & & \\0 & & & & 1 & & & & & & & \end{matrix}\right]
d2κy‾ du‾el2=12L  [000000000000000−10010000100−1 ]\frac{{d}^{2}{\kappa }_{\overline{y}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}^{2}}=\frac{1}{{2L}^{\, \, }}\left[ \begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & -1\\0 & & & & & & & & & & & \\0 & & & & & & & & & 1 & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & 1 & & & & & & \\0 & & & & & & & & & & & \\0 & & & -1 & \, & & & & & & & \end{matrix}\right]
d2κz‾ du‾el2=12L  [0000000000000001 0−1  000000−1 01  ]\frac{{d}^{2}{\kappa }_{\overline{z}}^{\, }}{d{\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}}_{el}^{2}}=\frac{1}{{2L}^{\, \, }}\left[ \begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & 1 & \, \\0 & & & & & & & & & -1 & \, \, & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & & & & & & & & & \\0 & & & -1 & \, & & & & & & & \\0 & & 1 & \, & \, & & & & & & & \end{matrix}\right]
03

MATLAB routine

The first two derivatives of ε0\varepsilon_0, β0yˉ\beta_{0\bar y}, β0zˉ\beta_{0\bar z}, κxˉ\kappa_{\bar x}, κyˉ\kappa_{\bar y}, and κzˉ\kappa_{\bar z}, as presented in Chapter 24, are computed in the subprogram deriv.m. The corrected expressions of the curvatures, as shown here, are introduced in the subprogram deriv_corrected.m.

The two subprograms, deriv.m and deriv_corrected.m, generally give the same results. Moreover, using the deriv.m subprogram, convergence can be better, but it provides good results only if the load is far enough from the critical buckling load, when the second-order terms of the curvatures become very important. More about buckling in Chapter 34.

deriv_corrected.m
%*** deriv_corrected ***
u1=s(1); v1=s(2); w1=s(3); psix1=s(4);   psiy1=s(5);   psiz1=s(6); 
u2=s(7); v2=s(8); w2=s(9); psix2=s(10);  psiy2=s(11);  psiz2=s(12); 

ep =((v1 - v2)^2 + (w1 - w2)^2 + (L - u1 + u2)^2)^(1/2)/L - 1;
ep1=[ -(L - u1 + u2)
           (v1 - v2)
           (w1 - w2)
                   0
                   0
                   0
       (L - u1 + u2)
          -(v1 - v2)
          -(w1 - w2)
                   0
                   0
                   0]/L2;
ep2=[[    1,    0,    0, 0, 0, 0,   -1,    0,    0, 0, 0, 0]
     [    0,    1,    0, 0, 0, 0,    0,   -1,    0, 0, 0, 0]
     [    0,    0,    1, 0, 0, 0,    0,    0,   -1, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]
     [   -1,    0,    0, 0, 0, 0,    1,    0,    0, 0, 0, 0]
     [    0,   -1,    0, 0, 0, 0,    0,    1,    0, 0, 0, 0]
     [    0,    0,   -1, 0, 0, 0,    0,    0,    1, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]
     [    0,    0,    0, 0, 0, 0,    0,    0,    0, 0, 0, 0]]/L2;

bey =- psiz1/2 - psiz2/2 - atan((v1 - v2)/(L - u1 + u2));
bey1=[    -(v1 - v2)/L^2
      -(L - u1 + u2)/L^2
                       0
                       0
                       0
                    -1/2
           (v1 - v2)/L^2
       (L - u1 + u2)/L^2
                       0
                       0
                       0
                    -1/2];
bey2=[[      0,     -1, 0, 0, 0, 0,      0,      1, 0, 0, 0, 0]
      [      1,      0, 0, 0, 0, 0,     -1,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      1, 0, 0, 0, 0,      0,     -1, 0, 0, 0, 0]
      [     -1,      0, 0, 0, 0, 0,      1,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]
      [      0,      0, 0, 0, 0, 0,      0,      0, 0, 0, 0, 0]]/L2;
 
bez =atan((w1 - w2)/(L - u1 + u2)) - psiy2/2 - psiy1/2;
bez1=[    (w1 - w2)/L^2
                      0
      (L - u1 + u2)/L^2
                      0
                   -1/2
                      0
         -(w1 - w2)/L^2
                      0
     -(L - u1 + u2)/L^2
                      0
                   -1/2
                      0];
bez2 =[[  0, 0,  1, 0, 0, 0,  0, 0, -1, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [ -1, 0,  0, 0, 0, 0,  1, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  0, 0, -1, 0, 0, 0,  0, 0,  1, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  1, 0,  0, 0, 0, 0, -1, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]
       [  0, 0,  0, 0, 0, 0,  0, 0,  0, 0, 0, 0]]/L2;
 
 
% Shear convention: beta_z = atan(w'/(1+u')) + psi_y.
bez=-bez;
bez1=-bez1;
bez2=-bez2;

cax =- (psix1 - psix2)/L + ((psiy1*psiz2)/2 - (psiy2*psiz1)/2)/L;
  
cax1 =[     0
            0
            0
         -1/L
  psiz2/(2*L)
 -psiy2/(2*L)
            0
            0
            0
          1/L
 -psiz1/(2*L)
  psiy1/(2*L)];
 
cax2 = [[ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,  1/(2*L)]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0, -1/(2*L),        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,       0, 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,        0,-1/(2*L), 0, 0, 0, 0,        0,        0]
        [ 0, 0, 0, 0,  1/(2*L),       0, 0, 0, 0, 0,        0,        0]];
 
cay =-(psiy1 - psiy2)/L - ((psix1*psiz2)/2 - (psix2*psiz1)/2)/L;
 
cay1 =[     0
            0
            0
 -psiz2/(2*L)
         -1/L
  psix2/(2*L)
            0
            0
            0
  psiz1/(2*L)
          1/L
 -psix1/(2*L)];
 
cay2 = [[ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,-1/(2*L)]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,  1/(2*L), 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,  1/(2*L), 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,       0, 0,        0, 0, 0, 0,        0, 0,       0]
        [ 0, 0, 0,-1/(2*L), 0,        0, 0, 0, 0,        0, 0,       0]];
 
caz =- (psiz1 - psiz2)/L + ((psix1*psiy2)/2 - (psix2*psiy1)/2)/L;
  
caz1 =[     0
            0
            0
  psiy2/(2*L)
 -psix2/(2*L)
         -1/L
            0
            0
            0
 -psiy1/(2*L)
  psix1/(2*L)
          1/L];
 
caz2 = [[ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,  1/(2*L), 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,-1/(2*L),        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,-1/(2*L), 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,  1/(2*L),       0, 0, 0, 0, 0,       0,        0, 0]
        [ 0, 0, 0,        0,       0, 0, 0, 0, 0,       0,        0, 0]];

The routine deriv_corrected.m is included in the MATLAB package available in Chapter 24.