How to solve Euler–Bernoulli beam equation in FiPy?

Question

To understand how FiPy works, I want to solve the Euler–Bernoulli beam equation with fixed endpoints:

w''''(x) = q(x,t),    w(0) = w(1) = 0,  w'(0) = w'(1) = 0.

For simplicity, let q(x,t) = sin(x).

How can I define and solve it in FiPy? How to specify the source term sin(x) with respect to the only independent variable in the equation?

from fipy import CellVariable, Grid1D, DiffusionTerm, ExplicitDiffusionTerm
from fipy.tools import numerix

nx = 50
dx = 1/nx
mesh = Grid1D(nx=nx, dx=dx)

w = CellVariable(name="deformation",mesh=mesh,value=0.0)

valueLeft = 0.0
valueRight = 0.0

w.constrain(valueLeft, mesh.facesLeft)
w.constrain(valueRight, mesh.facesRight)
w.faceGrad.constrain(valueLeft, mesh.facesLeft)
w.faceGrad.constrain(valueRight, mesh.facesRight)

# does not work:
eqX = DiffusionTerm((1.0, 1.0)) == numerix.sin(x)
eqX.solve(var=w)

Answer 1

Here is what seems to be a working version of your problem

from fipy import CellVariable, Grid1D, DiffusionTerm
from fipy.tools import numerix
from fipy.solvers.pysparse.linearPCGSolver import LinearPCGSolver
from fipy import Viewer
import numpy as np


L = 1.
nx = 500
dx = L / nx
mesh = Grid1D(nx=nx, dx=dx)

w = CellVariable(name="deformation",mesh=mesh,value=0.0)

valueLeft = 0.0
valueRight = 0.0

w.constrain(valueLeft, mesh.facesLeft)
w.constrain(valueRight, mesh.facesRight)
w.faceGrad.constrain(valueLeft, mesh.facesLeft)
w.faceGrad.constrain(valueRight, mesh.facesRight)

x = mesh.x
k_0 = 0
k_1 = -1
k_2 = 2 + np.cos(L) - 3 * np.sin(L)
k_3 = -1 + 2 * np.sin(L) - np.cos(L)
w_analytical = numerix.sin(x) + k_3 * x**3 + k_2 * x**2 + k_1 * x + k_0
w_analytical.name = 'analytical'

# does not work:
eqX = DiffusionTerm((1.0, 1.0)) == numerix.sin(x)
eqX.solve(var=w, solver=LinearPCGSolver(iterations=20000))
Viewer([w_analytical, w]).plot()
raw_input('stopped')

After running this, the FiPy solution seems to be quite close to the analytical result.

The two important changes from the OP's implementation.

• Using mesh.x which is the correct way to refer to the spatial variable for use in FiPy equations.

• Specifying the solver and number of iterations. The problem seems to be slow to converge so needed a lot of iterations. From my experience, fourth order spatial equations often need good preconditioners to converge quickly. You might try using the Trilinos solver package with Fipy to make this work better as it has a wider range of available preconditioners.

• Use an explicit float for L to avoid integer maths in Python 2.7 (edit from comment)