notebook.community

Edit and run

Airy stress functions method

The method of Airy stress functions is a method to find the stress distribution in an elastic solid. The method uses a Airy stress function $\varphi$, a scalar potential, and can be used when:

  • The problem is two dimensional.
  • The body forces are derivable from a potential function $V$, i.e.,
$$\mathbf{b} = -\nabla V\, .$$

Computation of stresses

The stresses are computed from the Airy functions as

\begin{align} \sigma_{xx} &= \frac{\partial^2 \varphi}{\partial y^2}\\ \sigma_{yy} &= \frac{\partial^2 \varphi}{\partial x^2}\\ \sigma_{xy} &= -\frac{\partial^2 \varphi}{\partial x \partial y}\, , \end{align}

in Cartesian coordinates, and

\begin{align} \sigma_{rr} &= \frac{1}{r}\frac{\partial \varphi}{\partial r} + \frac{1}{r^2}\frac{\partial^2 \varphi}{\partial r^2}\\ \sigma_{\theta\theta} &= \frac{\partial^2 \varphi}{\partial r^2}\\ \sigma_{r\theta} &= -\frac{\partial}{\partial r}\left(\frac{1}{r}\frac{\partial \varphi}{\partial \theta}\right)\, , \end{align}

in polar coordinates.

Airy stress functions are biharmonic

In 2D, the compatibility equations and equilibrium imply

$$\nabla^2 \sigma_v = \frac{4}{1 + \kappa} \nabla^2 V\, ,$$

where $\sigma_v = \mathrm{Tr}\,(\sigma)$ is the volumetric component of the stress tensor, and

$$\kappa = \begin{cases} 3 - 4\nu\quad \text{plane strain}\\ \frac{3 - \nu}{1 + \nu}\quad \text{plane strain} \end{cases}$$

Then, plugging the Airy stress function into this equation leads to

$$\nabla^2 \nabla^2 \varphi = -\frac{2(\kappa - 1)}{(1 + \kappa)} \nabla^2 V$$

or, just

$$\nabla^4 \varphi = -\frac{2(\kappa - 1)}{(1 + \kappa)} \nabla^2 V$$

for short.

This means that the Airy stress functions are biharmonic, that is, they satisfy the above equation.





In [1]:



    


from sympy import *
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams



    











In [2]:



    


x, y, z, t = symbols('x y z t')
f, g, h = symbols('f g h', cls=Function)
init_printing()



    











In [3]:



    


%matplotlib notebook
plt.style.use("seaborn-notebook")
plt.rcParams["image.cmap"] = "RdYlBu_r"

The following cell contains the routines that we will use in the remaining part of the notebook.





In [4]:



    


def poly_airy(x, y, order, coeff="c"):
    """Return a polynomial of a given order in x and y
    
    The coefficients of the polynomial are returned as a set in the
    variable ``c``.
    """
    num = (order**2 + 3*order + 2)/2
    c = symbols('%s0:%d' % (coeff, num))
    terms = [x**(order - k - j) * y**k
             for j in range(0, order - 1)
             for k in range(0, order + 1 - j)]
    phi = sum(c[k]*term for k, term in enumerate(terms))
    return phi, c


def airy_stress(phi, x, y):
    """Return stresses from the Airy stress function
    
    The coordinate system is Cartesian
    """
    return diff(phi, y, 2), diff(phi, x, 2), -diff(phi, y, 1, x, 1)


def conds2eqs(conds, var=(x, y)):
    """Turn the list of boundary conditions into equations to be solved
    
    In general, the final list can have repeated equations. Thus, we
    are removing those before returning it.
    """
    eqlist = [coef for cond in conds for coef in Poly(cond, var).coeffs()]
    return list(set(eqlist))


def plot_contours(x_coord, y_coord, field, fig=None, figsize=(7, 2), lines=True):
    """Plot the contours for a solution"""
    if fig is None:
        fig = plt.figure(figsize=figsize)
    f_max = max(np.max(field), -np.min(field))
    plt.contourf(x_coord, y_coord, field, vmin=-f_max, vmax=f_max)
    plt.colorbar()
    if lines:
        plt.contour(x_coord, y_coord, field, colors="k")
    plt.axis("image")
    return fig

Solutions

Cantilever beam with load at the end

A cantilever beam is a beam that is anchored at one end. In this case, the beam has length $L$ and height $h$. The force $F$ is applied at $x=L$.





In [5]:



    


L, h, F =  symbols("L h F")

To choose the Airy polynomial we need to fix the maximum power. This seem to be arbitrary, but we can use the intuition that we have grasped in the study of resistance of materials. Taking that into account, we should have stresses with powers of two, then the Airy stress functions should have at least powers of 4.

Picking 4 as maximum power will let us recover the desired solution, but in general we don't know this beforehand, and we should use a trial and error procedure. You could try with a maximum power of 10 for this example and it should still give you the same answer. It will take a while, though.





In [6]:



    


phi, c = poly_airy(x, y, 4)

phi



    


















    
Out[6]:







$\displaystyle c_{0} x^{4} + c_{1} x^{3} y + c_{10} x y + c_{11} y^{2} + c_{2} x^{2} y^{2} + c_{3} x y^{3} + c_{4} y^{4} + c_{5} x^{3} + c_{6} x^{2} y + c_{7} x y^{2} + c_{8} y^{3} + c_{9} x^{2}$
















In [7]:



    


s11, s22, s12 = airy_stress(phi, x, y)
display(s11, s22, s12)



    


















    







$\displaystyle 2 \left(c_{11} + c_{2} x^{2} + 3 c_{3} x y + 6 c_{4} y^{2} + c_{7} x + 3 c_{8} y\right)$









    







$\displaystyle 2 \left(6 c_{0} x^{2} + 3 c_{1} x y + c_{2} y^{2} + 3 c_{5} x + c_{6} y + c_{9}\right)$









    







$\displaystyle - 3 c_{1} x^{2} - c_{10} - 4 c_{2} x y - 3 c_{3} y^{2} - 2 c_{6} x - 2 c_{7} y$

The stress boundary conditions for this problem are:

\begin{align} &\sigma_{22}(x, y=h/2) = 0\\ &\sigma_{12}(x, y=h/2) = 0\\ &\sigma_{22}(x, y=-h/2) = 0\\ &\sigma_{12}(x, y=-h/2) = 0 \end{align}

Additionally, we have weak boundary conditions. That means, that we are not satisfying a distribution of tractions point-wise but in an average sense, namely

\begin{align} &\int_{-h/2}^{h/2} \sigma_{11} \mathrm{d}y = 0\\ &\int_{-h/2}^{h/2} \sigma_{12} \mathrm{d}y = F\\ &\int_{-h/2}^{h/2} y\sigma_{11} \mathrm{d}y = 0\, . \end{align}

With the weak conditions we are imposing that the system satisfy global equilibrium





In [8]:



    


t1 = s22.subs(y, h/2)
t2 = s12.subs(y, h/2)
t3 = s22.subs(y, -h/2)
t4 = s12.subs(y, -h/2)
t5 = s11.subs(x, L)
t6 = s12.subs(x, L)



    











In [9]:



    


conds = [t1, t2, t3, t4,
         integrate(t5, (y, -h/2, h/2)),
         integrate(t6, (y, -h/2, h/2)) - F,
         integrate(y*t5, (y, -h/2, h/2)),
         diff(phi, x, 4) + diff(phi, y, 4) + 2*diff(phi, x, 2, y, 2)]



    











In [10]:



    


eqlist = conds2eqs(conds, (x, y))
eqlist



    


















    
Out[10]:







$\displaystyle \left[ 12 c_{0}, \  2 L^{2} c_{2} h + 2 L c_{7} h + 2 c_{11} h + c_{4} h^{3}, \  24 c_{0} + 8 c_{2} + 24 c_{4}, \  \frac{L c_{3} h^{3}}{2} + \frac{c_{8} h^{3}}{2}, \  - c_{10} - \frac{3 c_{3} h^{2}}{4} - c_{7} h, \  \frac{c_{2} h^{2}}{2} + c_{6} h + 2 c_{9}, \  - F - 3 L^{2} c_{1} h - 2 L c_{6} h - c_{10} h - \frac{c_{3} h^{3}}{4}, \  - 2 c_{2} h - 2 c_{6}, \  3 c_{1} h + 6 c_{5}, \  - 3 c_{1} h + 6 c_{5}, \  2 c_{2} h - 2 c_{6}, \  - c_{10} - \frac{3 c_{3} h^{2}}{4} + c_{7} h, \  - 3 c_{1}, \  \frac{c_{2} h^{2}}{2} - c_{6} h + 2 c_{9}\right]$
















In [11]:



    


sol = solve(eqlist, c)



    











In [12]:



    


sol



    


















    
Out[12]:







$\displaystyle \left\{ c_{0} : 0, \  c_{1} : 0, \  c_{10} : - \frac{3 F}{2 h}, \  c_{11} : 0, \  c_{2} : 0, \  c_{3} : \frac{2 F}{h^{3}}, \  c_{4} : 0, \  c_{5} : 0, \  c_{6} : 0, \  c_{7} : 0, \  c_{8} : - \frac{2 F L}{h^{3}}, \  c_{9} : 0\right\}$
















In [13]:



    


s11.subs(sol)



    


















    
Out[13]:







$\displaystyle - \frac{12 F L y}{h^{3}} + \frac{12 F x y}{h^{3}}$
















In [14]:



    


s12.subs(sol)



    


















    
Out[14]:







$\displaystyle \frac{3 F}{2 h} - \frac{6 F y^{2}}{h^{3}}$
















In [15]:



    


s22.subs(sol)



    


















    
Out[15]:







$\displaystyle 0$
















In [16]:



    


x_vec, y_vec = np.mgrid[0:10:100j, -1:1:100j]
s11_fun = lambdify((x, y, h, L, F), s11.subs(sol), "numpy")
s12_fun = lambdify((x, y, h, L, F), s12.subs(sol), "numpy")



    











In [17]:



    


s11_vec = s11_fun(x_vec, y_vec, 2, 10, 1)
plot_contours(x_vec, y_vec, s11_vec);



    


















    




















    
























In [18]:



    


s12_vec = s12_fun(x_vec, y_vec, 2, 10, 1)
plot_contours(x_vec, y_vec, s12_vec, lines=False);

Beam under pure bending





In [19]:



    


L, h, S =  symbols("L h S")



    











In [20]:



    


phi, c = poly_airy(x, y, 3)

phi



    


















    
Out[20]:







$\displaystyle c_{0} x^{3} + c_{1} x^{2} y + c_{2} x y^{2} + c_{3} y^{3} + c_{4} x^{2} + c_{5} x y + c_{6} y^{2}$
















In [21]:



    


s11, s22, s12 = airy_stress(phi, x, y)



    











In [22]:



    


conds = [s22.subs(y, h/2),
         s12.subs(y, h/2),
         s22.subs(y, -h/2),
         s12.subs(y, -h/2),
         s11.subs(x, L/2) - S*y,
         s12.subs(x, L/2),
         s11.subs(x, -L/2) - S*y,
         s12.subs(x, -L/2)]

conds



    


















    
Out[22]:







$\displaystyle \left[ 6 c_{0} x + c_{1} h + 2 c_{4}, \  - 2 c_{1} x - c_{2} h - c_{5}, \  6 c_{0} x - c_{1} h + 2 c_{4}, \  - 2 c_{1} x + c_{2} h - c_{5}, \  L c_{2} - S y + 6 c_{3} y + 2 c_{6}, \  - L c_{1} - 2 c_{2} y - c_{5}, \  - L c_{2} - S y + 6 c_{3} y + 2 c_{6}, \  L c_{1} - 2 c_{2} y - c_{5}\right]$
















In [23]:



    


eqlist = conds2eqs(conds)

eqlist



    


















    
Out[23]:







$\displaystyle \left[ - c_{1} h + 2 c_{4}, \  - L c_{1} - c_{5}, \  c_{2} h - c_{5}, \  - 2 c_{2}, \  - S + 6 c_{3}, \  - L c_{2} + 2 c_{6}, \  c_{1} h + 2 c_{4}, \  - c_{2} h - c_{5}, \  - 2 c_{1}, \  L c_{2} + 2 c_{6}, \  6 c_{0}, \  L c_{1} - c_{5}\right]$
















In [24]:



    


sol = solve(eqlist, c)



    











In [25]:



    


sol



    


















    
Out[25]:







$\displaystyle \left\{ c_{0} : 0, \  c_{1} : 0, \  c_{2} : 0, \  c_{3} : \frac{S}{6}, \  c_{4} : 0, \  c_{5} : 0, \  c_{6} : 0\right\}$
















In [26]:



    


s11.subs(sol)



    


















    
Out[26]:







$\displaystyle S y$
















In [27]:



    


s12.subs(sol)



    


















    
Out[27]:







$\displaystyle 0$
















In [28]:



    


s22.subs(sol)



    


















    
Out[28]:







$\displaystyle 0$
















In [29]:



    


s11_fun = lambdify((x, y, S), s11.subs(sol), "numpy")



    











In [30]:



    


s11_vec = s11_fun(x_vec, y_vec, 1)
plot_contours(x_vec, y_vec, s11_vec);

Simply supported beam with uniform force





In [31]:



    


L, h, q =  symbols("L h q")



    











In [32]:



    


phi, d = poly_airy(x, y, 5)

phi



    


















    
Out[32]:







$\displaystyle c_{0} x^{5} + c_{1} x^{4} y + c_{10} y^{4} + c_{11} x^{3} + c_{12} x^{2} y + c_{13} x y^{2} + c_{14} y^{3} + c_{15} x^{2} + c_{16} x y + c_{17} y^{2} + c_{2} x^{3} y^{2} + c_{3} x^{2} y^{3} + c_{4} x y^{4} + c_{5} y^{5} + c_{6} x^{4} + c_{7} x^{3} y + c_{8} x^{2} y^{2} + c_{9} x y^{3}$
















In [33]:



    


s11, s22, s12 = airy_stress(phi, x, y)
display(s11, s22, s12)



    


















    







$\displaystyle 2 \left(6 c_{10} y^{2} + c_{13} x + 3 c_{14} y + c_{17} + c_{2} x^{3} + 3 c_{3} x^{2} y + 6 c_{4} x y^{2} + 10 c_{5} y^{3} + c_{8} x^{2} + 3 c_{9} x y\right)$









    







$\displaystyle 2 \left(10 c_{0} x^{3} + 6 c_{1} x^{2} y + 3 c_{11} x + c_{12} y + c_{15} + 3 c_{2} x y^{2} + c_{3} y^{3} + 6 c_{6} x^{2} + 3 c_{7} x y + c_{8} y^{2}\right)$









    







$\displaystyle - 4 c_{1} x^{3} - 2 c_{12} x - 2 c_{13} y - c_{16} - 6 c_{2} x^{2} y - 6 c_{3} x y^{2} - 4 c_{4} y^{3} - 3 c_{7} x^{2} - 4 c_{8} x y - 3 c_{9} y^{2}$
















In [34]:



    


conds = [s22.subs(y, h/2) + q,
         s12.subs(y, h/2),
         s22.subs(y, -h/2),
         s12.subs(y, -h/2),
         integrate(s11.subs(x, L/2), (y,-h/2, h/2)),
         integrate(s12.subs(x, L/2), (y,-h/2, h/2)) - q*L/2,
         integrate(y*s11.subs(x, L/2), (y,-h/2, h/2)),
         diff(phi, x, 4) + diff(phi, y, 4) + 2*diff(phi, x, 2, y, 2)]
conds



    


















    
Out[34]:







$\displaystyle \left[ 20 c_{0} x^{3} + 6 c_{1} h x^{2} + 6 c_{11} x + c_{12} h + 2 c_{15} + \frac{3 c_{2} h^{2} x}{2} + \frac{c_{3} h^{3}}{4} + 12 c_{6} x^{2} + 3 c_{7} h x + \frac{c_{8} h^{2}}{2} + q, \  - 4 c_{1} x^{3} - 2 c_{12} x - c_{13} h - c_{16} - 3 c_{2} h x^{2} - \frac{3 c_{3} h^{2} x}{2} - \frac{c_{4} h^{3}}{2} - 3 c_{7} x^{2} - 2 c_{8} h x - \frac{3 c_{9} h^{2}}{4}, \  20 c_{0} x^{3} - 6 c_{1} h x^{2} + 6 c_{11} x - c_{12} h + 2 c_{15} + \frac{3 c_{2} h^{2} x}{2} - \frac{c_{3} h^{3}}{4} + 12 c_{6} x^{2} - 3 c_{7} h x + \frac{c_{8} h^{2}}{2}, \  - 4 c_{1} x^{3} - 2 c_{12} x + c_{13} h - c_{16} + 3 c_{2} h x^{2} - \frac{3 c_{3} h^{2} x}{2} + \frac{c_{4} h^{3}}{2} - 3 c_{7} x^{2} + 2 c_{8} h x - \frac{3 c_{9} h^{2}}{4}, \  \frac{h^{3} \left(2 L c_{4} + 4 c_{10}\right)}{4} + h \left(\frac{L^{3} c_{2}}{4} + \frac{L^{2} c_{8}}{2} + L c_{13} + 2 c_{17}\right), \  - \frac{L q}{2} + \frac{h^{3} \left(- L c_{3} - c_{9}\right)}{4} + h \left(- \frac{L^{3} c_{1}}{2} - \frac{3 L^{2} c_{7}}{4} - L c_{12} - c_{16}\right), \  \frac{c_{5} h^{5}}{4} + \frac{h^{3} \left(\frac{L^{2} c_{3}}{2} + L c_{9} + 2 c_{14}\right)}{4}, \  24 \left(c_{10} + c_{4} x + 5 c_{5} y\right) + 24 \left(5 c_{0} x + c_{1} y + c_{6}\right) + 8 \left(3 c_{2} x + 3 c_{3} y + c_{8}\right)\right]$
















In [35]:



    


eqlist = conds2eqs(conds)

eqlist



    


















    
Out[35]:







$\displaystyle \left[ - 4 c_{1}, \  c_{12} h + 2 c_{15} + \frac{c_{3} h^{3}}{4} + \frac{c_{8} h^{2}}{2} + q, \  3 c_{2} h - 3 c_{7}, \  6 c_{11} + \frac{3 c_{2} h^{2}}{2} - 3 c_{7} h, \  \frac{L^{2} c_{3} h^{3}}{8} + \frac{L c_{9} h^{3}}{4} + \frac{c_{14} h^{3}}{2} + \frac{c_{5} h^{5}}{4}, \  6 c_{11} + \frac{3 c_{2} h^{2}}{2} + 3 c_{7} h, \  - 2 c_{12} - \frac{3 c_{3} h^{2}}{2} + 2 c_{8} h, \  20 c_{0}, \  - \frac{L^{3} c_{1} h}{2} - \frac{3 L^{2} c_{7} h}{4} - L c_{12} h - \frac{L c_{3} h^{3}}{4} - \frac{L q}{2} - c_{16} h - \frac{c_{9} h^{3}}{4}, \  24 c_{10} + 24 c_{6} + 8 c_{8}, \  - 3 c_{2} h - 3 c_{7}, \  \frac{L^{3} c_{2} h}{4} + \frac{L^{2} c_{8} h}{2} + L c_{13} h + \frac{L c_{4} h^{3}}{2} + c_{10} h^{3} + 2 c_{17} h, \  - c_{13} h - c_{16} - \frac{c_{4} h^{3}}{2} - \frac{3 c_{9} h^{2}}{4}, \  6 c_{1} h + 12 c_{6}, \  24 c_{1} + 24 c_{3} + 120 c_{5}, \  - 6 c_{1} h + 12 c_{6}, \  120 c_{0} + 24 c_{2} + 24 c_{4}, \  c_{13} h - c_{16} + \frac{c_{4} h^{3}}{2} - \frac{3 c_{9} h^{2}}{4}, \  - 2 c_{12} - \frac{3 c_{3} h^{2}}{2} - 2 c_{8} h, \  - c_{12} h + 2 c_{15} - \frac{c_{3} h^{3}}{4} + \frac{c_{8} h^{2}}{2}\right]$
















In [36]:



    


sol = solve(eqlist, d)



    











In [37]:



    


sol



    


















    
Out[37]:







$\displaystyle \left\{ c_{0} : 0, \  c_{1} : 0, \  c_{10} : 0, \  c_{11} : 0, \  c_{12} : - \frac{3 q}{4 h}, \  c_{13} : 0, \  c_{14} : - \frac{L^{2} q}{4 h^{3}} + \frac{q}{10 h}, \  c_{15} : - \frac{q}{4}, \  c_{16} : 0, \  c_{17} : 0, \  c_{2} : 0, \  c_{3} : \frac{q}{h^{3}}, \  c_{4} : 0, \  c_{5} : - \frac{q}{5 h^{3}}, \  c_{6} : 0, \  c_{7} : 0, \  c_{8} : 0, \  c_{9} : 0\right\}$
















In [38]:



    


factor(s11.subs(sol))



    


















    
Out[38]:







$\displaystyle \frac{q y \left(- 15 L^{2} + 6 h^{2} + 60 x^{2} - 40 y^{2}\right)}{10 h^{3}}$
















In [39]:



    


factor(s12.subs(sol))



    


















    
Out[39]:







$\displaystyle - \frac{3 q x \left(- h + 2 y\right) \left(h + 2 y\right)}{2 h^{3}}$
















In [40]:



    


factor(s22.subs(sol))



    


















    
Out[40]:







$\displaystyle \frac{q \left(- h + y\right) \left(h + 2 y\right)^{2}}{2 h^{3}}$
















In [41]:



    


x_vec, y_vec = np.mgrid[-5:5:100j, -1:1:100j]
s11_fun = lambdify((x, y, L, h, q), s11.subs(sol), "numpy")
s12_fun = lambdify((x, y, L, h, q), s12.subs(sol), "numpy")
s22_fun = lambdify((x, y, L, h, q), s22.subs(sol), "numpy")



    











In [42]:



    


s11_vec = s11_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s11_vec);



    


















    




















    
























In [43]:



    


s12_vec = s12_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s12_vec);



    


















    




















    
























In [44]:



    


s22_vec = s22_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s22_vec);

Simply supported beam with triangular force





In [45]:



    


L, h, q =  symbols("L h q")



    











In [46]:



    


phi, d = poly_airy(x, y, 6)

phi



    


















    
Out[46]:







$\displaystyle c_{0} x^{6} + c_{1} x^{5} y + c_{10} x^{2} y^{3} + c_{11} x y^{4} + c_{12} y^{5} + c_{13} x^{4} + c_{14} x^{3} y + c_{15} x^{2} y^{2} + c_{16} x y^{3} + c_{17} y^{4} + c_{18} x^{3} + c_{19} x^{2} y + c_{2} x^{4} y^{2} + c_{20} x y^{2} + c_{21} y^{3} + c_{22} x^{2} + c_{23} x y + c_{24} y^{2} + c_{3} x^{3} y^{3} + c_{4} x^{2} y^{4} + c_{5} x y^{5} + c_{6} y^{6} + c_{7} x^{5} + c_{8} x^{4} y + c_{9} x^{3} y^{2}$
















In [47]:



    


s11, s22, s12 = airy_stress(phi, x, y)
display(s11, s22, s12)



    


















    







$\displaystyle 2 \left(3 c_{10} x^{2} y + 6 c_{11} x y^{2} + 10 c_{12} y^{3} + c_{15} x^{2} + 3 c_{16} x y + 6 c_{17} y^{2} + c_{2} x^{4} + c_{20} x + 3 c_{21} y + c_{24} + 3 c_{3} x^{3} y + 6 c_{4} x^{2} y^{2} + 10 c_{5} x y^{3} + 15 c_{6} y^{4} + c_{9} x^{3}\right)$









    







$\displaystyle 2 \left(15 c_{0} x^{4} + 10 c_{1} x^{3} y + c_{10} y^{3} + 6 c_{13} x^{2} + 3 c_{14} x y + c_{15} y^{2} + 3 c_{18} x + c_{19} y + 6 c_{2} x^{2} y^{2} + c_{22} + 3 c_{3} x y^{3} + c_{4} y^{4} + 10 c_{7} x^{3} + 6 c_{8} x^{2} y + 3 c_{9} x y^{2}\right)$









    







$\displaystyle - 5 c_{1} x^{4} - 6 c_{10} x y^{2} - 4 c_{11} y^{3} - 3 c_{14} x^{2} - 4 c_{15} x y - 3 c_{16} y^{2} - 2 c_{19} x - 8 c_{2} x^{3} y - 2 c_{20} y - c_{23} - 9 c_{3} x^{2} y^{2} - 8 c_{4} x y^{3} - 5 c_{5} y^{4} - 4 c_{8} x^{3} - 6 c_{9} x^{2} y$
















In [48]:



    


conds = [s22.subs(y, h/2) + (x + L/2)*q/L,
         s12.subs(y, h/2),
         s22.subs(y, -h/2),
         s12.subs(y, -h/2),
         integrate(s11.subs(x, L/2), (y,-h/2, h/2)),
         integrate(s12.subs(x, L/2), (y,-h/2, h/2)) - q*L/2,
         integrate(y*s11.subs(x, L/2), (y,-h/2, h/2)),
         diff(phi, x, 4) + diff(phi, y, 4) + 2*diff(phi, x, 2, y, 2)]
conds



    


















    
Out[48]:







$\displaystyle \left[ 30 c_{0} x^{4} + 10 c_{1} h x^{3} + \frac{c_{10} h^{3}}{4} + 12 c_{13} x^{2} + 3 c_{14} h x + \frac{c_{15} h^{2}}{2} + 6 c_{18} x + c_{19} h + 3 c_{2} h^{2} x^{2} + 2 c_{22} + \frac{3 c_{3} h^{3} x}{4} + \frac{c_{4} h^{4}}{8} + 20 c_{7} x^{3} + 6 c_{8} h x^{2} + \frac{3 c_{9} h^{2} x}{2} + \frac{q \left(\frac{L}{2} + x\right)}{L}, \  - 5 c_{1} x^{4} - \frac{3 c_{10} h^{2} x}{2} - \frac{c_{11} h^{3}}{2} - 3 c_{14} x^{2} - 2 c_{15} h x - \frac{3 c_{16} h^{2}}{4} - 2 c_{19} x - 4 c_{2} h x^{3} - c_{20} h - c_{23} - \frac{9 c_{3} h^{2} x^{2}}{4} - c_{4} h^{3} x - \frac{5 c_{5} h^{4}}{16} - 4 c_{8} x^{3} - 3 c_{9} h x^{2}, \  30 c_{0} x^{4} - 10 c_{1} h x^{3} - \frac{c_{10} h^{3}}{4} + 12 c_{13} x^{2} - 3 c_{14} h x + \frac{c_{15} h^{2}}{2} + 6 c_{18} x - c_{19} h + 3 c_{2} h^{2} x^{2} + 2 c_{22} - \frac{3 c_{3} h^{3} x}{4} + \frac{c_{4} h^{4}}{8} + 20 c_{7} x^{3} - 6 c_{8} h x^{2} + \frac{3 c_{9} h^{2} x}{2}, \  - 5 c_{1} x^{4} - \frac{3 c_{10} h^{2} x}{2} + \frac{c_{11} h^{3}}{2} - 3 c_{14} x^{2} + 2 c_{15} h x - \frac{3 c_{16} h^{2}}{4} - 2 c_{19} x + 4 c_{2} h x^{3} + c_{20} h - c_{23} - \frac{9 c_{3} h^{2} x^{2}}{4} + c_{4} h^{3} x - \frac{5 c_{5} h^{4}}{16} - 4 c_{8} x^{3} + 3 c_{9} h x^{2}, \  \frac{3 c_{6} h^{5}}{8} + \frac{h^{3} \left(L^{2} c_{4} + 2 L c_{11} + 4 c_{17}\right)}{4} + h \left(\frac{L^{4} c_{2}}{8} + \frac{L^{3} c_{9}}{4} + \frac{L^{2} c_{15}}{2} + L c_{20} + 2 c_{24}\right), \  - \frac{L q}{2} - \frac{c_{5} h^{5}}{16} + \frac{h^{3} \left(- \frac{3 L^{2} c_{3}}{4} - L c_{10} - c_{16}\right)}{4} + h \left(- \frac{5 L^{4} c_{1}}{16} - \frac{L^{3} c_{8}}{2} - \frac{3 L^{2} c_{14}}{4} - L c_{19} - c_{23}\right), \  \frac{h^{5} \left(2 L c_{5} + 4 c_{12}\right)}{16} + \frac{h^{3} \left(\frac{L^{3} c_{3}}{4} + \frac{L^{2} c_{10}}{2} + L c_{16} + 2 c_{21}\right)}{4}, \  24 \left(15 c_{0} x^{2} + 5 c_{1} x y + c_{13} + c_{2} y^{2} + 5 c_{7} x + c_{8} y\right) + 8 \left(3 c_{10} y + c_{15} + 6 c_{2} x^{2} + 9 c_{3} x y + 6 c_{4} y^{2} + 3 c_{9} x\right) + 24 \left(c_{11} x + 5 c_{12} y + c_{17} + c_{4} x^{2} + 5 c_{5} x y + 15 c_{6} y^{2}\right)\right]$
















In [49]:



    


eqlist = conds2eqs(conds)

eqlist



    


















    
Out[49]:







$\displaystyle \left[ \frac{L^{4} c_{2} h}{8} + \frac{L^{3} c_{9} h}{4} + \frac{L^{2} c_{15} h}{2} + \frac{L^{2} c_{4} h^{3}}{4} + \frac{L c_{11} h^{3}}{2} + L c_{20} h + c_{17} h^{3} + 2 c_{24} h + \frac{3 c_{6} h^{5}}{8}, \  10 c_{1} h + 20 c_{7}, \  12 c_{13} + 3 c_{2} h^{2} - 6 c_{8} h, \  120 c_{1} + 72 c_{3} + 120 c_{5}, \  24 c_{10} + 120 c_{12} + 24 c_{8}, \  - 5 c_{1}, \  - 4 c_{2} h - 4 c_{8}, \  24 c_{13} + 8 c_{15} + 24 c_{17}, \  - \frac{5 L^{4} c_{1} h}{16} - \frac{L^{3} c_{8} h}{2} - \frac{3 L^{2} c_{14} h}{4} - \frac{3 L^{2} c_{3} h^{3}}{16} - \frac{L c_{10} h^{3}}{4} - L c_{19} h - \frac{L q}{2} - \frac{c_{16} h^{3}}{4} - c_{23} h - \frac{c_{5} h^{5}}{16}, \  24 c_{11} + 120 c_{7} + 24 c_{9}, \  - 3 c_{14} - \frac{9 c_{3} h^{2}}{4} + 3 c_{9} h, \  - 3 c_{14} - \frac{9 c_{3} h^{2}}{4} - 3 c_{9} h, \  - \frac{c_{10} h^{3}}{4} + \frac{c_{15} h^{2}}{2} - c_{19} h + 2 c_{22} + \frac{c_{4} h^{4}}{8}, \  24 c_{2} + 48 c_{4} + 360 c_{6}, \  \frac{c_{11} h^{3}}{2} - \frac{3 c_{16} h^{2}}{4} + c_{20} h - c_{23} - \frac{5 c_{5} h^{4}}{16}, \  - \frac{3 c_{10} h^{2}}{2} - 2 c_{15} h - 2 c_{19} - c_{4} h^{3}, \  \frac{12 L c_{14} h + 24 L c_{18} + 3 L c_{3} h^{3} + 6 L c_{9} h^{2} + 4 q}{4 L}, \  - \frac{c_{11} h^{3}}{2} - \frac{3 c_{16} h^{2}}{4} - c_{20} h - c_{23} - \frac{5 c_{5} h^{4}}{16}, \  12 c_{13} + 3 c_{2} h^{2} + 6 c_{8} h, \  - \frac{3 c_{10} h^{2}}{2} + 2 c_{15} h - 2 c_{19} + c_{4} h^{3}, \  \frac{L^{3} c_{3} h^{3}}{16} + \frac{L^{2} c_{10} h^{3}}{8} + \frac{L c_{16} h^{3}}{4} + \frac{L c_{5} h^{5}}{8} + \frac{c_{12} h^{5}}{4} + \frac{c_{21} h^{3}}{2}, \  4 c_{2} h - 4 c_{8}, \  - 3 c_{14} h + 6 c_{18} - \frac{3 c_{3} h^{3}}{4} + \frac{3 c_{9} h^{2}}{2}, \  360 c_{0} + 48 c_{2} + 24 c_{4}, \  - 10 c_{1} h + 20 c_{7}, \  30 c_{0}, \  \frac{c_{10} h^{3}}{4} + \frac{c_{15} h^{2}}{2} + c_{19} h + 2 c_{22} + \frac{c_{4} h^{4}}{8} + \frac{q}{2}\right]$
















In [50]:



    


sol = solve(eqlist, d)



    











In [51]:



    


sol



    


















    
Out[51]:







$\displaystyle \left\{ c_{0} : 0, \  c_{1} : 0, \  c_{10} : \frac{q}{2 h^{3}}, \  c_{11} : 0, \  c_{12} : - \frac{q}{10 h^{3}}, \  c_{13} : 0, \  c_{14} : - \frac{q}{4 L h}, \  c_{15} : 0, \  c_{16} : \frac{L q}{4 h^{3}} + \frac{q}{10 L h}, \  c_{17} : 0, \  c_{18} : - \frac{q}{12 L}, \  c_{19} : - \frac{3 q}{8 h}, \  c_{2} : 0, \  c_{20} : 0, \  c_{21} : \frac{q \left(- 35 L^{2} + 6 h^{2}\right)}{120 h^{3}}, \  c_{22} : - \frac{q}{8}, \  c_{23} : - \frac{q \left(15 L^{2} + h^{2}\right)}{80 L h}, \  c_{24} : 0, \  c_{3} : \frac{q}{3 L h^{3}}, \  c_{4} : 0, \  c_{5} : - \frac{q}{5 L h^{3}}, \  c_{6} : 0, \  c_{7} : 0, \  c_{8} : 0, \  c_{9} : 0\right\}$
















In [52]:



    


factor(s11.subs(sol))



    


















    
Out[52]:







$\displaystyle \frac{q y \left(- 35 L^{3} + 30 L^{2} x + 6 L h^{2} + 60 L x^{2} - 40 L y^{2} + 12 h^{2} x + 40 x^{3} - 80 x y^{2}\right)}{20 L h^{3}}$
















In [53]:



    


factor(s12.subs(sol))



    


















    
Out[53]:







$\displaystyle - \frac{q \left(- h + 2 y\right) \left(h + 2 y\right) \left(15 L^{2} + 60 L x + h^{2} + 60 x^{2} - 20 y^{2}\right)}{80 L h^{3}}$
















In [54]:



    


factor(s22.subs(sol))



    


















    
Out[54]:







$\displaystyle \frac{q \left(L + 2 x\right) \left(- h + y\right) \left(h + 2 y\right)^{2}}{4 L h^{3}}$
















In [55]:



    


x_vec, y_vec = np.mgrid[-5:5:100j, -1:1:100j]
s11_fun = lambdify((x, y, L, h, q), s11.subs(sol), "numpy")
s12_fun = lambdify((x, y, L, h, q), s12.subs(sol), "numpy")
s22_fun = lambdify((x, y, L, h, q), s22.subs(sol), "numpy")



    











In [56]:



    


s11_vec = s11_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s11_vec);



    


















    




















    
























In [57]:



    


s12_vec = s12_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s12_vec);



    


















    




















    
























In [58]:



    


s22_vec = s22_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s22_vec);

Cantilever beam with variying thickness





In [59]:



    


L, h, q =  symbols("L h q")



    











In [60]:



    


phi, d = poly_airy(x, y, 5)

phi



    


















    
Out[60]:







$\displaystyle c_{0} x^{5} + c_{1} x^{4} y + c_{10} y^{4} + c_{11} x^{3} + c_{12} x^{2} y + c_{13} x y^{2} + c_{14} y^{3} + c_{15} x^{2} + c_{16} x y + c_{17} y^{2} + c_{2} x^{3} y^{2} + c_{3} x^{2} y^{3} + c_{4} x y^{4} + c_{5} y^{5} + c_{6} x^{4} + c_{7} x^{3} y + c_{8} x^{2} y^{2} + c_{9} x y^{3}$
















In [61]:



    


s11, s22, s12 = airy_stress(phi, x, y)
display(s11, s22, s12)



    


















    







$\displaystyle 2 \left(6 c_{10} y^{2} + c_{13} x + 3 c_{14} y + c_{17} + c_{2} x^{3} + 3 c_{3} x^{2} y + 6 c_{4} x y^{2} + 10 c_{5} y^{3} + c_{8} x^{2} + 3 c_{9} x y\right)$









    







$\displaystyle 2 \left(10 c_{0} x^{3} + 6 c_{1} x^{2} y + 3 c_{11} x + c_{12} y + c_{15} + 3 c_{2} x y^{2} + c_{3} y^{3} + 6 c_{6} x^{2} + 3 c_{7} x y + c_{8} y^{2}\right)$









    







$\displaystyle - 4 c_{1} x^{3} - 2 c_{12} x - 2 c_{13} y - c_{16} - 6 c_{2} x^{2} y - 6 c_{3} x y^{2} - 4 c_{4} y^{3} - 3 c_{7} x^{2} - 4 c_{8} x y - 3 c_{9} y^{2}$
















In [62]:



    


conds = [s22.subs(y, h/2) + q,
         s12.subs(y, h/2),
         s22.subs(y, -h/2),
         s12.subs(y, -h/2),
         integrate(s11.subs(x, 0), (y,-h/2, h/2)),
         integrate(s12.subs(x, 0), (y,-h/2, h/2)) + q*L,
         integrate(y*s11.subs(x, 0), (y,-h/2, h/2)) - q*L**2/2,
         diff(phi, x, 4) + diff(phi, y, 4) + 2*diff(phi, x, 2, y, 2)]
conds



    


















    
Out[62]:







$\displaystyle \left[ 20 c_{0} x^{3} + 6 c_{1} h x^{2} + 6 c_{11} x + c_{12} h + 2 c_{15} + \frac{3 c_{2} h^{2} x}{2} + \frac{c_{3} h^{3}}{4} + 12 c_{6} x^{2} + 3 c_{7} h x + \frac{c_{8} h^{2}}{2} + q, \  - 4 c_{1} x^{3} - 2 c_{12} x - c_{13} h - c_{16} - 3 c_{2} h x^{2} - \frac{3 c_{3} h^{2} x}{2} - \frac{c_{4} h^{3}}{2} - 3 c_{7} x^{2} - 2 c_{8} h x - \frac{3 c_{9} h^{2}}{4}, \  20 c_{0} x^{3} - 6 c_{1} h x^{2} + 6 c_{11} x - c_{12} h + 2 c_{15} + \frac{3 c_{2} h^{2} x}{2} - \frac{c_{3} h^{3}}{4} + 12 c_{6} x^{2} - 3 c_{7} h x + \frac{c_{8} h^{2}}{2}, \  - 4 c_{1} x^{3} - 2 c_{12} x + c_{13} h - c_{16} + 3 c_{2} h x^{2} - \frac{3 c_{3} h^{2} x}{2} + \frac{c_{4} h^{3}}{2} - 3 c_{7} x^{2} + 2 c_{8} h x - \frac{3 c_{9} h^{2}}{4}, \  c_{10} h^{3} + 2 c_{17} h, \  L q - c_{16} h - \frac{c_{9} h^{3}}{4}, \  - \frac{L^{2} q}{2} + \frac{c_{14} h^{3}}{2} + \frac{c_{5} h^{5}}{4}, \  24 \left(c_{10} + c_{4} x + 5 c_{5} y\right) + 24 \left(5 c_{0} x + c_{1} y + c_{6}\right) + 8 \left(3 c_{2} x + 3 c_{3} y + c_{8}\right)\right]$
















In [63]:



    


eqlist = conds2eqs(conds)
eqlist



    


















    
Out[63]:







$\displaystyle \left[ - 4 c_{1}, \  c_{12} h + 2 c_{15} + \frac{c_{3} h^{3}}{4} + \frac{c_{8} h^{2}}{2} + q, \  3 c_{2} h - 3 c_{7}, \  6 c_{11} + \frac{3 c_{2} h^{2}}{2} - 3 c_{7} h, \  6 c_{11} + \frac{3 c_{2} h^{2}}{2} + 3 c_{7} h, \  - 2 c_{12} - \frac{3 c_{3} h^{2}}{2} + 2 c_{8} h, \  20 c_{0}, \  24 c_{10} + 24 c_{6} + 8 c_{8}, \  - \frac{L^{2} q}{2} + \frac{c_{14} h^{3}}{2} + \frac{c_{5} h^{5}}{4}, \  c_{10} h^{3} + 2 c_{17} h, \  - 3 c_{2} h - 3 c_{7}, \  - c_{13} h - c_{16} - \frac{c_{4} h^{3}}{2} - \frac{3 c_{9} h^{2}}{4}, \  6 c_{1} h + 12 c_{6}, \  24 c_{1} + 24 c_{3} + 120 c_{5}, \  L q - c_{16} h - \frac{c_{9} h^{3}}{4}, \  - 6 c_{1} h + 12 c_{6}, \  120 c_{0} + 24 c_{2} + 24 c_{4}, \  c_{13} h - c_{16} + \frac{c_{4} h^{3}}{2} - \frac{3 c_{9} h^{2}}{4}, \  - 2 c_{12} - \frac{3 c_{3} h^{2}}{2} - 2 c_{8} h, \  - c_{12} h + 2 c_{15} - \frac{c_{3} h^{3}}{4} + \frac{c_{8} h^{2}}{2}\right]$
















In [64]:



    


sol = solve(eqlist, d)



    











In [65]:



    


sol



    


















    
Out[65]:







$\displaystyle \left\{ c_{0} : 0, \  c_{1} : 0, \  c_{10} : 0, \  c_{11} : 0, \  c_{12} : - \frac{3 q}{4 h}, \  c_{13} : 0, \  c_{14} : \frac{L^{2} q}{h^{3}} + \frac{q}{10 h}, \  c_{15} : - \frac{q}{4}, \  c_{16} : \frac{3 L q}{2 h}, \  c_{17} : 0, \  c_{2} : 0, \  c_{3} : \frac{q}{h^{3}}, \  c_{4} : 0, \  c_{5} : - \frac{q}{5 h^{3}}, \  c_{6} : 0, \  c_{7} : 0, \  c_{8} : 0, \  c_{9} : - \frac{2 L q}{h^{3}}\right\}$
















In [66]:



    


factor(s11.subs(sol))



    


















    
Out[66]:







$\displaystyle \frac{q y \left(30 L^{2} - 60 L x + 3 h^{2} + 30 x^{2} - 20 y^{2}\right)}{5 h^{3}}$
















In [67]:



    


factor(s12.subs(sol))



    


















    
Out[67]:







$\displaystyle - \frac{3 q \left(- L + x\right) \left(- h + 2 y\right) \left(h + 2 y\right)}{2 h^{3}}$
















In [68]:



    


factor(s22.subs(sol))



    


















    
Out[68]:







$\displaystyle \frac{q \left(- h + y\right) \left(h + 2 y\right)^{2}}{2 h^{3}}$
















In [69]:



    


x_vec, y_vec = np.mgrid[0:10:100j, -1:1:100j]
s11_fun = lambdify((x, y, L, h, q), s11.subs(sol), "numpy")
s12_fun = lambdify((x, y, L, h, q), s12.subs(sol), "numpy")
s22_fun = lambdify((x, y, L, h, q), s22.subs(sol), "numpy")



    











In [70]:



    


s11_vec = s11_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s11_vec);



    


















    




















    
























In [71]:



    


s12_vec = s12_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s12_vec);



    


















    




















    
























In [72]:



    


s22_vec = s22_fun(x_vec, y_vec, 10, 2, 1)
plot_contours(x_vec, y_vec, s22_vec);

Cantilever beam with uniform force and triangular shape

We want to check that the stress function

$$\phi = C\left[r^2(\alpha - \theta) + \frac{1}{2}r^2 \sin(2\theta) - r^2\cos^2\theta \tan \alpha\right]$$

satisfy the boundary conditions for this problem.





In [73]:



    


r, theta, C, alpha, q = symbols("r  theta C alpha q")
phi = C*(r**2*(alpha - theta) + r**2*sin(2*theta)/2 - r**2*cos(theta)**2*tan(alpha))
phi



    


















    
Out[73]:







$\displaystyle C \left(r^{2} \left(\alpha - \theta\right) + \frac{r^{2} \sin{\left(2 \theta \right)}}{2} - r^{2} \cos^{2}{\left(\theta \right)} \tan{\left(\alpha \right)}\right)$
















In [74]:



    


Srr = simplify(1/r**2*diff(phi, theta, 2) + 1/r*diff(phi, r))
Srr



    


















    
Out[74]:







$\displaystyle C \left(2 \alpha - 2 \theta - 2 \sin^{2}{\left(\theta \right)} \tan{\left(\alpha \right)} - \sin{\left(2 \theta \right)}\right)$
















In [75]:



    


Stt = simplify(diff(phi, r, 2))
Stt



    


















    
Out[75]:







$\displaystyle C \left(2 \alpha - 2 \theta + \sin{\left(2 \theta \right)} - 2 \cos^{2}{\left(\theta \right)} \tan{\left(\alpha \right)}\right)$
















In [76]:



    


Srt = -simplify(diff(1/r*diff(phi, theta), r))
Srt



    


















    
Out[76]:







$\displaystyle - C \left(\sin{\left(2 \theta \right)} \tan{\left(\alpha \right)} + \cos{\left(2 \theta \right)} - 1\right)$

In the upper edge we have that

\begin{align} \sigma_{\theta\theta} &= -q\\ \sigma_{r\theta} &= 0 \end{align}

and in the lower one

\begin{align} \sigma_{\theta\theta} &= 0\\ \sigma_{r\theta} &= 0 \end{align}

The constant $C$ can be inferred from the first boundary condition, i.e.





In [77]:



    


sol = solve(Stt.subs(theta, 0) + q, C)
sol



    


















    
Out[77]:







$\displaystyle \left[ - \frac{q}{2 \alpha - 2 \tan{\left(\alpha \right)}}\right]$

And the solution is then





In [78]:



    


Srr2 = Srr.subs(C, sol[0])
Stt2 = Stt.subs(C, sol[0])
Srt2 = Srt.subs(C, sol[0])

display(Srr2)
display(Stt2)
display(Srt2)



    


















    







$\displaystyle - \frac{q \left(2 \alpha - 2 \theta - 2 \sin^{2}{\left(\theta \right)} \tan{\left(\alpha \right)} - \sin{\left(2 \theta \right)}\right)}{2 \alpha - 2 \tan{\left(\alpha \right)}}$









    







$\displaystyle - \frac{q \left(2 \alpha - 2 \theta + \sin{\left(2 \theta \right)} - 2 \cos^{2}{\left(\theta \right)} \tan{\left(\alpha \right)}\right)}{2 \alpha - 2 \tan{\left(\alpha \right)}}$









    







$\displaystyle \frac{q \left(\sin{\left(2 \theta \right)} \tan{\left(\alpha \right)} + \cos{\left(2 \theta \right)} - 1\right)}{2 \alpha - 2 \tan{\left(\alpha \right)}}$

And we need to check that the other boundary conditions are satisfied.





In [79]:



    


Srt.subs(theta, 0)



    


















    
Out[79]:







$\displaystyle 0$
















In [80]:



    


trigsimp(Stt.subs(theta, alpha))



    


















    
Out[80]:







$\displaystyle 0$

SymPy does not detect that the other expression is zero, but we can help rewriting the expression





In [81]:



    


simplify((Srt.subs(theta, alpha)).rewrite(sin))



    


















    
Out[81]:







$\displaystyle 0$

If we want to write the solution in Cartesian coordinates we can change the frame by

$$ [\sigma'] = [Q^T] [\sigma] [Q]$$

with

$$Q = \begin{bmatrix} \cos(\theta) &\sin(\theta)\\ -\sin(\theta) &\cos(\theta) \end{bmatrix}$$ 




In [82]:



    


Q = Matrix([
    [cos(theta), sin(theta)],
    [-sin(theta), cos(theta)]])
S_polar = Matrix([
    [Srr2, Srt2],
    [Srt2, Stt2]])



    











In [83]:



    


S_cart = simplify(Q.T *  S_polar * Q)
S_cart



    


















    
Out[83]:







$\displaystyle \left[\begin{matrix}\frac{q \left(- 2 \alpha + 2 \theta + \sin{\left(2 \theta \right)}\right)}{2 \left(\alpha - \tan{\left(\alpha \right)}\right)} & \frac{q \sin^{2}{\left(\theta \right)}}{\alpha - \tan{\left(\alpha \right)}}\\\frac{q \sin^{2}{\left(\theta \right)}}{\alpha - \tan{\left(\alpha \right)}} & \frac{q \left(- 2 \alpha + 2 \theta - \sin{\left(2 \theta \right)} + 2 \tan{\left(\alpha \right)}\right)}{2 \left(\alpha - \tan{\left(\alpha \right)}\right)}\end{matrix}\right]$
















In [84]:



    


angle = np.pi/10
t_vec, r_vec = np.mgrid[0:angle:100j, 0:5:100j]
x_vec = r_vec * np.cos(t_vec)
y_vec = r_vec * np.sin(t_vec)
s11_fun = lambdify((r, theta, alpha, q), S_cart[0, 0], "numpy")
s12_fun = lambdify((r, theta, alpha, q), S_cart[0, 1], "numpy")
s22_fun = lambdify((r, theta, alpha, q), S_cart[1, 1], "numpy")



    











In [85]:



    


s11_vec = s11_fun(r_vec, t_vec, angle, 1)
plot_contours(x_vec, y_vec, s11_vec)
plt.gca().invert_yaxis();



    


















    




















    























In [86]:



    


s12_vec = s12_fun(r_vec, t_vec, angle, 1)
plot_contours(x_vec, y_vec, s12_vec)
plt.gca().invert_yaxis();



    


















    




















    























In [87]:



    


s22_vec = s22_fun(r_vec, t_vec, angle, 1)
plot_contours(x_vec, y_vec, s22_vec)
plt.gca().invert_yaxis();



    


















    




















    























In [ ]:

The following cell changes the style of the notebook





In [88]:



    


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



    


















    
Out[88]:

Content source: nicoguaro/notebooks_examples

Similar notebooks:

notebook.community | gallery | about