Finite differences for linear elasticity

Definitions



In [1]:

    
from sympy import *
from continuum_mechanics.vector import lap, sym_grad
from continuum_mechanics.solids import navier_cauchy, strain_stress



In [2]:

    
init_printing()



In [3]:

    
x, y = symbols("x y")
lamda, mu, h = symbols("lamda mu h")



In [4]:

    
def construct_poly(pts, terms, var="u"):
    npts = len(pts)
    u = symbols("{}:{}".format(var, npts))
    vander = Matrix(npts, npts,
                    lambda i, j: (terms[j]).subs({x: pts[i][0],
                                                  y: pts[i][1]}))
    inv_vander = simplify(vander.inv())
    shape_funs = simplify(inv_vander.T * Matrix(terms))
    poly = sum(Matrix(u).T * shape_funs)
    return poly

Nine-point stencil

We can compute the finite difference for elasticity applying the Navier-Cauchy operator to a polynomial interpolator for our stencil.



In [5]:

    
pts = [[0, 0],
       [h, 0],
       [0, h],
       [-h, 0],
       [0, -h],
       [h, h],
       [-h, h],
       [-h, -h],
       [h, -h]]
terms = [S(1), x, y, x**2, x*y, y**2, x**2*y, x*y**2, x**2*y**2]

We can construct the interpolators for horizontal and vertical components of the displacement vector.



In [6]:

    
U = construct_poly(pts, terms, "u")
V = construct_poly(pts, terms, "v")
disp = Matrix([U, V, 0])

Let's take a look at one of the components



In [7]:

    
disp[0]









    Out[7]:





$$\frac{u_{0}}{h^{4}} \left(h^{4} - h^{2} \left(x^{2} + y^{2}\right) + x^{2} y^{2}\right) + \frac{u_{1} x}{2 h^{4}} \left(h^{3} + h^{2} x - h y^{2} - x y^{2}\right) + \frac{u_{2} y}{2 h^{4}} \left(h^{3} + h^{2} y - h x^{2} - x^{2} y\right) + \frac{u_{3} x}{2 h^{4}} \left(- h^{3} + h^{2} x + h y^{2} - x y^{2}\right) + \frac{u_{4} y}{2 h^{4}} \left(- h^{3} + h^{2} y + h x^{2} - x^{2} y\right) + \frac{u_{5} x y}{4 h^{4}} \left(h^{2} + h \left(x + y\right) + x y\right) + \frac{u_{6} x y}{4 h^{4}} \left(- h^{2} + h \left(x - y\right) + x y\right) + \frac{u_{7} x y}{4 h^{4}} \left(h^{2} - h \left(x + y\right) + x y\right) + \frac{u_{8} x y}{4 h^{4}} \left(- h^{2} + h \left(- x + y\right) + x y\right)$$

This expression is quite lengthy to manipulate by hand, but we can obtain the finite difference for the Navier operator straightforward.



In [8]:

    
simplify(navier_cauchy(disp, [lamda, mu]).subs({x:0, y:0}))









    Out[8]:





$$\left[\begin{matrix}\frac{1}{4 h^{2}} \left(\mu \left(- 8 u_{0} + 4 u_{2} + 4 u_{4} - v_{5} + v_{6} - v_{7} + v_{8}\right) + \left(\lambda + 2 \mu\right) \left(- 8 u_{0} + 4 u_{1} + 4 u_{3} + v_{5} - v_{6} + v_{7} - v_{8}\right)\right)\\\frac{1}{4 h^{2}} \left(\mu \left(- u_{5} + u_{6} - u_{7} + u_{8} - 8 v_{0} + 4 v_{1} + 4 v_{3}\right) + \left(\lambda + 2 \mu\right) \left(u_{5} - u_{6} + u_{7} - u_{8} - 8 v_{0} + 4 v_{2} + 4 v_{4}\right)\right)\\0\end{matrix}\right]$$

To impose Neuman boundary conditions we need to compute the stresses.



In [9]:

    
strain = (sym_grad(disp)).subs({x:0, y:0})
stress = strain_stress(strain, [lamda, mu])

The tractions are the projection of the stress tensor. For a uniform grid these would be horizontal or vertical.



In [10]:

    
t1 = stress * Matrix([1, 0, 0])
simplify(2*h*t1)









    Out[10]:





$$\left[\begin{matrix}\lambda \left(u_{1} - u_{3} + v_{2} - v_{4}\right) + 2 \mu \left(u_{1} - u_{3}\right)\\\mu \left(u_{2} - u_{4} + v_{1} - v_{3}\right)\\0\end{matrix}\right]$$



In [11]:

    
t2 = stress * Matrix([0, 1, 0])
simplify(2*h*t2)









    Out[11]:





$$\left[\begin{matrix}\mu \left(u_{2} - u_{4} + v_{1} - v_{3}\right)\\\lambda \left(u_{1} - u_{3} + v_{2} - v_{4}\right) + 2 \mu \left(v_{2} - v_{4}\right)\\0\end{matrix}\right]$$



In [ ]:



In [12]:

    
from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()









    Out[12]:



In [ ]: