3.1.3 Example: Simply supported beam with varying cross-sectional area
In this example we will study a simply supported beam with varying cross-sectional area and refer to the figure in our previous example (3.1.2 Example: Simply supported beam with constant cross-sectional area). The difference with the previous example is that the second area moments a function of \( X \) due to the varying cross-sectional area.
We denote the second area moment at \( X=0 \) as \( I_0 \), and let the second area moment \( I \) vary in the following manner: $$ \begin{equation} I(X)=\frac{I_0}{1+(X/L)^n},\ n=2, 4, 6, \dots \tag{3.42} \end{equation} $$
Consider a rectangular cross-sectional area with a constant width \( b \) and a varying height \( h \): $$ \begin{equation} h(X)=\frac{h_0}{[1+(X/L)^n]^{1/3}} \tag{3.43} \end{equation} $$
where \( h_0 \) denotes \( h \) at \( X=0 \). From equation (3.26) in example (3.1.2 Example: Simply supported beam with constant cross-sectional area) we have: $$ \begin{align} \frac{d^2U}{dX^2}+\frac{P}{EI}U &=-\frac{q}{2EI}(L^2-X^2) \tag{3.44}\\ U(-L)&=U(L)=0,\quad \frac{dU(0)}{dX}=0 \tag{3.45} \end{align} $$
We start by computing the moment distribution \( M(X) \) in the beam, which is related with the displacement by the following expression: $$ \begin{equation} \frac{d^2U}{dX^2}=-\frac{M}{EI} \tag{3.46} \end{equation} $$
By substitution of equation (3.46) in (3.44) and two subsequent differentiations we get: $$ \begin{equation} \frac{d^2M}{dX^2}+\frac{P}{EI}@, M=-q \tag{3.47} \end{equation} $$
To represent the resulting second order ODE on a more convenient and generic form we introduce the dimensionless variables: $$ \begin{equation} x=\frac{X}{L},\ m=\frac{M}{qL^2} \tag{3.48} \end{equation} $$ and let $$ \begin{equation} P=\frac{EI_0}{L^2} \tag{3.49} \end{equation} $$
such that equation (3.47) becomes: $$ \begin{equation} m''(x)+(1+x^n)\cdot m(x)=-1 \tag{3.50} \end{equation} $$
with the boundary conditions: $$ \begin{equation} m(-1)=m(1)=0,\ m'(0)=0 \tag{3.51} \end{equation} $$
To our knowledge, no analytical solution exists for equation (3.50), even if it is linear and analytical solutions exists for the previous example in (3.1.2 Example: Simply supported beam with constant cross-sectional area).
Once the moment distribution has been computed, the dimensionless displacement \( u(x) \) may retrieved from: $$ \begin{equation} u(x)=m(x)-\frac{1}{2}(1-x^2) \tag{3.52} \end{equation} $$
The displacement in physical dimensions \( U \) may subsequently computed from the dimensionless \( u \): $$ \begin{equation} U=\frac{qL^4}{EI_0} u \tag{3.53} \end{equation} $$
Further, the differential equation for \( u(x) \) may be found by substitution of equation (3.52) into (3.50): $$ \begin{equation} u''(x)+(1+x^n)\cdot u(x)=-\frac{1}{2}(1-x^2)\cdot (1+x^n) \tag{3.54} \end{equation} $$
By exploiting the symmetry in the problem, we solve equation (3.50) numerically with a shooting method. The boundary conditions are: $$ \begin{equation} m(1)=0,\qquad m'(0)=0 \tag{3.55} \end{equation} $$
The dimensionless displacement \( u(x) \) is subsequently computed from
equation (3.52). We have implemented the outlined approach
in beam_deflect_shoot_varying.py in such a way that one easily may
choose between various RK-solvers (euler, heun, and rk4) from the
module ODEschemes.
We represent the second order ODE (3.50) as a system of first order ODEs by introducing the following conventions \( m(x)=y_1(x) \) and \( m'(x)=y_2(x) \): $$ \begin{align} y_1' & =y_2 \tag{3.56}\\ y_2' & =-\left (1+(1+x^n)\cdot y_1 \right ) \nonumber \end{align} $$ and the boundary conditions become: $$ \begin{equation} y_2(0)=y_1(1)=0 \tag{3.57} \end{equation} $$
In this case \( y_2(0)\equiv m'(0) \) is given, which leaves \( m(0)\equiv y_1(0) \) to be guessed such that the boundary condition at the other end \( m(1)\equiv y_1(1)=0 \) is fulfilled. To express the approach algorithmically as previously, we let \( s=y_1(0) \) and let the boundary value error function be \( \phi (s)=y_1(1;s) \).
- Let the two initial values \( s^0=0 \) and \( s^1=1 \) for simplicity
- Compute the corresponding values for the error function \( \phi ^0 \) and \( \phi ^1 \) by solving the initial value problem (in the current example (3.56)).
- Find the correct initial value \( s^* \) by linear interpolation:
The simplified expression for \( s \) in equation (3.58) may be found from (3.14) by setting \( s^0=0 \) and \( s^1=1 \).
One should be aware of that the choice in (3.58) is not necessarily always a good choice even though the ODE is linear. If the solution is of the kind \( e^{\alpha x} \) when both \( \alpha \) and \( x \) are large, we may end up outside the allowable range even for double precision which is approximately \( 10^{308} \). In our example above this is not a problem and we may use (3.58) and we have implemented the approach in beam_deflect_shoot_varying.py as shown below where both the moment and deflection is computed for \( n=2 \).
# src-ch2/beam_deflect_shoot_varying.py;ODEschemes.py @ git@lrhgit/tkt4140/src/src-ch2/ODEschemes.py;
from ODEschemes import euler, heun, rk4
from numpy import cos, sin
import numpy as np
from matplotlib.pyplot import *
# change some default values to make plots more readable
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
def f(y, x):
"""Governing differential equation on beam bending as a system of 1. order equations.
Args:
y(array): an array containg y and its derivatives up to second order. (RHS)
x(array): space array
Returns:
dydx(array): dydx. RHS of reduced system of 1st order equations
"""
yout = np.zeros_like(y)
yout[:] = [y[1],-(1+(1+x**n)*y[0])]
return yout
# === main program starts here ===
N = 10 # number of elements
L = 1.0 # half the length of the beam
x = np.linspace(0,L,N+1) # vector of
theta = 1 # PL**2/EI
theta2 = theta**2
h0=1 # height of beam at x = 0
n=2 # polynomial order to go into h(x)
h=h0/(1+(x/L)**n)**(1/3.) # height of beam assuming constant width
solvers = [euler, heun, rk4]
solver = solvers[2]
s = [0, 1] # guessed values
# === shoot ===
y01 = [s[0], 0] # initial values
y02 = [s[1], 0]
m0 = solver(f, y01, x)
m1 = solver(f, y02, x)
phi0 = m0[-1,0]
phi1 = m1[-1,0]
sfinal = phi0/(phi0 - phi1) # find correct initial value of moment m(x=0) using secant method
y03 = [sfinal, 0]
mfinal = solver(f,y03,x)
m = mfinal[:, 0] # extract moment from solution data
u = m -0.5*(1 - x**2) # final solution of dimensionless deflection
ua = (1/theta2)*(cos(theta*x)/cos(theta)-1)-(1-x**2)/2 #analytical solution with constant stifness
legendList=[] # empty list to append legends as plots are generated
plot(x, m, 'b',)
legendList.append('m')
plot(x, u, 'r',)
legendList.append('u')
## Add the labels
legend(legendList,loc='best',frameon=False)
ylabel('m, u')
xlabel('x')
grid(b=True, which='both', axis='both',linestyle='-')
show()
Figure 36: Moment \( m \) and deflection \( u \) .
