import sympy
import numpy
from matplotlib import pyplot
from sympy import symbols
from IPython.display import display, Latex
sympy.init_printing(use_latex=True)
from sympy.utilities.lambdify import lambdify
class Kinematics(object):
def __init__(self, var = symbols('x'), l = 1.0, dimension = 2):
self.var, self.l, self.dimension = var, l, dimension
def basis_functions(self):
var = self.var
return [1-3*var**2+2*var**3, var-2*var**2+var**3, 3*var**2-2*var**3, -var**2+var**3]
def phi_matrix(self):
phi, l, eye = self.basis_functions(), self.l, sympy.eye(self.dimension)
return sympy.Matrix([phi[0] * eye, l * phi[1] * eye, phi[2] * eye, l * phi[3] * eye]).T
x = symbols('x')
kin = Kinematics(var = x)
S = kin.phi_matrix()
e contains the nodal coordinates for a single element r is the product of the shape functions in S and the nodal coordinates e
e = sympy.Matrix(sympy.symbols(['r[%d]'%i for i in range(8)]))
r = S * e
rx = r.diff(x)
rxx = r.diff(x,2)
Latex('$ \\small e^{\\rm T} = '+sympy.latex(e.T)
+'\\\\ \\small s = ' + sympy.latex(S)
+'\\\\ \\small r = S e = ' + sympy.latex(r)
+'\\\\ \\frac{\partial r}{\partial x} = ' + sympy.latex(rx)
+'\\\\ \\frac{\partial^2 r}{\partial^2 x} = ' + sympy.latex(rxx)
+ '$')
elements = [0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0, 1.0, 0.6, 1.0, 0.0], [0.0, 0.0, 1.0, -1.0, 1.0, 0.0, 1.0, -1.0]
lam_f = lambdify([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], x], r)
def eval_element(nodal_coordinates):
nc = nodal_coordinates
x_coords, y_coords = [], []
for t in numpy.linspace(0, 1, 50):
xc, yc = lam_f(nc[0], nc[1], nc[2], nc[3], nc[4], nc[5], nc[6], nc[7], t)
x_coords.append(xc)
y_coords.append(yc)
return x_coords, y_coords
coords = [eval_element(element) for element in elements]
%matplotlib inline
pyplot.figure(figsize=(17,3))
pyplot.subplot(1,3,1)
pyplot.plot(coords[0][0], coords[0][1])
pyplot.xlim([-0.1,1.2])
pyplot.ylim([-1,1])
pyplot.title('First a boring straight element from (0,0) to (1, 0)')
pyplot.subplot(1,3,2)
pyplot.plot(coords[1][0], coords[1][1])
pyplot.xlim([-0.1,1.2])
pyplot.ylim([-1,1])
pyplot.title('A beam element from (0,0) to (1, 0.6)')
pyplot.subplot(1,3,3)
pyplot.plot(coords[2][0], coords[2][1])
pyplot.xlim([-0.1,1.2])
pyplot.ylim([-1,1])
pyplot.title('A beam element from (0,0) to (1, 0) with some sway');
class Dynamics(object):
def __init__(self):
self.g = -9.81
self.l = 1.0 # Length of beam element
# The following are arbitrary values
self.rho = 1000.0 # Density (mass per volume)
self.A = 0.4 # Area of cross section
self.I = 0.03 # Inertia of cross section
self.E = 1000.0 # Youngs modulus
def gravitational_force(self):
# Taken from Equation 18 in "Application of the Absolute Nodal Coordinate Formulation to multibody system dynamics"
# by "J.L. Escalona, H.A. Hussein and A.A. Shabana, May 1997
return self.rho * self.A * sympy.integrate(self.g * S[1,:], (x, 0, self.l))
def mass_matrix(self):
# Taken from Equation 8 in "Application of the Absolute Nodal Coordinate Formulation to multibody system dynamics"
# by "J.L. Escalona, H.A. Hussein and A.A. Shabana, May 1997
return sympy.integrate(self.rho * S.T * S, (x, 0, self.l))
def strain_energy_longitudinal_deformation(self):
# Strain energy based on the continuum mechanics formulation in
# Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation"
# by "Berzeri, M and Shabana, Ahmed A", 2000
# Strain energy due to longitudinal deformation
return 0.5 * self.E * self.A * sympy.integrate(((rx.T * rx)[0] - 1)**2, (x, 0, self.l))
def strain_force_longitudinal_deformation(self):
# Strain energy based on the continuum mechanics formulation in
# Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation"
# by "Berzeri, M and Shabana, Ahmed A", 2000
# Strain energy due to longitudinal deformation
seld = self.strain_energy_longitudinal_deformation()
return sympy.Matrix([seld.diff(ee) for ee in e])
def strain_energy_curvature_function(self):
# Strain energy based on the continuum mechanics formulation in
# Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation"
# by "Berzeri, M and Shabana, Ahmed A", 2000
# Strain energy due to curvature
# Values for the Gauss-Legendre quadrature
# https://en.wikipedia.org/wiki/Gaussian_quadrature
x_vals = list(map(lambda val: self.l*(val+1.0)/2.0, [0.0,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585]))
weights = [0.4179591836734694, 0.3818300505051189, 0.3818300505051189, 0.2797053914892766, 0.2797053914892766, 0.1294849661688697, 0.1294849661688697]
def cross_prod_2d(a, b):
return sympy.Matrix([a[0] * b[1] - a[1] * b[0]])
num = (cross_prod_2d(rx, rxx).T * cross_prod_2d(rx, rxx))[0,0]
den = (rx.T*rx * rx.T*rx * rx.T*rx)[0,0]
expr = self.E * self.I * num / den
def integrate_expr(nodal_coordinates):
nc_subs = list(zip(e, nodal_coordinates))
expr_ = expr.subs(nc_subs)
res = 0.0
for x_val, weight in zip(x_vals, weights):
res += weight * expr_.subs({x : x_val})
return res
return integrate_expr
def strain_force_curvature_function(self):
# Strain energy based on the continuum mechanics formulation in
# Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation"
# by "Berzeri, M and Shabana, Ahmed A", 2000
# The force of strain due to curvature
# Values for the Gauss-Legendre quadrature
# https://en.wikipedia.org/wiki/Gaussian_quadrature
x_vals = list(map(lambda val: self.l*(val+1.0)/2.0, [0.0,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585]))
weights = [0.4179591836734694, 0.3818300505051189, 0.3818300505051189, 0.2797053914892766, 0.2797053914892766, 0.1294849661688697, 0.1294849661688697]
def cross_prod_2d(a, b):
return sympy.Matrix([a[0] * b[1] - a[1] * b[0]])
num = (cross_prod_2d(rx, rxx).T * cross_prod_2d(rx, rxx))[0,0]
den = (rx.T*rx * rx.T*rx * rx.T*rx)[0,0]
expr = self.E * self.I * num / den
dexpr_de = sympy.Matrix([expr.diff(ee) for ee in e])
def integrate_expr(nodal_coordinates):
nc_subs = list(zip(e, nodal_coordinates))
dexpr_de_ = dexpr_de.subs(nc_subs)
res = sympy.zeros(8,1)
for x_val, weight in zip(x_vals, weights):
res += weight * dexpr_de_.subs({x : x_val})
return res
return integrate_expr
dyn = Dynamics()
dyn.gravitational_force()
Latex('$\\tiny'+sympy.latex(dyn.mass_matrix())+'$')
lam_f = lambdify(e, dyn.strain_energy_longitudinal_deformation())
def element_strain_energy_longitudinal_deformation(nodal_coordinates):
nc = nodal_coordinates
return lam_f(nc[0], nc[1], nc[2], nc[3], nc[4], nc[5], nc[6], nc[7])
fnc1 = dyn.strain_energy_curvature_function()
fnc2 = element_strain_energy_longitudinal_deformation
sympy.Matrix([[fnc1(element), fnc2(element)] for element in elements])
lam_f = lambdify(e, dyn.strain_force_longitudinal_deformation())
def element_strain_force_longitudinal_deformation(nodal_coordinates):
nc = nodal_coordinates
return lam_f(nc[0], nc[1], nc[2], nc[3], nc[4], nc[5], nc[6], nc[7])
fnc1 = dyn.strain_force_curvature_function()
fnc2 = element_strain_force_longitudinal_deformation
res = sympy.zeros(8,6)
for i in range(3):
res[:,i] = fnc1(elements[i])
res[:,i+3] = fnc2(elements[i])
Latex('$\\small'+sympy.latex(res)+'$')