
## 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: $$$$I(X)=\frac{I_0}{1+(X/L)^n},\ n=2, 4, 6, \dots \tag{3.42}$$$$

Consider a rectangular cross-sectional area with a constant width $$b$$ and a varying height $$h$$: $$$$h(X)=\frac{h_0}{[1+(X/L)^n]^{1/3}} \tag{3.43}$$$$

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: $$$$\frac{d^2U}{dX^2}=-\frac{M}{EI} \tag{3.46}$$$$

By substitution of equation (3.46) in (3.44) and two subsequent differentiations we get: $$$$\frac{d^2M}{dX^2}+\frac{P}{EI}@, M=-q \tag{3.47}$$$$

To represent the resulting second order ODE on a more convenient and generic form we introduce the dimensionless variables: $$$$x=\frac{X}{L},\ m=\frac{M}{qL^2} \tag{3.48}$$$$ and let $$$$P=\frac{EI_0}{L^2} \tag{3.49}$$$$

such that equation (3.47) becomes: $$$$m''(x)+(1+x^n)\cdot m(x)=-1 \tag{3.50}$$$$

with the boundary conditions: $$$$m(-1)=m(1)=0,\ m'(0)=0 \tag{3.51}$$$$

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: $$$$u(x)=m(x)-\frac{1}{2}(1-x^2) \tag{3.52}$$$$

The displacement in physical dimensions $$U$$ may subsequently computed from the dimensionless $$u$$: $$$$U=\frac{qL^4}{EI_0} u \tag{3.53}$$$$

Further, the differential equation for $$u(x)$$ may be found by substitution of equation (3.52) into (3.50): $$$$u''(x)+(1+x^n)\cdot u(x)=-\frac{1}{2}(1-x^2)\cdot (1+x^n) \tag{3.54}$$$$

By exploiting the symmetry in the problem, we solve equation (3.50) numerically with a shooting method. The boundary conditions are: $$$$m(1)=0,\qquad m'(0)=0 \tag{3.55}$$$$

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: $$$$y_2(0)=y_1(1)=0 \tag{3.57}$$$$

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)$$.

Simplified shooting method for linear ODEs
1. Let the two initial values $$s^0=0$$ and $$s^1=1$$ for simplicity
2. 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)).
3. Find the correct initial value $$s^*$$ by linear interpolation:
$$$$s^*=\frac{\phi^0}{\phi^0-\phi^1} \tag{3.58}$$$$

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')

Figure 36: Moment $$m$$ and deflection $$u$$ .