Linear Shell solution

Init symbols for sympy



In [17]:

    
from sympy import *
from geom_util import *
from sympy.vector import CoordSys3D
import matplotlib.pyplot as plt
import sys
sys.path.append("../")

%matplotlib inline

%reload_ext autoreload
%autoreload 2
%aimport geom_util



In [18]:

    
# Any tweaks that normally go in .matplotlibrc, etc., should explicitly go here
%config InlineBackend.figure_format='retina'
plt.rcParams['figure.figsize'] = (12, 12)

plt.rc('text', usetex=True)
    
plt.rc('font', family='serif')

init_printing()



In [19]:

    
N = CoordSys3D('N')
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3", real = True, positive=True)
A,K,rho = symbols("A K rho")

Tymoshenko theory

$u_1 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u\left( \alpha_1 \right)+\alpha_3\gamma \left( \alpha_1 \right) $

$u_2 \left( \alpha_1, \alpha_2, \alpha_3 \right)=0 $

$u_3 \left( \alpha_1, \alpha_2, \alpha_3 \right)=w\left( \alpha_1 \right) $

$ \left( \begin{array}{c} u_1 \\ \frac { \partial u_1 } { \partial \alpha_1} \\ \frac { \partial u_1 } { \partial \alpha_2} \\ \frac { \partial u_1 } { \partial \alpha_3} \\ u_2 \\ \frac { \partial u_2 } { \partial \alpha_1} \\ \frac { \partial u_2 } { \partial \alpha_2} \\ \frac { \partial u_2 } { \partial \alpha_3} \\ u_3 \\ \frac { \partial u_3 } { \partial \alpha_1} \\ \frac { \partial u_3 } { \partial \alpha_2} \\ \frac { \partial u_3 } { \partial \alpha_3} \\ \end{array} \right) = T \cdot \left( \begin{array}{c} u \\ \frac { \partial u } { \partial \alpha_1} \\ \gamma \\ \frac { \partial \gamma } { \partial \alpha_1} \\ w \\ \frac { \partial w } { \partial \alpha_1} \\ \end{array} \right) $



In [20]:

    
T=zeros(12,6)
T[0,0]=1
T[0,2]=alpha3
T[1,1]=1
T[1,3]=alpha3
T[3,2]=1

T[8,4]=1
T[9,5]=1
T









    Out[20]:





$$\left[\begin{matrix}1 & 0 & \alpha_{3} & 0 & 0 & 0\\0 & 1 & 0 & \alpha_{3} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 1\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$



In [21]:

    
B=Matrix([[0, 1/(A*(K*alpha3 + 1)), 0, 0, 0, 0, 0, 0, K/(K*alpha3 + 1), 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1/(A*(K*alpha3 + 1)), 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [-K/(K*alpha3 + 1), 0, 0, 0, 0, 0, 0, 0, 0, 1/(A*(K*alpha3 + 1)), 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])
B









    Out[21]:





$$\left[\begin{array}{cccccccccccc}0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$



In [22]:

    
E=zeros(6,9)
E[0,0]=1
E[1,4]=1
E[2,8]=1
E[3,1]=1
E[3,3]=1
E[4,2]=1
E[4,6]=1
E[5,5]=1
E[5,7]=1
E









    Out[22]:





$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0\end{matrix}\right]$$



In [23]:

    
simplify(E*B*T)









    Out[23]:





$$\left[\begin{matrix}0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & \frac{\alpha_{3}}{A \left(K \alpha_{3} + 1\right)} & \frac{K}{K \alpha_{3} + 1} & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & \frac{1}{K \alpha_{3} + 1} & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)}\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$



In [24]:

    
mu = Symbol('mu')
la = Symbol('lambda')
C_tensor = getIsotropicStiffnessTensor(mu, la)
C = convertStiffnessTensorToMatrix(C_tensor)
C









    Out[24]:





$$\left[\begin{matrix}\lambda + 2 \mu & \lambda & \lambda & 0 & 0 & 0\\\lambda & \lambda + 2 \mu & \lambda & 0 & 0 & 0\\\lambda & \lambda & \lambda + 2 \mu & 0 & 0 & 0\\0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & \mu & 0\\0 & 0 & 0 & 0 & 0 & \mu\end{matrix}\right]$$



In [25]:

    
S=T.T*B.T*E.T*C*E*B*T*A*(1+alpha3*K)**2
S=simplify(S)
S









    Out[25]:





$$\left[\begin{matrix}A K^{2} \mu & 0 & - A K \mu & 0 & 0 & - K \mu\\0 & \frac{1}{A} \left(\lambda + 2 \mu\right) & 0 & \frac{\alpha_{3}}{A} \left(\lambda + 2 \mu\right) & K \left(\lambda + 2 \mu\right) & 0\\- A K \mu & 0 & A \mu & 0 & 0 & \mu\\0 & \frac{\alpha_{3}}{A} \left(\lambda + 2 \mu\right) & 0 & \frac{\alpha_{3}^{2}}{A} \left(\lambda + 2 \mu\right) & K \alpha_{3} \left(\lambda + 2 \mu\right) & 0\\0 & K \left(\lambda + 2 \mu\right) & 0 & K \alpha_{3} \left(\lambda + 2 \mu\right) & A K^{2} \left(\lambda + 2 \mu\right) & 0\\- K \mu & 0 & \mu & 0 & 0 & \frac{\mu}{A}\end{matrix}\right]$$



In [26]:

    
h=Symbol('h')
S_in = integrate(S*(1-alpha3*K+(alpha3**2)*K),(alpha3, -h/2, h/2))
S_in









    Out[26]:





$$\left[\begin{matrix}\frac{A \mu}{12} K^{3} h^{3} + A K^{2} h \mu & 0 & - \frac{A \mu}{12} K^{2} h^{3} - A K h \mu & 0 & 0 & - \frac{K^{2} \mu}{12} h^{3} - K h \mu\\0 & \frac{h^{3}}{12 A} \left(K \lambda + 2 K \mu\right) + \frac{h}{A} \left(\lambda + 2 \mu\right) & 0 & \frac{h^{3}}{12 A} \left(- K \lambda - 2 K \mu\right) & \frac{h^{3}}{4} \left(\frac{K^{2} \lambda}{3} + \frac{2 \mu}{3} K^{2}\right) + h \left(K \lambda + 2 K \mu\right) & 0\\- \frac{A \mu}{12} K^{2} h^{3} - A K h \mu & 0 & \frac{A K}{12} h^{3} \mu + A h \mu & 0 & 0 & \frac{K \mu}{12} h^{3} + h \mu\\0 & \frac{h^{3}}{12 A} \left(- K \lambda - 2 K \mu\right) & 0 & \frac{h^{5}}{80 A} \left(K \lambda + 2 K \mu\right) + \frac{h^{3}}{12 A} \left(\lambda + 2 \mu\right) & \frac{h^{3}}{4} \left(- \frac{K^{2} \lambda}{3} - \frac{2 \mu}{3} K^{2}\right) & 0\\0 & \frac{h^{3}}{4} \left(\frac{K^{2} \lambda}{3} + \frac{2 \mu}{3} K^{2}\right) + h \left(K \lambda + 2 K \mu\right) & 0 & \frac{h^{3}}{4} \left(- \frac{K^{2} \lambda}{3} - \frac{2 \mu}{3} K^{2}\right) & \frac{h^{3}}{4} \left(\frac{A \lambda}{3} K^{3} + \frac{2 A}{3} K^{3} \mu\right) + h \left(A K^{2} \lambda + 2 A K^{2} \mu\right) & 0\\- \frac{K^{2} \mu}{12} h^{3} - K h \mu & 0 & \frac{K \mu}{12} h^{3} + h \mu & 0 & 0 & \frac{K h^{3} \mu}{12 A} + \frac{h \mu}{A}\end{matrix}\right]$$



In [11]:

    
E,nu=symbols('E nu')
lambda_elastic=E*nu/((1+nu)*(1-2*nu))
mu_elastic=E/(2*(1+nu))
S_ins=simplify(S_in.subs(A,1).subs(la,lambda_elastic).subs(mu,mu_elastic))
S_ins









    Out[11]:





$$\left[\begin{matrix}\frac{E K}{2 \left(\nu + 1\right)} \left(- \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{2 \left(\nu + 1\right)} \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & 0 & \frac{E}{2 \left(\nu + 1\right)} \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right)\\0 & \frac{E}{K \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{K^{2} \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(- K h - \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & \frac{E}{\left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0\\\frac{E}{2 \left(\nu + 1\right)} \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{2 K \left(\nu + 1\right)} \left(- \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & 0 & \frac{E}{2 K \left(\nu + 1\right)} \left(- \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right)\\0 & \frac{E}{K^{2} \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(- K h - \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{K^{3} \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(K h + \log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & \frac{E}{K \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(- K h - \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0\\0 & \frac{E}{\left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{K \left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(- K h - \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & \frac{E K}{\left(\nu + 1\right) \left(2 \nu - 1\right)} \left(- \nu + 1\right) \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0\\\frac{E}{2 \left(\nu + 1\right)} \left(\log{\left (- \frac{K h}{2} + 1 \right )} - \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & \frac{E}{2 K \left(\nu + 1\right)} \left(- \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right) & 0 & 0 & \frac{E}{2 K \left(\nu + 1\right)} \left(- \log{\left (- \frac{K h}{2} + 1 \right )} + \log{\left (\frac{K h}{2} + 1 \right )}\right)\end{matrix}\right]$$



In [12]:

    
a11=E/(1-nu**2)
a44=5*E/(12*(1+nu))

AM=Matrix([[a11,0],[0,a44]])
strainT=Matrix([[1,alpha3,0],[0,0,1]])
AT=strainT.T*AM*strainT
integrate(AT,(alpha3, -h/2, h/2))









    Out[12]:





$$\left[\begin{matrix}\frac{E h}{- \nu^{2} + 1} & 0 & 0\\0 & - \frac{E h^{3}}{4 \left(3 \nu^{2} - 3\right)} & 0\\0 & 0 & \frac{5 E h}{12 \nu + 12}\end{matrix}\right]$$



In [27]:

    
M=Matrix([[rho, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, rho, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, rho, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])
M=T.T*M*T*A*(1+alpha3*K)
M









    Out[27]:





$$\left[\begin{matrix}A \rho \left(K \alpha_{3} + 1\right) & 0 & A \alpha_{3} \rho \left(K \alpha_{3} + 1\right) & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\A \alpha_{3} \rho \left(K \alpha_{3} + 1\right) & 0 & A \alpha_{3}^{2} \rho \left(K \alpha_{3} + 1\right) & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & A \rho \left(K \alpha_{3} + 1\right) & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$



In [28]:

    
M_in = integrate(M,(alpha3, -h/2, h/2))
M_in









    Out[28]:





$$\left[\begin{matrix}A h \rho & 0 & \frac{A K}{12} h^{3} \rho & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\\frac{A K}{12} h^{3} \rho & 0 & \frac{A \rho}{12} h^{3} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & A h \rho & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

Cartesian coordinates



In [29]:

    
import fem.geometry as g
import fem.model as m
import fem.material as mat
import fem.shell.shellsolver as s
import fem.shell.mesh1D as me
import plot

stiffness_matrix_func = lambdify([A, K, mu, la, h], S_in, "numpy")
mass_matrix_func = lambdify([A, K, rho, h], M_in, "numpy")


def stiffness_matrix(material, geometry, x1, x2, x3):
    A,K = geometry.get_A_and_K(x1,x2,x3)
    return stiffness_matrix_func(A, K, material.mu(), material.lam(), thickness)

def mass_matrix(material, geometry, x1, x2, x3):
    A,K = geometry.get_A_and_K(x1,x2,x3)
    return mass_matrix_func(A, K, material.rho, thickness)



def generate_layers(thickness, layers_count, material):
    layer_top = thickness / 2
    layer_thickness = thickness / layers_count
    layers = set()
    for i in range(layers_count):
        layer = m.Layer(layer_top - layer_thickness, layer_top, material, i)
        layers.add(layer)
        layer_top -= layer_thickness
    return layers


def solve(geometry, thickness, linear, N_width, N_height):
    layers_count = 1
    layers = generate_layers(thickness, layers_count, mat.IsotropicMaterial.steel())
    model = m.Model(geometry, layers, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
    mesh = me.Mesh1D.generate(width, layers, N_width, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
    lam, vec = s.solve(model, mesh, stiffness_matrix, mass_matrix)
    
    return lam, vec, mesh, geometry



width = 2
curvature = 0.8
thickness = 0.05

corrugation_amplitude = 0.05
corrugation_frequency = 20

# geometry = g.CorrugatedCylindricalPlate(width, curvature, corrugation_amplitude, corrugation_frequency)
geometry = g.CylindricalPlate(width, curvature)
# geometry = g.Plate(width)

N_width = 100
N_height = 4


lam, vec, mesh, geometry = solve(geometry, thickness, False, N_width, N_height)
results = s.convert_to_results(lam, vec, mesh, geometry)

results_index = 0
    
plot.plot_init_and_deformed_geometry_in_cartesian(results[results_index], 0, width, -thickness / 2, thickness / 2, 0, geometry.to_cartesian_coordinates)
to_print = 20
if (len(results) < to_print):
    to_print = len(results)
    
for i in range(to_print):
    print(results[i].rad_per_sec_to_Hz(results[i].freq))









    



/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)






    












    



118.45066894862332
274.2375759657385
522.2831935950938
691.967748051237
862.397740568122
1171.1989757756553
1597.2421140530585
1648.692881098383
2068.9424765831905
2567.1828944531217
3061.780917484792
3128.464302821149
3727.516058592909
4348.34772984477
4527.867084461731
5024.295504229749
5723.174999997173
5990.050171185151
6450.437088698879
7195.497443789092