November, 2018
We want to find a Ritz approximation of vibration frequencies/modes of a cantilever beam. This is described by the following eigenvalue problem.
$$ \frac{\mathrm{d}^4w}{\mathrm{d}x^4} + \beta^4 w\, ,\quad 0 < x < L,\quad EI>0\, , $$with
$$ w(0) = w'(0) = 0,\quad \left(\frac{\mathrm{d}^2w}{\mathrm{d}x^2}\right)_{x=L} = 0,\quad \left(\frac{\mathrm{d}^3 w}{\mathrm{d}x^3}\right)_{x=L} = 0\, , $$and
$$\beta \equiv \left(\frac{\mu \omega^2}{EI}\right)^{1/4}\, .$$
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
In [2]:
%matplotlib notebook
init_printing()
# Graphics setup
gray = '#757575'
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["text.color"] = gray
plt.rcParams["font.size"] = 12
plt.rcParams["xtick.color"] = gray
plt.rcParams["ytick.color"] = gray
plt.rcParams["axes.labelcolor"] = gray
plt.rcParams["axes.edgecolor"] = gray
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
plt.rcParams["figure.figsize"] = 4, 3
The exact solution for this problem is
$$w_n(x) = A_1\left[(\cosh\beta_n x - \cos\beta_n x) + \frac{\cos\beta_n L + \cosh\beta_n L}{\sin\beta_n L + \sinh\beta_n L}(\sin\beta_n x - \sinh\beta_n x)\right]\, ,$$where $\beta_n$ is the $n$th root of $\cosh(\beta_n L)\cos(\beta_n L) + 1$.
In [3]:
x = symbols('x')
beta, L = symbols("beta L")
val1 = 1.8750909912
val2 = 4.6941049111
In [4]:
def plot_expr(expr, x, rango=(0, 1), ax=None, linestyle="solid"):
"""Plot SymPy expressions of a single variable"""
expr_num = lambdify(x, expr, "numpy")
x0 = rango[0]
x1 = rango[1]
x_num = np.linspace(0, 1, 101)
if ax is None:
plt.figure()
ax = plt.gca()
ax.plot(x_num, expr_num(x_num), linestyle=linestyle)
The quadratic functional for this problem is
$$J[u] = \int\limits_0^L \left[\left(\frac{\mathrm{d}^2 w}{\mathrm{d}x^2}\right)^2 + \beta^4 w^2\right]\mathrm{d}x\, ,$$and the weak problem $B(v, w) = \beta^2 A(v, w)$, with
$$ B(v, u) = \int\limits_0^L \frac{\mathrm{d}^2 v}{\mathrm{d}x^2}\frac{\mathrm{d}^2 u}{\mathrm{d}x^2}\mathrm{d}x\, ,\quad A(v, u) = \int\limits_0^L vw\mathrm{d}x\, . $$
In [5]:
def quad_fun(x, w, L):
F1 = diff(w, x, 2)**2
F2 = w**2
U = integrate(F1, (x, 0, L))
T = integrate(F2, (x, 0, L))
return U, T
In [6]:
def ritz_conventional(x, L, nterms):
a = symbols("a0:%i"%(nterms))
w = sum(a[k]*x**(k + 2) for k in range(nterms))
U, T = quad_fun(x, w, L)
K = Matrix(nterms, nterms, lambda i, j: U.diff(a[i], a[j]))
M = Matrix(nterms, nterms, lambda i, j: T.diff(a[i], a[j]))
return K, M
In [7]:
K, M = ritz_conventional(x, 1, 2)
In [8]:
Kaux = M.inv() * K
vals = list(Kaux.eigenvals())
vals
Out[8]:
In [9]:
nvals = [N(val**0.25) for val in vals]
nvals
Out[9]:
We can write the problem as minimizing the functional
$$J(\psi, w) = \int\limits_0^L\left[\left(\frac{\mathrm{d} \psi}{\mathrm{d}x}\right)^2 + \beta^4 w^2\right]\mathrm{d}x\, ,$$subject to
$$G(\psi, w) \equiv \psi + \frac{\mathrm{d}w}{\mathrm{d}x} = 0\, .$$The Lagrangian is given by
$$L(\psi, w, \lambda) = \int\limits_0^L\left[\left(\frac{\mathrm{d} \psi}{\mathrm{d}x}\right)^2 + \beta^4 w^2\right]\mathrm{d}x + \int\limits_0^L \lambda\left(\psi + \frac{\mathrm{d}w}{\mathrm{d}x}\right)\mathrm{d}x\, , $$where $\lambda$ is the Lagrange multiplier, which in this case represents the shear force.
In [10]:
def lagran(x, w, psi, lamda, L):
F1 = diff(psi, x)**2 + lamda*(psi + diff(w, x))
F2 = w**2
U = integrate(F1, (x, 0, L))
T = integrate(F2, (x, 0, L))
return U, T
In [11]:
def ritz_multiplier(x, L, nterms):
a = symbols("a0:%i"%(nterms))
b = symbols("b0:%i"%(nterms))
c = symbols("c0:%i"%(nterms))
var = a + b + c
w = sum(a[k]*x**(k + 1) for k in range(nterms))
psi = sum(b[k]*x**(k + 1) for k in range(nterms))
lamda = sum(c[k]*x**k for k in range(nterms))
U, T = lagran(x, w, psi, lamda, L)
K = Matrix(3*nterms, 3*nterms, lambda i, j: U.diff(var[i], var[j]))
M = Matrix(3*nterms, 3*nterms, lambda i, j: T.diff(var[i], var[j]))
return K, M
In [12]:
K, M = ritz_multiplier(x, 1, 2)
In [13]:
K
Out[13]:
In [14]:
M
Out[14]:
In [15]:
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [N(1/val**0.25) for val in vals if val != 0]
nvals
Out[15]:
In [16]:
def augmented(x, w, psi, K, L):
F1 = diff(psi, x)**2 + S(K)/2*(psi + diff(w, x))**2
F2 = w**2
U = integrate(F1, (x, 0, L))
T = integrate(F2, (x, 0, L))
return U, T
In [17]:
def ritz_penalty(x, K, L, nterms):
a = symbols("a0:%i"%(nterms))
b = symbols("b0:%i"%(nterms))
var = a + b
w = sum(a[k]*x**(k + 1) for k in range(nterms))
psi = sum(b[k]*x**(k + 1) for k in range(nterms))
U, T = augmented(x, w, psi, K, L)
K = Matrix(2*nterms, 2*nterms, lambda i, j: U.diff(var[i], var[j]))
M = Matrix(2*nterms, 2*nterms, lambda i, j: T.diff(var[i], var[j]))
return K, M
In [18]:
K, M = ritz_penalty(x, 100, 1, 2)
In [19]:
K
Out[19]:
In [20]:
M
Out[20]:
In [21]:
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [re(N(1/val**0.25)) for val in vals if val != 0]
nvals
Out[21]:
The mixed formulation involves rewriting a given higher order equation as a pair of lower order equations by introducing secondary dependent variables. The original equation can be decomposed into
$$ M(x)= \frac{\mathrm{d}^2 w}{\mathrm{d}x^2}\, ,\quad \frac{\mathrm{d}^2M(x)}{\mathrm{d}x^2} = -\beta^4 w\, ,\quad 0<x<L\, . $$The functional in this case is
$$ I(w, M) = \int\limits_0^L\left(\frac{\mathrm{d}w}{\mathrm{d}x}\frac{\mathrm{d}M}{\mathrm{d}x} + \frac{M^2}{2} - \beta^4 w^2\right)\mathrm{d}x $$
In [22]:
def mixed_fun(x, w, M, L):
F1 = diff(w, x)*diff(M, x) + M**2/2
F2 = -w**2
U = integrate(F1, (x, 0, L))
T = integrate(F2, (x, 0, L))
return U, T
In [23]:
def ritz_mixed(x, L, nterms):
a = symbols("a0:%i"%(nterms))
b = symbols("b0:%i"%(nterms))
var = a + b
w = sum(a[k]*x**(k + 1) for k in range(nterms))
M = sum(b[k]*(x - L)**(k + 1) for k in range(nterms))
display(w, M)
U, T = mixed_fun(x, w, M, L)
Kmat = Matrix(2*nterms, 2*nterms, lambda i, j: U.diff(var[i], var[j]))
Mmat = Matrix(2*nterms, 2*nterms, lambda i, j: T.diff(var[i], var[j]))
return Kmat, Mmat
In [24]:
K, M = ritz_mixed(x, 1, 2)
In [25]:
K
Out[25]:
In [26]:
M
Out[26]:
In [27]:
Maux = K.inv() * M
vals = list(Maux.eigenvals())
nvals = [re(N(1/val**0.25)) for val in vals if val != 0]
nvals
Out[27]:
In [ ]: