Solid Mechanics

Local Stiffness Matrix for Combined Beam and Spar Element With Axial and Lateral Linear Resistance

Don Warrington FORTRAN 77, Solid Mechanics Leave a comment

Our objective is to develop an element with the following characteristics:

  • Two-dimensional, two node element
  • Euler-Bernoulli beam theory
  • Axial stiffness (“spar” type element)
  • “Beam on elastic foundation” characteristic
  • Axial elastic resistance

Although it is doubtless possible to start with a single weak-form equation and develop the stiffness matrix, it is more convenient to develop the axial and bending local stiffness matrices separately and then to put them together with superposition.

Both spar and beam elements generally use two nodes, one at each end. For this derivation all of the constants (beam elastic modulus, moment of inertia, cross-sectional area and spring constants) will be assumed to be uniform the full length of the element. If one desires to model non-uniform beams, one can either develop an element with the desired non-uniformity or use more elements, and we see both in finite element analysis.

Let us start with the bending portion. The weak form of the equation for the fourth-order Euler-Bernoulli beam element is
\int_{x_{{1}}}^{x_{{2}}}\!{\it E1}\,{\it XI1}\,\left({\frac{d^{2}}{d{x}^{2}}}v(x)\right){\frac{d^{2}}{d{x}^{2}}}w(x)+gv(x)w(x)-v(x)t{dx}=0

where E1 is the Young’s modulus of the material and XI1 is the moment of inertia of the beam. The variable g represents the continuous spring constant along the length of the beam relative to the displacement of that beam, the “beam on elastic foundation.” The variable t is a uniform load along the beam. The equations were derived using Maple with the idea of the results used on FORTRAN 77, thus the naming convention of some of the variables. An explanation of the weak form, its derivation and the significance of w\left(x\right) and v\left(x\right) can be found in a finite element text such as this.

At this point we need to select appropriate weighting functions for the equation. For beam elements we choose weighting functions to satisfy the Hermite interplation of the two primary variables at local nodes 1 and 2, to wit

\Delta_{{1}}=C_{{1}}+C_{{2}}x_{{1}}+C_{{3}}{x_{{1}}}^{2}+C_{{4}}{x_{{1}}}^{3}
\Delta_{{2}}=-C_{{2}}-2\, C_{{3}}x_{{1}}-3\, C_{{4}}{x_{{1}}}^{2}
\Delta_{{3}}=C_{{1}}+C_{{2}}x_{{2}}+C_{{3}}{x_{{2}}}^{2}+C_{{4}}{x_{{2}}}^{3}
\Delta_{{4}}=-C_{{2}}-2\, C_{{3}}x_{{2}}-3\, C_{{4}}{x_{{2}}}^{2}

where \Delta_{1},\,\Delta_{3} are the “displacements” for nodes
1 and 2 and \Delta_{2},\,\Delta_{4} are the first derivative slopes
at these nodes.

This can be expressed in matrix form as follows:

\left[\begin{array}{cccc} 1 & x_{{1}} & {x_{{1}}}^{2} & {x_{{1}}}^{3}\\ \noalign{\medskip}0 & -1 & -2\, x_{{1}} & -3\,{x_{{1}}}^{2}\\ \noalign{\medskip}1 & x_{{2}} & {x_{{2}}}^{2} & {x_{{2}}}^{3}\\ \noalign{\medskip}0 & -1 & -2\, x_{{2}} & -3\,{x_{{2}}}^{2} \end{array}\right]\left[\begin{array}{c} C_{{1}}\\ \noalign{\medskip}C_{{2}}\\ \noalign{\medskip}C_{{3}}\\ \noalign{\medskip}C_{{4}} \end{array}\right]=\left[\begin{array}{c} \Delta_{{1}}\\ \noalign{\medskip}\Delta_{{2}}\\ \noalign{\medskip}\Delta_{{3}}\\ \noalign{\medskip}\Delta_{{4}} \end{array}\right]

Inverting the matrix, we have

\left[\begin{array}{c} C_{{1}}\\ \noalign{\medskip}C_{{2}}\\ \noalign{\medskip}C_{{3}}\\ \noalign{\medskip}C_{{4}} \end{array}\right]=\left[\begin{array}{cccc} {\frac{\left(3\, x_{{1}}-x_{{2}}\right){x_{{2}}}^{2}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & {\frac{x_{{1}}{x_{{2}}}^{2}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}} & {\frac{\left(x_{{1}}-3\, x_{{2}}\right){x_{{1}}}^{2}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & {\frac{{x_{{1}}}^{2}x_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}-6\,{\frac{x_{{1}}x_{{2}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & -{\frac{\left(2\, x_{{1}}+x_{{2}}\right)x_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}} & 6\,{\frac{x_{{1}}x_{{2}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & -{\frac{\left(x_{{1}}+2\, x_{{2}}\right)x_{{1}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}3\,{\frac{x_{{1}}+x_{{2}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & {\frac{x_{{1}}+2\, x_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}} & -3\,{\frac{x_{{1}}+x_{{2}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}} & {\frac{2\, x_{{1}}+x_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}-2\,\left({x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}\right)^{-1} & -\left({x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}\right)^{-1} & 2\,\left({x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}\right)^{-1} & -\left({x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}\right)^{-1} \end{array}\right]\left[\begin{array}{c} \Delta_{{1}}\\ \noalign{\medskip}\Delta_{{2}}\\ \noalign{\medskip}\Delta_{{3}}\\ \noalign{\medskip}\Delta_{{4}} \end{array}\right]
Multiplying the result, we have for the coefficients
\left[\begin{array}{c} C_{{1}}\\ \noalign{\medskip}C_{{2}}\\ \noalign{\medskip}C_{{3}}\\ \noalign{\medskip}C_{{4}} \end{array}\right]=\left[\begin{array}{c} {\frac{\left(3\, x_{{1}}-x_{{2}}\right){x_{{2}}}^{2}\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{x_{{1}}{x_{{2}}}^{2}\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}+{\frac{\left(x_{{1}}-3\, x_{{2}}\right){x_{{1}}}^{2}\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{{x_{{1}}}^{2}x_{{2}}\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}-6\,{\frac{x_{{1}}x_{{2}}\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}-{\frac{\left(2\, x_{{1}}+x_{{2}}\right)x_{{2}}\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}+6\,{\frac{x_{{1}}x_{{2}}\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}-{\frac{\left(x_{{1}}+2\, x_{{2}}\right)x_{{1}}\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}3\,{\frac{\left(x_{{1}}+x_{{2}}\right)\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{\left(x_{{1}}+2\, x_{{2}}\right)\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}-3\,{\frac{\left(x_{{1}}+x_{{2}}\right)\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{\left(2\, x_{{1}}+x_{{2}}\right)\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}\\ \noalign{\medskip}-2\,{\frac{\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}-{\frac{\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}+2\,{\frac{\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}-{\frac{\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}} \end{array}\right]

The weighting function in its complete form is thus

w={\frac{\left(3\, x_{{1}}-x_{{2}}\right){x_{{2}}}^{2}\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{x_{{1}}{x_{{2}}}^{2}\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}+{\frac{\left(x_{{1}}-3\, x_{{2}}\right){x_{{1}}}^{2}\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}+{\frac{{x_{{1}}}^{2}x_{{2}}\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}+

-6\,{\frac{x_{{1}}x_{{2}}\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}x-{\frac{\left(2\, x_{{1}}+x_{{2}}\right)x_{{2}}\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}x+6\,{\frac{x_{{1}}x_{{2}}\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}x-{\frac{\left(x_{{1}}+2\, x_{{2}}\right)x_{{1}}\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}x

3\,{\frac{\left(x_{{1}}+x_{{2}}\right)\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}{x}^{2}+{\frac{\left(x_{{1}}+2\, x_{{2}}\right)\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}{x}^{2}-3\,{\frac{\left(x_{{1}}+x_{{2}}\right)\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}{x}^{2}+{\frac{\left(2\, x_{{1}}+x_{{2}}\right)\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}{x}^{2}

-2\,{\frac{\Delta_{{1}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}{x}^{3}-{\frac{\Delta_{{2}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}{x}^{3}+2\,{\frac{\Delta_{{3}}}{{x_{{1}}}^{3}-3\,{x_{{1}}}^{2}x_{{2}}+3\, x_{{1}}{x_{{2}}}^{2}-{x_{{2}}}^{3}}}{x}^{3}-{\frac{\Delta_{{4}}}{{x_{{2}}}^{2}-2\, x_{{1}}x_{{2}}+{x_{{1}}}^{2}}}{x}^{3}

This breaks down in to weighting functions for each independent variable as follows:
\Phi_{1}=1-3\,{\frac{{\it \bar{x}}^{2}}{{\it he}^{2}}}+2\,{\frac{{\it \bar{x}}^{3}}{{\it he}^{3}}}

\Phi_{2}=2\,{\frac{{\it \bar{x}}^{2}}{{\it he}}}-{\it \bar{x}}-{\frac{{\it \bar{x}}^{3}}{{\it he}^{2}}}

\Phi_{3}=3\,{\frac{{\it \bar{x}}^{2}}{{\it he}^{2}}}-2\,{\frac{{\it \bar{x}}^{3}}{{\it he}^{3}}}

\Phi_{4}=-{\frac{{\it \bar{x}}^{3}}{{\it he}^{2}}}+{\frac{{\it \bar{x}}^{2}}{{\it he}}}

additionally assuming that

\bar{x}=x-x_{1}
he=x_{2}-x_{1}

If we substitute these weighting functions into the weak form of the governing equations, perform the appropriate substitution, differentiation, integration and algebra, the first term results in the following stiffness matrix:

M=\left[\begin{array}{cccc} 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}} & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & -12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}} & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}\\ \noalign{\medskip}-6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}} & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}\\ \noalign{\medskip}-12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}} & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}} & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}\\ \noalign{\medskip}-6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}} & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}} & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}} \end{array}\right]

The second term (for the “elastic foundation”) yields the following stiffness matrix:
N=\left[\begin{array}{cccc} {\frac{13}{35}}\,{\it he}\, g & -{\frac{11}{210}}\,{\it he}^{2}g & {\frac{9}{70}}\,{\it he}\, g & {\frac{13}{420}}\,{\it he}^{2}g\\ \noalign{\medskip}-{\frac{11}{210}}\,{\it he}^{2}g & {\frac{1}{105}}\,{\it he}^{3}g & -{\frac{13}{420}}\,{\it he}^{2}g & -{\frac{1}{140}}\,{\it he}^{3}g\\ \noalign{\medskip}{\frac{9}{70}}\,{\it he}\, g & -{\frac{13}{420}}\,{\it he}^{2}g & {\frac{13}{35}}\,{\it he}\, g & {\frac{11}{210}}\,{\it he}^{2}g\\ \noalign{\medskip}{\frac{13}{420}}\,{\it he}^{2}g & -{\frac{1}{140}}\,{\it he}^{3}g & {\frac{11}{210}}\,{\it he}^{2}g & {\frac{1}{105}}\,{\it he}^{3}g \end{array}\right]

The combined local stiffness matrix for bending only is

K_{b}=\left[\begin{array}{cccc} 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{13}{35}}\,{\it he}\, g & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{11}{210}}\,{\it he}^{2}g & -12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{9}{70}}\,{\it he}\, g & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{13}{420}}\,{\it he}^{2}g\\ \noalign{\medskip}-6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{11}{210}}\,{\it he}^{2}g & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}+{\frac{1}{105}}\,{\it he}^{3}g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{13}{420}}\,{\it he}^{2}g & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}-{\frac{1}{140}}\,{\it he}^{3}g\\ \noalign{\medskip}-12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{9}{70}}\,{\it he}\, g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{13}{420}}\,{\it he}^{2}g & 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{13}{35}}\,{\it he}\, g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{11}{210}}\,{\it he}^{2}g\\ \noalign{\medskip}-6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{13}{420}}\,{\it he}^{2}g & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}-{\frac{1}{140}}\,{\it he}^{3}g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{11}{210}}\,{\it he}^{2}g & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}+{\frac{1}{105}}\,{\it he}^{3}g \end{array}\right]

The FORTRAN 77 code for this is as follows:

K(1,1) = 12/he**3*E1*XI1+13.E0/35.E0*he*g
K(1,2) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g
K(1,3) = -12/he**3*E1*XI1+9.E0/70.E0*he*g
K(1,4) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g
K(2,1) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g
K(2,2) = 4/he*E1*XI1+he**3*g/105
K(2,3) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g
K(2,4) = 2/he*E1*XI1-he**3*g/140
K(3,1) = -12/he**3*E1*XI1+9.E0/70.E0*he*g
K(3,2) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g
K(3,3) = 12/he**3*E1*XI1+13.E0/35.E0*he*g
K(3,4) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g
K(4,1) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g
K(4,2) = 2/he*E1*XI1-he**3*g/140
K(4,3) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g
K(4,4) = 4/he*E1*XI1+he**3*g/105

The vector for the last term is

T=\left[\begin{array}{c} 1/2\,{\it he}\, t\\ \noalign{\medskip}-1/12\,{\it he}^{2}t\\ \noalign{\medskip}1/2\,{\it he}\, t\\ \noalign{\medskip}1/12\,{\it he}^{2}t \end{array}\right]

and the FORTRAN for this is

te(1,1) = he*t/2
te(2,1) = -he**2*t/12
te(3,1) = he*t/2
te(4,1) = he**2*t/12

Now let us turn to the spar element part of the stiffness matrix. The weak form equation for this is
{\it E1}\, A\int_{0}^{{\it he}}\!\left({\frac{d}{dx}}w(x)\right){\frac{d}{dx}}y(x){dx}+w(x){\frac{d}{dx}}y(x)+\int_{0}^{{\it he}}\! w(x)cy(x){dx}-\int\! w(x)q{dx}=0

Here A is the cross-sectional area of the beam, c is a distributed axial spring constant along the spar, and q is a distributed axial force along the element. To integrate from 0 to he is no different than doing so from x_{1} to x{}_{2}, only the coordinates change.

In this case we select linear weighting functions, to wit

W_{1}=1-{\frac{x}{{\it he}}}
W_{2}={\frac{x}{{\it he}}}

If as before we do the substitutions and integrations, we end up with a local stiffness matrix for the spar element only as follows:
K_{s}=\left[\begin{array}{cc} {\frac{{\it E1}\, A}{{\it he}}}+1/3\, c{\it he} & -{\frac{{\it E1}\, A}{{\it he}}}+1/6\, c{\it he}\\ \noalign{\medskip}-{\frac{{\it E1}\, A}{{\it he}}}+1/6\, c{\it he} & {\frac{{\it E1}\, A}{{\it he}}}+1/3\, c{\it he} \end{array}\right]

FORTRAN code for this is

K1(1,1) = E1*A/he+c*he/3
K1(1,2) = -E1*A/he+c*he/6
K1(2,1) = -E1*A/he+c*he/6
K1(2,2) = E1*A/he+c*he/3

The right hand side vector is as follows:
T=\left[\begin{array}{c} 1/2\,{\it he}\, q\\ \noalign{\medskip}1/2\,{\it he}\, q \end{array}\right]

and the code for this is

fe1(1,1) = he*q/2
fe1(2,1) = he*q/2

Now we need to combine these. We note that there are three variables:

  • x displacement (spar element only)
  • y displacement (beam element only)
  • rotation (beam element only)

We thus construct a 6\times6 element with the rows and columns in the above order, repeated twice each way for the two nodes. Doing this results in the following local stiffness matrix:
K=\left[\begin{array}{cccccc} {\frac{{\it E1}\, A}{{\it he}}}+1/3\, c{\it he} & 0 & 0 & -{\frac{{\it E1}\, A}{{\it he}}}+1/6\, c{\it he} & 0 & 0\\ \noalign{\medskip}0 & 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{13}{35}}\,{\it he}\, g & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{11}{210}}\,{\it he}^{2}g & 0 & -12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{9}{70}}\,{\it he}\, g & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{13}{420}}\,{\it he}^{2}g\\ \noalign{\medskip}0 & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{11}{210}}\,{\it he}^{2}g & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}+{\frac{1}{105}}\,{\it he}^{3}g & 0 & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{13}{420}}\,{\it he}^{2}g & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}-{\frac{1}{140}}\,{\it he}^{3}g\\ \noalign{\medskip}-{\frac{{\it E1}\, A}{{\it he}}}+1/6\, c{\it he} & 0 & 0 & {\frac{{\it E1}\, A}{{\it he}}}+1/3\, c{\it he} & 0 & 0\\ \noalign{\medskip}0 & -12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{9}{70}}\,{\it he}\, g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}-{\frac{13}{420}}\,{\it he}^{2}g & 0 & 12\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{3}}}+{\frac{13}{35}}\,{\it he}\, g & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{11}{210}}\,{\it he}^{2}g\\ \noalign{\medskip}0 & -6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{13}{420}}\,{\it he}^{2}g & 2\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}-{\frac{1}{140}}\,{\it he}^{3}g & 0 & 6\,{\frac{{\it E1}\,{\it XI1}}{{\it he}^{2}}}+{\frac{11}{210}}\,{\it he}^{2}g & 4\,{\frac{{\it E1}\,{\it XI1}}{{\it he}}}+{\frac{1}{105}}\,{\it he}^{3}g \end{array}\right]

or in code

K2(1,1) = E1*A/he+c*he/3
K2(1,2) = 0
K2(1,3) = 0
K2(1,4) = -E1*A/he+c*he/6
K2(1,5) = 0
K2(1,6) = 0
K2(2,1) = 0
K2(2,2) = 12/he**3*E1*XI1+13.E0/35.E0*he*g
K2(2,3) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g
K2(2,4) = 0
K2(2,5) = -12/he**3*E1*XI1+9.E0/70.E0*he*g
K2(2,6) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g
K2(3,1) = 0
K2(3,2) = -6/he**2*E1*XI1-11.E0/210.E0*he**2*g
K2(3,3) = 4/he*E1*XI1+he**3*g/105
K2(3,4) = 0
K2(3,5) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g
K2(3,6) = 2/he*E1*XI1-he**3*g/140
K2(4,1) = -E1*A/he+c*he/6
K2(4,2) = 0
K2(4,3) = 0
K2(4,4) = E1*A/he+c*he/3
K2(4,5) = 0
K2(4,6) = 0
K2(5,1) = 0
K2(5,2) = -12/he**3*E1*XI1+9.E0/70.E0*he*g
K2(5,3) = 6/he**2*E1*XI1-13.E0/420.E0*he**2*g
K2(5,4) = 0
K2(5,5) = 12/he**3*E1*XI1+13.E0/35.E0*he*g
K2(5,6) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g
K2(6,1) = 0
K2(6,2) = -6/he**2*E1*XI1+13.E0/420.E0*he**2*g
K2(6,3) = 2/he*E1*XI1-he**3*g/140
K2(6,4) = 0
K2(6,5) = 6/he**2*E1*XI1+11.E0/210.E0*he**2*g
K2(6,6) = 4/he*E1*XI1+he**3*g/105

As long as all of the elements line up along the x-axis, we are done. But we know that this cannot always be the case. So we need to effect a rotation of the local stiffness matrix. Since each element can be either oriented differently, of different length or both, we need
to rotate the local stiffness matrix before inserting it into the global one. The rotation matrix is
G=\left[\begin{array}{cccccc} {\it cosine} & {\it sine} & 0 & 0 & 0 & 0\\ \noalign{\medskip}-{\it sine} & {\it cosine} & 0 & 0 & 0 & 0\\ \noalign{\medskip}0 & 0 & 1 & 0 & 0 & 0\\ \noalign{\medskip}0 & 0 & 0 & {\it cosine} & {\it sine} & 0\\ \noalign{\medskip}0 & 0 & 0 & -{\it sine} & {\it cosine} & 0\\ \noalign{\medskip}0 & 0 & 0 & 0 & 0 & 1 \end{array}\right]

where sine and cosine are the angles of the elements from the x-axis. To effect a rotation, we need to first premultiply the matrix $K$ by the inverse of G and then postmultiply the result by G. That process is somewhat simplified by the fact that G is orthogonal; thus, its inverse and transpose are identical. Going through that process, the rotated local stiffness matrix is (in code only; we have overwhelmed WordPress’ LaTex conversion capability):

Kglobal(1,1) = cosine**2*(E1*A/he+c*he/3)+sine**2*(12/he**3*E1*XI1
#+13.E0/35.E0*he*g)
Kglobal(1,2) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+
#13.E0/35.E0*he*g)*cosine
Kglobal(1,3) = -sine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)
Kglobal(1,4) = cosine**2*(-E1*A/he+c*he/6)+sine**2*(-12/he**3*E1*X
#I1+9.E0/70.E0*he*g)
Kglobal(1,5) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI
#1+9.E0/70.E0*he*g)*cosine
Kglobal(1,6) = -sine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)
Kglobal(2,1) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+
#13.E0/35.E0*he*g)*cosine
Kglobal(2,2) = sine**2*(E1*A/he+c*he/3)+cosine**2*(12/he**3*E1*XI1
#+13.E0/35.E0*he*g)
Kglobal(2,3) = cosine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)
Kglobal(2,4) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI
#1+9.E0/70.E0*he*g)*cosine
Kglobal(2,5) = sine**2*(-E1*A/he+c*he/6)+cosine**2*(-12/he**3*E1*X
#I1+9.E0/70.E0*he*g)
Kglobal(2,6) = cosine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)
Kglobal(3,1) = -sine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)
Kglobal(3,2) = cosine*(-6/he**2*E1*XI1-11.E0/210.E0*he**2*g)
Kglobal(3,3) = 4/he*E1*XI1+he**3*g/105
Kglobal(3,4) = -sine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)
Kglobal(3,5) = cosine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)
Kglobal(3,6) = 2/he*E1*XI1-he**3*g/140
Kglobal(4,1) = cosine**2*(-E1*A/he+c*he/6)+sine**2*(-12/he**3*E1*X
#I1+9.E0/70.E0*he*g)
Kglobal(4,2) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI
#1+9.E0/70.E0*he*g)*cosine
Kglobal(4,3) = -sine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)
Kglobal(4,4) = cosine**2*(E1*A/he+c*he/3)+sine**2*(12/he**3*E1*XI1
#+13.E0/35.E0*he*g)
Kglobal(4,5) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+
#13.E0/35.E0*he*g)*cosine
Kglobal(4,6) = -sine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)
Kglobal(5,1) = cosine*(-E1*A/he+c*he/6)*sine-sine*(-12/he**3*E1*XI
#1+9.E0/70.E0*he*g)*cosine
Kglobal(5,2) = sine**2*(-E1*A/he+c*he/6)+cosine**2*(-12/he**3*E1*X
#I1+9.E0/70.E0*he*g)
Kglobal(5,3) = cosine*(6/he**2*E1*XI1-13.E0/420.E0*he**2*g)
Kglobal(5,4) = cosine*(E1*A/he+c*he/3)*sine-sine*(12/he**3*E1*XI1+
#13.E0/35.E0*he*g)*cosine
Kglobal(5,5) = sine**2*(E1*A/he+c*he/3)+cosine**2*(12/he**3*E1*XI1
#+13.E0/35.E0*he*g)
Kglobal(5,6) = cosine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)
Kglobal(6,1) = -sine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)
Kglobal(6,2) = cosine*(-6/he**2*E1*XI1+13.E0/420.E0*he**2*g)
Kglobal(6,3) = 2/he*E1*XI1-he**3*g/140
Kglobal(6,4) = -sine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)
Kglobal(6,5) = cosine*(6/he**2*E1*XI1+11.E0/210.E0*he**2*g)
Kglobal(6,6) = 4/he*E1*XI1+he**3*g/105

The use of “sine” and “cosine” for the trigonometric functions makes it possible to compute these once for each matrix, thus speeding up computations.

One possible application of such a element is with driven piles or deep foundations in soil; the element can be used for both axial and flexural loads. The biggest problem is that the soil response is never linear, so they cannot be used in a “straightforward” fashion, but iteratively.

The Importance of Causality

Don Warrington Mathematics, Solid Mechanics Leave a comment

In his book Introduction to the Differential Equations of Physics, German physicist Ludwig Hopf opens with the following statement:

Any differential equation expresses a relation between derivatives or between derivatives and given functions of the variables.  It thus establishes a relation between the increments of certain quantities and these quantities themselves.  This property of a differential equation makes it the natural expression of the principle of causality which is the foundation of exact natural science.  the ancient Greeks established laws of nature in which certain relation between numbers (harmony of spheres) or certain shapes of bodies played a privileged role.  The law was supposed to state something about a process as a whole, or about the complete shape of a body.  In more recent times (Galileo, Newton, etc.) a different concept has been adopted.  We do not try to establish a relation between all phases of a process immediately, but only between one phase and the next.  A law of this type may express, for example, how a certain state will develop in the immediate future, or it may describe the influence of the state of a certain particle on the particles in the immediate neighbourhood.  Thus we have a procedure for the description of a law of nature in terms of small (mathematically speaking, infinitesimal) differences of time and space.  The increments with which the law is concerned appear as derivatives, i.e., as the limits of the quotient of the increments of the variables which describe the process over the increment of space or time in wihch this development takes place.  A law of nature of this form is the expression of the relation between one state and the neighbouring (in time or space) states and therefore represents a special form of the principle of causality.

The whole issue of causality is an important one for both scientific and theological reasons, and I want to touch on one of each.

Every event that takes place in the universe is a result of an event before it.  Those events in turn are the results of those which have gone before.  All of these events form a chain which leads back to the first cause.  The need for the first cause is one of St. Thomas Aquinas’ proofs of God’s existence:

The second way is from the nature of the efficient cause. In the world of sense we find there is an order of efficient causes. There is no case known (neither is it, indeed, possible) in which a thing is found to be the efficient cause of itself; for so it would be prior to itself, which is impossible. Now in efficient causes it is not possible to go on to infinity, because in all efficient causes following in order, the first is the cause of the intermediate cause, and the intermediate is the cause of the ultimate cause, whether the intermediate cause be several, or only one. Now to take away the cause is to take away the effect. Therefore, if there be no first cause among efficient causes, there will be no ultimate, nor any intermediate cause. But if in efficient causes it is possible to go on to infinity, there will be no first efficient cause, neither will there be an ultimate effect, nor any intermediate efficient causes; all of which is plainly false. Therefore it is necessary to admit a first efficient cause, to which everyone gives the name of God.

Although, as Hopf points out, our understanding of how that causality actually works in the physical universe is different from the Greeks (and Thomas Aquinas worked in a Greek concept of natural philosophy) the truth of the importance of causality is undiminished.

To determine what comes after is a major reason for differential equations, which contain three elements: the equation itself, the initial conditions and the boundary conditions.  Once we have these, we can predict the behaviour of a system.  In some cases we can do so with a “simple” equation, others require discretisation and numerical modelling.  And that leads to our second point.

It’s interesting that Hopf speaks of “exact natural science.”  Today much of science and engineering is driven by probabalistic considerations, which in turn lead to statistical analysis.  Probability and statistics is a very useful tool, but not a substitute for the understanding of the actual mechanisms by which things work.  The actual mechanisms (physical laws, etc.) are what cause the phenomena which we record as statistics, not the other way around.  The fact that there are variations in these should not blind us to the core reality.

The advent of computers with broad-based number crunching abilities has only inflated our overconfidence in such methods.  It is essential, however, that we understand the why of phenomena as well as the what.  We must both be able to quantify the results and the correct causes of what is going on around us.  Two recent debacles illustrate this.

The first is the climate change fiasco we’ve been treated to of late.  Removing the dissimulation (as opposed to simulation) of some involved in the science, the core problem is that we do not as of yet have a model of global climate sufficiently comprehensive so that we can dispense with reliance on the statistics and project what will happen with a reasonable degree of confidence.  Part of the problem is the core problem in chaos theory: minor variations in initial conditions lead to major variations in the results.  But without such a model we are bereft with a definitive “why” as much as “what.”

The second is our financial collapse.  The models developed of the elaborate credit structure were fine as far as they went.  But ultimately they were divorced from sustainable reality because they did not take in to consideration all of the factors, many of which were obvious to those with raw experience.

The issue of causality is one that is central to our understanding of the universe.

Direct Derivation of the Equation of Motion for an Undamped Oscillating System in Phase Angle Form

Don Warrington Mathematics, Solid Mechanics 1 Comment

The equations of motion for linear vibrating systems are well known and widely used in both mechanical and electrical devices. However, when students are introduced to these, they are frequently presented with solutions which are either essentially underived or inadequately so.

This brief presentation will attempt to address this deficiency and hopefully show the derivation of the equation of motion for an undamped oscillating system in a more rigourous way.

Consider a simple spring/mass system without a forcing function. The equation of motion can be expressed as

1

where x(t) is displacement as a function of time, m is the mass of the system, and k is the spring constant. The negative sign on the right hand side of the equation is not an accident, as the spring force always opposes the motion of the mass, and is the result of using a mechanical engineers’ “free body diagram” method to develop the equation.

Solutions to this equation generally run in two forms. The first is a sum of sines and cosines:

2

But it’s more common to see it in the form of

3

The latter is simpler and easier to apply; however, it is seldom derived as much as assumed. So how can it be obtained from the original equation?

Let us begin by considering the original differential equation. With its constant coefficients, the most straightforward solution would be a solution where the derivative (and we, of course, would derive it twice) would be itself. This is the case where the function is exponential, so let us assume the equation to be in the form of

4

(I had an interesting fluid mechanics/heat transfer teacher who would say about this step that “you just write the answer down,” which we as his students found exasperating, but this method minimises that.)

Substituting this into the original equation of motion and diving out the identical exponentials yields

5

Solving for α yields

6

The right hand term is the natural frequency of the system, more generally expressed as a real number:

7

Thus for simplicity the solution can be written as

8

At this point it is not clear which of these two solutions is correct, so let us write the general solution as

9

Because of the complex exponential definition of sines and cosines, we see the beginning of a solution in simply one or the other, but at this point the coefficients are in the way.

These coefficients are determined from the initial conditions. Let us consider these at t=0:

10

Substituting these into our assumed general solution yields

11

The coefficients then solve to

12

It is noteworthy that the two coefficients are complex conjugates of each other.

Since the general solution is written in exponential form, it makes sense that, if the coefficients are to be removed so we can enable a direct solution to a sine or cosine, they too should be in exponential form. Converting the two coefficients to polar form yields

13

Substituting these coefficients into the general solution, we have

14

Factoring out the radical and recognising that the arctangent is an odd function,

15

The quantity in brackets is the complex exponential definition of the cosine, since the two exponents are negatives of each other. The solution can thus be written as

16

If we define

17

the solution is

18

which can be rewritten in a number of ways.

If the dampening is added, the problem can be solved in the same way, but the algebra is a little more complicated, and we will end up additionally with a real exponential (decay) in the final solution.

This derivation demonstrates the power of complex analysis as applied to differential equations even in a simple way.

More examples of this kind of thing are here and here.

Taking the Last Voyage with Newton and Pascal

Don Warrington Fluid Mechanics, Solid Mechanics 1 Comment

He’s not widely known outside of the fields he specialised in, but Adhémar Jean Claude Barré de Saint-Venant (1797-1886, usually known in the Anglophone world as simply Saint-Venant) was one of the premier scientists, engineers and mathematicians of the nineteenth century.  His accomplishments were many and include the following:

  • saint-venant-torsionSuccessful derivation of the Navier-Stokes Equations for a viscous flow before Stokes; these equations are the basis of computational fluid dynamics and the analysis of things that fly.
  • Systematisation and development of methods in the theory of elasticity of solids, including his semi-inverse methods for torsion, important in things such as automobile crank shafts.
  • Methods for the analysis of wave mechanics in bars, which we see in many places, from musical instruments to driven foundation piles.

Saint-Venant was born into a royalist, aristocratic, traditionally Roman Catholic family at a time when it was not safe to be any of these: the French Revolution, at that point stumbling from the Reign of Terror to control of France–and most of Europe–by Napoleon Bonaparte.  It was about the latter where Saint-Venant made a statement about himself that got him into trouble with the “new” Europe.  As described in S. Timoshenko’s History of Strength of Materials:

The political events of 1814 had a great effect on Saint-Venant’s career.  In March of that year, the armies of the allies were approaching Paris and the students of the École Polytechnique were mobilized.  On March 30, 1814, they were moving their guns to the Paris fortification when Saint-Venant, who was the first sargeant of the detachment, stepped out from the ranks with the exclamation: “My conscience forbids me to fight for an usurper…” His schoolmates resented that action very much and Saint-Venant was proclaimed a deserter and never allowed to resume his study at the École Polytechnique.

Saint-Venant’s statement of conscience was at once a political and religious statement, and “progressives” of his day didn’t miss either.  The French, then and now innocent of anti-discrimination legislation or sentiments, made his life miserable. The École Polytechnique was and is France’s premier technical institute of higher learning; getting kicked out of it was the equivalent of, say, being expelled from Princeton or MIT.  He worked in the powder industry for nine years, then was admitted to the École des Ponts et Chausées, where his fellow students shunned him.  He graduated first in his class anyway and began his illustrious career in technical things both theoretical and practical.

In spite of his difficulties within France, his reputation outside of her was another matter.  When François Napoleon Moigno wrote his book on statics, he discovered the following:

He (Moigno) wanted the portion on the statics of elastic bodies to be written by an expert in the theory of elasticity, but every time he asked for the collaboration of an English or a German scientist, he was given the same answer: “You have there, close to you, the authority par excellence, M. de Saint-Venant, consult him, listen to him, follow him.” One of them, M. Ettingshausen, added: “Your Academy of Sciences makes a mistake, a great mistake when it does not open its doors to a mathematician who is so highly placed in the opinion of the most competent judges.” In conclusion Moigno observes: “Fatally belittled in France of which he is the purest mathematical glory, M. de Saint-Venant enjoys a reputation in foreign countries which we dare to call grandiose.”

The French finally broke down and admitted Saint-Venant into the Academy of Sciences in 1868.  He continued his work, much of it from his home, up until the time of his death.  When the President of the Academy announced that passing, he made the following statement:

Old age was kind to our great colleague.  He died, advanced in years, without infirmities, occupied up to the last hour with problems which were dear to him and supported in the great passage by the hopes which had supported Pascal and Newton.

Europeans of the time would not have missed the import of the last statement: Pascal and Newton were Christians, and Saint-Venant was being identified with them as one also.  It was also a statement that Saint-Venant, for all of his achievements and interests which have enriched the world, also had an eternal goal as well.

There’s no evidence that Saint-Venant was ostentatious in his faith walk; descriptions of his life show the contrary.  And–shock to today’s atheist–there’s no evidence that it ever impeded the progress of his research or his thought.  As the statistician and eugenicist Karl Pearson, no friend of Christianity, noted:

The more I studied Saint-Venant’s work, the more new directions it seemed to me to open up for original investigation of the most valuable kind. It suggested innumerable unsolved problems in atomic physics, in impact, in plasticity and in a variety of other branches of elasticity, which do not seem beyond solution, and the solution of which if obtained would be of extreme importance. I felt convinced that a study of Saint-Venant’s researches would be a most valuable directive to the several young scientists, whose recent memoirs shew their interest in elasticity as well as their mathematical capacity. Many of the problems raised by Saint-Venant’s suggestive memoirs were quite beyond my powers of analysis, and I recognised that the most useful task I could undertake, was by a careful account of the memoirs themselves to lead the more competent on to their solution.

saint-venant-wave-equation

The biggest impediment he had to face was the blowback from his stand at the École Polytechnique, and that came from his secularist colleagues.  But, when the end came, all of his colleagues knew where he stood, in this life and the next one.

I spend a lot of time on this site and others talking about sea (and sometimes air) voyages.  And I’ve spent most of my career (and all the academic part of same) in the applied sciences.  But when I take my last voyage into eternity, I want to do it in the same hope of Newton and Pascal–and Saint-Venant and Euler for that matter–namely that which comes from following Jesus Christ out of the grave and into eternal life.

And you should too.

Note: I have a more complete account of Saint Venant’s story with a list of references here.