Consider a beam of length L with a straigth reference configuration and denote by $s$ the arc-length in this configuration. Without loss of generality, let us use a Caresian reference system $(O,\underline{e}_1,\underline{e}_2)$ with the origin at the left end of the beam and $\underline{e}_1$ oriented as the beam in the reference configuration. Hence the reference configuration can be wirtten with the following parametric representation
$$ \underline x_R(s) = s \,\underline{e_1} $$The current cunfiguration is given by
$$ \underline x(s) = (s+u(s)) \underline{e}_1 + v(s) \underline{e}_2 $$where $\underline u(s)= u(s) \underline{e}_1 + v(s) \underline{e}_2$ is the displacement vector, with the axial ($u$) and transversal ($v$) compontents wrt the reference configuration.
We can use two python packages, numpy and matplotlib to represent the beam.
First let us define of the reference configuration as two functions:
In [2]:
import numpy as np
def xR(s):
return s
def yR(s):
return 0.*s
s = np.linspace(0.,1.,20)
xR(s)
Out[2]:
We can represent it using matplotlib
In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(xR(s),yR(s))
Out[3]:
Then we can do it better with the current configuration
In [4]:
x = lambda s : xR(s) + .1*s
y = lambda s : yR(s) + .3*s**2
with plt.xkcd():
fig, ax = plt.subplots()
ax.plot(x(s),y(s),"r",color='b', ls='-', lw=3)
ax.plot(xR(s),yR(s),"r",color='b', ls='--', lw=3)
ax.quiver(xR(.7),yR(.7),x(.7)-xR(.7),y(.7)-yR(.7),color='r',scale_units='xy', angles='xy', scale=1)
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
ax.set_title("Reference and deformed configuration")
ax.set_xlim(-0.2, 1.2)
ax.set_ylim(-0.05, .35)
ax.text(xR(.7)-.2,yR(.7)-.03,r'$({x}_R(s),y_R(s))$')
ax.text(x(.7)-.3,y(.7)+.01,r'$({x}(s),y(s))$')
The unit tangent venctor to the current configuration is given by
\begin{equation} \underline{t}(s) = \dfrac{\underline{x}'(s)}{\Vert{\underline{x}'(s)}\Vert} = \cos\theta(s)\,\underline{e}_1 + \sin\theta(s) \,\underline{e}_2 \end{equation}where in the last equality we introduced the angle $\theta(s)$ between the tangent $\underline{t}(s)$ and $\underline{e}_1$.
We define the following deformations measures for the beam:
\begin{equation} \begin{cases} \epsilon(s) = \Vert{\underline{x}'(s)}\Vert - 1 = \sqrt{(1+u'(s))^2+v'(s)^2}-1&(extension)\\ \kappa(s) = \theta'(s) & (curvature) \end{cases} \end{equation}When $\epsilon(s)=0$ the beam is inextensible and $s$ is a arc-length also in the current configurations. Otherwise the arclength in the current configuration is given by $\tilde s$, such that $\mathrm d \tilde s ^2= \mathrm d x^2+ \mathrm d y^2 = \Vert{\underline{x}'(s)}\Vert^2 ds^2 $.
The definions of the extension and the tangent give the following expressions for the derivatives of the displacement field:
$$ \begin{cases} u'(s)&=(1+\epsilon(s))\cos\theta(s) - 1\\ v'(s)&=(1+\epsilon(s))\sin\theta(s) \end{cases} $$Let us assume that $\theta$ and $e$ are small. Then by a Taylor expansion
$$ \begin{cases} u'(s)&=(1+\epsilon(s))\cos\theta(s) - 1\simeq \epsilon(s)-\dfrac{\theta(s)^2}{2}\\ v'(s)&=(1+\epsilon(s))\sin\theta(s) \simeq \theta(s)+o(\theta^2) \end{cases} $$Retaing only the linear terms in the Taylor expansion above one get the fully linearized kinematics for a planar Euler-Bernoulli beam:
$$ \begin{cases} \epsilon(s) = u'(s)\\ \theta(s) = v'(s)\quad\Rightarrow\quad \kappa(s)=v''(s) \end{cases} $$We present below the linear model, assuming the fully linearized kinematics above. This model is valid for small loadings. Let us consider the case of beam of length $L$, bending stiffness $EI$ and axial stiffness $ES$, which submitted to a distributed force per unit length $\underline f(s)= f_u(s)\underline e_1 + f_v(s)\underline e_2$.
Let us consider the case of beam of length $L$, bending stiffness $EI$ and axial stiffness $ES$
We assume a linear constitutive model and that the strain energy density (per unit line) is quadratic in the deformation
$$ W(e,k) = \dfrac{EI}{2} \kappa^2 + \dfrac{ES}{2} \epsilon^2 $$Hence the total elastic energy of the beam is given by the following functional
$$ \mathcal{E}(u,v)=\int_0^L\dfrac{EI}{2} v''(s)^2 + \dfrac{ES}{2} u'(s)^2\,ds $$In the case fo conservative loadings, the potential energy, say $\mathcal{P}(u,v)$, is the difference between the elastic energy and the work of the external forces, say $\mathcal{F}(u,v)$.
$$ \mathcal{P}(u,v) = \mathcal{E}(u,v) - \mathcal{F}(u,v) $$Let us consider the case in which the end of the beam are clamped and the external forces are in the form of distributed force per unit length $\underline f(s)= f_u(s)\underline e_1 + f_v(s)\underline e_2 $. Hence the work of external forces is
$$ \mathcal{F}(u,v) = \int_0^L f_u(s) \,u(s)+f_v(s) \,v(s)\,ds $$If concentrated forces and/or coupled are present, their work contribution should be added to the functional above.
Using the principle of the minimum of the potential energy, the equilibrium configurations are found looking for the stationarity point of the energy functional $\mathcal{P}(u,v)$ among all admissible displacement field respecting the boundary conditions.
In a first step, we will suppose that displacements and rotations are imposed at both ends, with
$$ \begin{cases} u(0)= u_0, \, v(0)=v_0, \, v'(0)=\theta_0,\\ u(L)= u_L, \, v(0)=v_L, \, v'(0)=\theta_L \end{cases} $$Instead of solving this problem exactly (see the notes of 3A103), we adopt an approximate method, basis on the discretisation using a Galerking approximation
We look for approximate solutions in the form of linear polynomials for the axial displacement and cubic polynomial for the transverse displacement:
$$ \begin{cases} u(s) = c_0 +c_1 s\\ v(s) = c_2 + c_3 s + c_4 s^2+ c_5 s^3 \end{cases} $$They can be defined in python using sympy. We declare the variable $s$ and the coefficients $c_i$'s as sympy
symbols.
In [5]:
# I use sympy http://docs.sympy.org/
import sympy as sp
from sympy.interactive import printing
printing.init_printing()
s = sp.Symbol('s')
c0, c1, c2, c3, c4, c5 = sp.symbols('c0, c1, c2, c3, c4, c5')
u = c0 + c1*s
v = c2 + c3*s + c4*s**2 + c5*s**3
sp.Matrix([u, v])
Out[5]:
We denote by $q(s)=v'(s)$ the rotation field. To calculate it we can use the symbolic differntiation in sympy
In [6]:
q = v.diff(s)
q
Out[6]:
Let us now define the vector collecting the displacement and rotation at lefth and right ends of the beam
In [7]:
u0, v0, q0, u1, v1, q1 = sp.symbols('u0, v0, q0, u1, v1, q1')
U = sp.Matrix([u0, v0, q0, u1, v1, q1])
U
Out[7]:
The polynomial approximation allows us to express the displacement in any point in the of the end displacement $U$. To do it we have to resolve a linear system. We do it by using the capability of sympy.
In [8]:
L = sp.Symbol('L')
solu = sp.solve([u.subs({s:0})-u0, u.subs({s:L})-u1], (c0,c1))
solv = sp.solve([v.subs({s:0})-v0, q.subs({s:0})-q0, v.subs({s:L})-v1, q.subs({s:L})-q1], (c2,c3,c4,c5))
u = u.subs(solu)
v = v.subs(solv)
sp.Matrix([u, v])
Out[8]:
We can then define the following shape functions $S_{ij}$ that gives the displacement field in the beam as a function of the displacement and rotations at the ends: $$ \begin{bmatrix} u(s)\ v(s)
S(s)U $$ We can calculate the matrix S by using the derivative to extract the coefficients with respect to the components of $U$
In [9]:
S = sp.Matrix([[u.diff(Ui) for Ui in U],[v.diff(Ui) for Ui in U]])
S
Out[9]:
We used above what is called a python list comprehension, which is a syntax to write table, for example
In [10]:
[i**2 for i in [3,4,5]]
Out[10]:
We can represent the basis functions using matplotlib
In [11]:
sn = np.linspace(0.,1.,20)
for i in range(2):
Sn = np.array([S[0,:].subs({s:si, L:1}) for si in sn])
plt.plot(sn,Sn)
plt.xlabel('s')
plt.ylabel('shape functions')
Out[11]:
In [11]:
for i in range(6):
Sn = np.array([S[1,:].subs({s:si, L:1}) for si in sn])
plt.plot(sn,Sn)
plt.xlabel('s')
plt.ylabel('shape functions')
In [12]:
eps = u.diff(s)
kappa = v.diff(s,2)
[eps,kappa]
Out[12]:
In [14]:
B = sp.Matrix([[eps.diff(Ui) for Ui in U],[kappa.diff(Ui) for Ui in U]])
B
In [14]:
ES = sp.Symbol('ES')
Ke_ext = ES*sp.integrate(B[0,:].transpose()*B[0,:],(s,0,L))
Ke_ext
Out[14]:
Similarly for the bending energy $$ \mathcal E_{\mathrm{bend}} = \int_0^L\dfrac{EI}{2} v''(s)^2 = \dfrac{1}{2}K^e_{\mathrm{ben}} U\cdot U $$ where
In [15]:
EI = sp.Symbol('EI')
Ke_bend = EI*sp.integrate(B[1,:].transpose()*B[1,:],(s,0,L))
Ke_bend
Out[15]:
Hence the total elastic energy of the beam is given by $$ \mathcal E = \dfrac{1}{2}K_e U\cdot U $$ where $K_e$ is called stiffness matrix
In [16]:
Ke = Ke_ext + Ke_bend
Ke
Out[16]:
Using the polynomial approximation for the displacement the work of the external forces gives $$ \mathcal{F} = \int_0^L f_u(s) \,u(s)+f_v(s) \,v(s)\,ds = F^{(e)}\cdot U $$ where $$ F^{(e)}_i = \int_0^L f_u(s) \,S_{0i}(s)+f_v(s) \,S_{1i}(s)\,ds $$ $S$ Being the matrix of shape function define above.
For example, for the case of a beam of mass density per unit line $\rho S$ submitted its own weight
In [17]:
p = sp.Symbol("p")
f_u = 0
f_v = - p
where $p=\rho S\,g$ and
In [18]:
Fe = sp.Matrix([sp.integrate(f_u *S[0,i]+f_v*S[1,i],(s,0,L)) for i in range(6)])
Fe
Out[18]:
In [19]:
import numpy as np
Ken = np.array(Ke.subs({EI:1., ES:100., L:1.}))
Fen = np.array(Fe.subs({L:1., p:1.}))
Ken,Fen
Out[19]:
Or, if we want to print it out nicely
In [20]:
sp.Matrix(Ken),sp.Matrix(Fen)
Out[20]:
We can define nodes as points
In [21]:
n_elements = 4
xnodes = np.linspace(0,1,n_elements + 1)
ynodes = np.linspace(0,0,n_elements + 1)
nodes = np.array([xnodes,ynodes]).transpose()
n_nodes = xnodes.size
nodes
Out[21]:
And elements as lists of nodes
In [22]:
elements = np.array([[i,i+1] for i in range(n_elements)])
elements
Out[22]:
We show below how to plot the element and the nodes in python using matplotlib (http://matplotlib.org)
In [23]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.plot(xnodes,ynodes,'o-')
Out[23]:
Or you can have fun and make it prettier
In [24]:
with plt.xkcd():
fig, ax = plt.subplots()
ax.plot(xnodes,ynodes,'o-',lw=2, color='black', ms=10)
for i in range(n_nodes):
ax.text(xnodes[i], -.02, np.array([i,i+1,i+2]))
ax.text(xnodes[i], .01, i, color = 'b')
for e in range(n_elements):
ax.text(xnodes[e]+(xnodes[e+1]-xnodes[e])/2., .02, str(e), bbox=dict(facecolor='red'))
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
ax.set_title("Elements, nodes, and dofs")
ax.set_xlim(-0.2, 1.2)
ax.set_ylim(-0.05, 0.05)
Let us discretise each beam segment with the same polynomial function as above. We collect the displacement and rotations of all the nodes in the vector of the global degrees of freedom $$ U= \begin{bmatrix} u_1&v_1&q_1&u_2&v_2&q_2,\ldots,u_n&v_n&q_n \end{bmatrix}^T $$
For a generic element we can express the displacement in each element as the displacement and rotation of its end nodes $$ U^{(e)}= \begin{bmatrix} \tilde u_0&\tilde v_0& \tilde q_0& \tilde u_1& \tilde v_1& \tilde q_1 \end{bmatrix}^T $$ which are the local degrees of freedom. We look how to define the correspondance between the components of the two vectors.
We call dofmap the function associating to an element number $e$ and the local index of a component of $U^{(e)}$, say $i_\textrm{local}$, the global index of the corresponding component of $U$, say $i$: $$ i = \textrm{dofmap}(e,i_\textrm{local}) $$ such that $$ U_i = U^{(e)}_{i_\textrm{local}} $$ For the case of the beam above, one can easily check that $ i = 3*e + i_\textrm{local} $
In [25]:
def dof_map(e, i_local):
return 3*(e)+i_local
For example, the global indeces of the dof of the element $e=2$ are (python indeces start from $0$)
In [26]:
def dof_map_e(e):
return [dof_map(e, i_local) for i_local in range(6)]
dof_map_e(2)
Out[26]:
We can write the elastic energy of the beam in the following form
$$ \mathcal E = \sum_{el=1}\int_{el}\dfrac{ES}{2}e^2+\dfrac{EI}{2}k^2ds=\sum_{el=1}\dfrac{1}{2}K_{e} U^{(e)}\cdot U^{(e)} $$where $K^{e}$ is the local stifness matrix of the element $e$
We define the global stiffness matrix $K$ the matrix such that $$ \mathcal E = \sum_{el=1}\dfrac{1}{2}K^{(e)} U^{(e)}\cdot U^{(e)} = \dfrac{1}{2}K U\cdot U $$
Here $K^{(e)}$ is $6\times 6$, whilst $K$ is $3 n\times 3 n$.
The global stiffness is obained as follows
In [27]:
sp.Matrix(Ken)
Out[27]:
In [28]:
K = np.zeros([3*n_nodes,3*n_nodes])
for e in range(n_elements):
for i_local in range(6):
for j_local in range(6):
K[dof_map(e, i_local),dof_map(e, j_local)] += Ken[i_local,j_local]
sp.Matrix(K)
Out[28]:
We can also represent it graphically
In [29]:
plt.matshow(K, interpolation='none')
plt.colorbar()
Out[29]:
As for the elastic energy, we can calculate the work of external forces using a global force vector. $$ \mathcal{F} = \int_0^L f_u(s) \,u(s)+f_v(s) \,v(s)\,ds = \sum_e F^{(e)}\cdot U^{(e)}= F\cdot U $$ where the global vector is found by assembling the local ones.
For the case of the beam submitted to its own weight
In [30]:
F = np.zeros([3*n_nodes])
for e in range(n_elements):
for i_local in range(6):
F[dof_map(e, i_local)] += Fen[i_local]
sp.Matrix(F)
Out[30]:
In our case of a clamped 4-element beam, the first three and the last three component of $U$ are blocked.
In [31]:
blocked_dof = np.array([0, 1, 2, F.size-3, F.size-2, F.size-1])
blocked_dof
Out[31]:
Their values are all imposeed to zero. We collect these values in a vector
In [32]:
bc_values = np.array([0, 0, 0, 0, 0, 0])
bc_values
Out[32]:
It is useful also to define the following vector, saying if a giving dof is blocked or not
In [33]:
bc_flags = 0*F.astype("bool")
bc_flags[blocked_dof] = 1
bc_flags
Out[33]:
In [34]:
def bc_apply(K,F,blocked_dof,bc_values):
for (i, dof) in enumerate(blocked_dof):
Kbc = K
Fbc = F
Kbc[dof, :] = 0
Kbc[:, dof] = 0
Kbc[dof, dof] = 1
Fbc += - K[:,dof]*bc_values[i]
Fbc[dof] = bc_values[i]
return Kbc, Fbc
In [35]:
Kbc, Fbc = bc_apply(K, F, blocked_dof, bc_values)
sp.Matrix(Kbc), sp.Matrix(Fbc)
Out[35]:
To solve the linear systems we use the tools available in numpy (we could also program our LU or CG solver ...)
In [36]:
# from scipy.sparse import linalg
Usol = np.linalg.solve(Kbc,Fbc)
Usol
Out[36]:
This is a simple visualization of the solution
In [37]:
def plot_deformed(Usol):
x_c = nodes[:,0]+[Usol[3*e] for e in range(n_nodes)]
y_c = nodes[:,1]+[Usol[3*e+1] for e in range(n_nodes)]
return plt.plot(x_c, y_c, '-o')
plot_deformed(Usol)
Out[37]:
To impose th bcs, we solve a reduced system, where we eliminate the rows and columns corresponding to blocked dof. We call K_red and F_red the new matrix
In [38]:
K_red = np.delete(K, blocked_dof, axis = (0))
K_red = np.delete(K_red, blocked_dof, axis = (1))
F_red = np.delete(F, blocked_dof, axis = (0))
sp.Matrix(K_red),sp.Matrix(F_red)
Out[38]:
We solve the reduced linear system using the tools available in numpy (more on the next classes)
In [39]:
from scipy.sparse import linalg
Usol_red = np.linalg.solve(K_red,F_red)
Usol_red
Out[39]:
Of course we miss the points with blocked dofs ... But we can easily reintegrate the missing point into the vector
In [40]:
Usol_full = np.zeros(F.size)
Usol_full[bc_flags==0]=Usol_red
print Usol_full
plot_deformed(Usol_full)
Out[40]: