Finite Difference Method

November 26, 2015 - Engineering

This document uses the open source Jupyter Notebook and the Python programming language to demonstrate a finite difference numerical method for solving a differential equation.

%matplotlib inline

# Required libraries
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Seaborn aesthetics
sns.set_context("paper", font_scale=1.8, 
                rc={"lines.linewidth": 3, "figure.figsize": (10,6)})
sns.set_style("ticks")

Beam Deflection Problem

Loaded Beam
Figure 1. Uniformly-loaded beam (Chapra, p. 639).

The basic differential equation of the elastic curve for a simply supported, uniformly loaded beam (Fig. 1) is given as

\begin{equation} EI \frac{d^2y}{dx^2} = \frac{wLx}{2} - \frac{wx^2}{2} \tag{3} \end{equation}

where $E=$ the modulus of elasticity, and $I=$ the moment of inertia. The boundary conditions are $y(0)=y(L)=0$. Solve for the deflection of the beam using the finite-difference approach $(\Delta x = 0.6 m)$. The following parameter values apply: $E = 200$ GPa, $I = 30,000$ cm$^4$, $w = 15$ kN/m, and $L = 3$ m. Compare your numerical results to the analytical solution:

\begin{equation} y = \frac{wLx^3}{12EI} - \frac{wx^4}{24EI} - \frac{wL^3x}{24EI} \tag{4} \end{equation}

Analytical Solution

def y(x):
    """
    Vertical deflection, y, of the beam in 10E-5 m.
    """
    E = 200 * 10**9     # Pa
    I = 30000 / 100**4  # m^4
    w = 15 * 1000       # N/m
    L = 3               # m
    return ((w*L*x**3)/(12*E*I) - (w*x**4)/(24*E*I) - \
           (w*x*L**3)/(24*E*I)) * 10**5

Central Difference Method

def cdiff_beam(L,ya,yb,n):
    """
    Approximate vertical deflection, y, of the beam in 10E-5 m,
    using central difference method.
    Args:
        L  = length of the beam, m
        ya = vertical deflection at x = 0
        yb = vertical deflection at x = L
        n  = desired number of interior segments
    
    """
    E = 200 * 10**9     # Pa
    I = 30000 / 100**4  # m^4
    w = 15 * 1000       # N/m
    i = 0
    x = 0
    dx = L/n
    X = [x]
    A = np.zeros(shape=(n+1,n+1))
    b = np.zeros(shape=(n+1,1))
    while i < n-1:
        i = i + 1
        x = x + dx
        b[i] = dx**2 * (1/(2*E*I)) * (w*L*x - w*x**2)
        A[i][i] = -2
        if i < n:
            A[i][i+1] = 1
        A[i][i-1] = 1
        X = np.append(X, x)
    A[0][0] = 1
    A[n][n] = 1
    b[0] = ya
    b[n] = yb
    X = np.append(X, x+dx)
    Y = np.linalg.solve(A,b)
    # print(A); print(b)
    return X, Y

Plot of Solutions

x = np.linspace(0,3)
X,Y = cdiff_beam(3,0,0,8)

plt.plot(x,y(x), lw=7, alpha=0.6, label="Analytical Solution")
plt.plot(X,Y*10**5, "o-", lw=2, ms=8, color="#c44e52", label="Finite Difference Method")
plt.legend(bbox_to_anchor=[0.665, 0.81])
plt.xlabel("distance, $x$ (m)")
plt.ylabel("deflection, $y$ (10$^{-5}$ m)")
sns.despine()
plt.show()
Loaded Beam

Comments

comments powered by Disqus