Shells 2D

Init symbols for sympy



In [1]:

    
from sympy import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
x1, x2 = symbols("x_1 x_2")
alpha1, alpha2 = symbols("alpha_1 alpha_2")
R, L, ga, gv = symbols("R L g_a g_v")
init_printing()

Cylindrical coordinates



In [2]:

    
a1 = pi / 2 + (L / 2 - alpha1)/R

x = R * cos(a1)
y = R * sin(a1)

r = x*N.i + y*N.j

Curve in 2D coordinates system will be defined with the following vector $\vec{r}=\vec{r(\alpha_1)}$



In [3]:

    
r









    Out[3]:





$$(- R \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (R \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

Tangent to curve



In [4]:

    
v_temp = r.diff(alpha1)
dr_len = v_temp.magnitude()
v = v_temp / dr_len
v = trigsimp(v)
v_temp









    Out[4]:





$$(\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$



In [5]:

    
trigsimp(dr_len)









    Out[5]:





$$1$$

Normal to curve



In [6]:

    
n_temp = v.diff(alpha1)
k=trigsimp(n_temp.magnitude())
n = n_temp/k
q=1/(R*sqrt(1/R**2))
n = trigsimp(n).subs(q, 1)
n









    Out[6]:





$$(\sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$



In [7]:

    
v.dot(n)









    Out[7]:





$$0$$



In [8]:

    
n.dot(v)









    Out[8]:





$$0$$

Curvature



In [9]:

    
sympify(k)









    Out[9]:





$$\sqrt{\frac{1}{R^{2}}}$$

Derivative of base vectors

Let's find $\frac { d\vec{n} } { d\alpha_1}$ $\frac { d\vec{v} } { d\alpha_1}$ $\frac { d\vec{n} } { d\alpha_2}$ $\frac { d\vec{v} } { d\alpha_2}$



In [10]:

    
n.diff(alpha1)









    Out[10]:





$$(- \frac{1}{R} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \frac{1}{R} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

$ \frac { d\vec{n} } { d\alpha_1} = -\frac {1}{R} \vec{v} = -k \vec{v} $



In [11]:

    
v.diff(alpha1)









    Out[11]:





$$(\frac{1}{R} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \frac{1}{R} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

$ \frac { d\vec{v} } { d\alpha_1} = \frac {1}{R} \vec{n} = k \vec{n} $

Derivative of vectors

$ \vec{u} = u_v \vec{v} + u_n\vec{n} $

$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u_v\vec{v}) } { d\alpha_1} + \frac { d(u_n\vec{n}) } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_1} + \frac { du_v } { d\alpha_1} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_1} = \frac { du_n } { d\alpha_1} \vec{n} - u_n k \vec{v} + \frac { du_v } { d\alpha_1} \vec{v} + u_v k \vec{n}$

Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du_v } { d\alpha_1} - u_n k \right) \vec{v} + \left( \frac { du_n } { d\alpha_1} + u_v k \right) \vec{n}$

$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u_n\vec{n}) } { d\alpha_2} + \frac { d(u_v\vec{v}) } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + u_n \frac { d\vec{n} } { d\alpha_2} + \frac { du_v } { d\alpha_2} \vec{v} + u_v \frac { d\vec{v} } { d\alpha_2} = \frac { du_n } { d\alpha_2} \vec{n} + \frac { du_v } { d\alpha_2} \vec{v} $

Base Vectors $\vec{R}_1, \vec{R}_2$



In [12]:

    
R_alpha=r+alpha2*n
R_alpha









    Out[12]:





$$(- R \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + \alpha_{2} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (R \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - \alpha_{2} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$



In [13]:

    
R1=R_alpha.diff(alpha1)
R2=R_alpha.diff(alpha2)



In [14]:

    
trigsimp(R1)









    Out[14]:





$$(\left(1 - \frac{\alpha_{2}}{R}\right) \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\left(1 - \frac{\alpha_{2}}{R}\right) \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$



In [15]:

    
R2









    Out[15]:





$$(\sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

Let's find Jacobi matrix:

$ A = \left( \begin{array}{cc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} \end{array} \right)$

$ \left[ \begin{array}{cc} \vec{R}_1 & \vec{R}_2 \end{array} \right] = \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] \cdot \left( \begin{array}{cc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} \end{array} \right) = \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] \cdot A$

$ \left[ \begin{array}{cc} \vec{e}_1 & \vec{e}_2 \end{array} \right] = \left[ \begin{array}{cc} \vec{R}_1 & \vec{R}_2 \end{array} \right] \cdot A^{-1}$



In [16]:

    
m11=R1.dot(N.i)
m12=R2.dot(N.i)
m21=R1.dot(N.j)
m22=R2.dot(N.j)
A=Matrix([[m11, m12], [m21, m22]])
A









    Out[16]:





$$\left[\begin{matrix}\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - \frac{\alpha_{2}}{R} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\\\sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - \frac{\alpha_{2}}{R} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & - \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\end{matrix}\right]$$



In [17]:

    
A_inv = trigsimp(A**-1)
sympify(trigsimp(Matrix([R1, R2]).T*A_inv))









    Out[17]:





$$\left[\begin{matrix}(\sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} & (\frac{1}{2} \left(- \cos{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + 1\right))\mathbf{\hat{j}_{N}} + (\frac{1}{2} \left(\cos{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + 1\right))\mathbf{\hat{j}_{N}}\end{matrix}\right]$$



In [18]:

    
trigsimp(A.det())









    Out[18]:





$$-1 + \frac{\alpha_{2}}{R}$$

Metric tensor



In [19]:

    
g11=R1.dot(R1)
g12=R1.dot(R2)
g21=R2.dot(R1)
g22=R2.dot(R2)

G=Matrix([[g11, g12],[g21, g22]])
G=trigsimp(G)
G









    Out[19]:





$$\left[\begin{matrix}\frac{1}{R^{2}} \left(R - \alpha_{2}\right)^{2} & 0\\0 & 1\end{matrix}\right]$$



In [20]:

    
G_inv = G**-1

Derivative of base vectors



In [21]:

    
dR1dalpha1 = trigsimp(R1.diff(alpha1))
dR1dalpha1









    Out[21]:





$$(\frac{1}{R} \left(1 - \frac{\alpha_{2}}{R}\right) \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\frac{1}{R} \left(-1 + \frac{\alpha_{2}}{R}\right) \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

$ \frac { d\vec{R_1} } { d\alpha_1} = \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \vec{R_2} $



In [22]:

    
dR1dalpha2 = trigsimp(R1.diff(alpha2))
dR1dalpha2









    Out[22]:





$$(- \frac{1}{R} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \frac{1}{R} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

$ \frac { d\vec{R_1} } { d\alpha_2} = -\frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} $



In [23]:

    
dR2dalpha1 = trigsimp(R2.diff(alpha1))
dR2dalpha1









    Out[23]:





$$(- \frac{1}{R} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (- \frac{1}{R} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{j}_{N}}$$

$ \frac { d\vec{R_2} } { d\alpha_1} = -\frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} $



In [24]:

    
dR2dalpha2 = trigsimp(R2.diff(alpha2))
dR2dalpha2









    Out[24]:





$$\mathbf{\hat{0}}$$

$ \frac { d\vec{R_2} } { d\alpha_2} = \vec{0} $

Derivative of vectors

$ \vec{u} = u^1 \vec{R_1} + u^2\vec{R_2} $

$ \frac { d\vec{u} } { d\alpha_1} = \frac { d(u^1\vec{R_1}) } { d\alpha_1} + \frac { d(u^2\vec{R_2}) } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \vec{R_2} + \frac { du^2 } { d\alpha_1} \vec{R_2} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1}$

Then $ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du^1 } { d\alpha_1} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \right) \vec{R_1} + \left( \frac { du^2 } { d\alpha_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) \right) \vec{R_2}$

$ \frac { d\vec{u} } { d\alpha_2} = \frac { d(u^1\vec{R_1}) } { d\alpha_2} + \frac { d(u^2\vec{R_2}) } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_2} + \frac { du^2 } { d\alpha_2} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} $

Then $ \frac { d\vec{u} } { d\alpha_2} = \left( \frac { du^1 } { d\alpha_2} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}} \right) \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2}$

$\nabla_1 u^1 = \frac { \partial u^1 } { \partial \alpha_1} - u^2 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}}$

$\nabla_1 u^2 = \frac { \partial u^2 } { \partial \alpha_1} + u^1 \frac {1}{R} \left( 1-\frac{\alpha_2}{R} \right) $

$\nabla_2 u^1 = \frac { \partial u^1 } { \partial \alpha_2} - u^1 \frac {1}{R} \frac {1}{1-\frac{\alpha_2}{R}}$

$\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$

$ \nabla \vec{u} = \left( \begin{array}{cc} \nabla_1 u^1 & \nabla_1 u^2 \\ \nabla_2 u^1 & \nabla_2 u^2 \end{array} \right)$



In [25]:

    
u1=Function('u^1')
u2=Function('u^2')

u1_nabla1 = u1(alpha1, alpha2).diff(alpha1) - u2(alpha1, alpha2) / R * (S(1)/(1-alpha2/R))
u2_nabla1 = u2(alpha1, alpha2).diff(alpha1) + u1(alpha1, alpha2) / R * ( 1-alpha2/R) 
u1_nabla2 = u1(alpha1, alpha2).diff(alpha2) - u1(alpha1, alpha2) / R * (S(1)/(1-alpha2/R))
u2_nabla2 = u2(alpha1, alpha2).diff(alpha2)
# $\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$

grad_u = Matrix([[u1_nabla1, u2_nabla1],[u1_nabla2, u2_nabla2]])
grad_u









    Out[25]:





$$\left[\begin{matrix}\frac{\partial}{\partial \alpha_{1}} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )} - \frac{\operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )}}{R \left(1 - \frac{\alpha_{2}}{R}\right)} & \frac{\partial}{\partial \alpha_{1}} \operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )} + \frac{1}{R} \left(1 - \frac{\alpha_{2}}{R}\right) \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )}\\\frac{\partial}{\partial \alpha_{2}} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )} - \frac{\operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )}}{R \left(1 - \frac{\alpha_{2}}{R}\right)} & \frac{\partial}{\partial \alpha_{2}} \operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )}\end{matrix}\right]$$



In [26]:

    
q=Symbol('q')
grad_u_down=grad_u.subs(1-alpha2/R, q)*G.subs((R-alpha2)/R,q)
#grad_u_down=grad_u*G
expand(simplify(grad_u_down))#.subs((R-alpha2)/R, q)









    Out[26]:





$$\left[\begin{matrix}q^{2} \frac{\partial}{\partial \alpha_{1}} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )} - \frac{q}{R} \operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )} & \frac{\partial}{\partial \alpha_{1}} \operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )} + \frac{q}{R} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )}\\q^{2} \frac{\partial}{\partial \alpha_{2}} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )} - \frac{q}{R} \operatorname{u^{1}}{\left (\alpha_{1},\alpha_{2} \right )} & \frac{\partial}{\partial \alpha_{2}} \operatorname{u^{2}}{\left (\alpha_{1},\alpha_{2} \right )}\end{matrix}\right]$$

$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \end{array}

\right)

\left( \begin{array}{c} \left( 1-\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_1} - u^2 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \left( 1-\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_2} - u^1 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \frac { \partial u^2 } { \partial \alpha_1} + u^1 \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} \\ \frac { \partial u^2 } { \partial \alpha_2} \end{array} \right) $

$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \end{array}

\right)

\left( \begin{array}{cccccc} 0 & \left( 1-\frac{\alpha_2}{R} \right)^2 & 0 & -\frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & 0 \\ -\frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & \left( 1-\frac{\alpha_2}{R} \right)^2 & 0 & 0 & 0 \\ \frac {\left( 1-\frac{\alpha_2}{R} \right)}{R} & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \end{array} \right) $

Elasticity tensor(stiffness tensor)



In [51]:

    
from sympy import MutableDenseNDimArray
C_x = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
                el = Symbol(elem_index)
                C_x[i,j,k,l] = el
                
C_x









    Out[51]:





$$\left[\begin{matrix}\left[\begin{matrix}C^{1111} & C^{1112} & C^{1113}\\C^{1121} & C^{1122} & C^{1123}\\C^{1131} & C^{1132} & C^{1133}\end{matrix}\right] & \left[\begin{matrix}C^{1211} & C^{1212} & C^{1213}\\C^{1221} & C^{1222} & C^{1223}\\C^{1231} & C^{1232} & C^{1233}\end{matrix}\right] & \left[\begin{matrix}C^{1311} & C^{1312} & C^{1313}\\C^{1321} & C^{1322} & C^{1323}\\C^{1331} & C^{1332} & C^{1333}\end{matrix}\right]\\\left[\begin{matrix}C^{2111} & C^{2112} & C^{2113}\\C^{2121} & C^{2122} & C^{2123}\\C^{2131} & C^{2132} & C^{2133}\end{matrix}\right] & \left[\begin{matrix}C^{2211} & C^{2212} & C^{2213}\\C^{2221} & C^{2222} & C^{2223}\\C^{2231} & C^{2232} & C^{2233}\end{matrix}\right] & \left[\begin{matrix}C^{2311} & C^{2312} & C^{2313}\\C^{2321} & C^{2322} & C^{2323}\\C^{2331} & C^{2332} & C^{2333}\end{matrix}\right]\\\left[\begin{matrix}C^{3111} & C^{3112} & C^{3113}\\C^{3121} & C^{3122} & C^{3123}\\C^{3131} & C^{3132} & C^{3133}\end{matrix}\right] & \left[\begin{matrix}C^{3211} & C^{3212} & C^{3213}\\C^{3221} & C^{3222} & C^{3223}\\C^{3231} & C^{3232} & C^{3233}\end{matrix}\right] & \left[\begin{matrix}C^{3311} & C^{3312} & C^{3313}\\C^{3321} & C^{3322} & C^{3323}\\C^{3331} & C^{3332} & C^{3333}\end{matrix}\right]\end{matrix}\right]$$



In [53]:

    
C_x_1 = MutableDenseNDimArray.zeros(3, 3, 3, 3)

def getCIndecies(index):
    if (index == 0):
        return 0, 0
    elif (index == 1):
        return 1, 1
    elif (index == 2):
        return 2, 2
    elif (index == 3):
        return 0, 1
    elif (index == 4):
        return 0, 2
    elif (index == 5):
        return 1, 2

for i in range(3):
    for j in range(i, 3):        
        for k in range(3):
            for l in range(k, 3):
                elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
                el = Symbol(elem_index)
                C_x_1[i,j,k,l] = el
                C_x_1[i,j,l,k] = el
                C_x_1[j,i,k,l] = el
                C_x_1[j,i,l,k] = el
                if (i >= k or j >= l):
                    C_x_1[k,l,i,j] = el
                    C_x_1[k,l,j,i] = el
                    C_x_1[l,k,i,j] = el
                    C_x_1[l,k,j,i] = el
                
C_x_1









    Out[53]:





$$\left[\begin{matrix}\left[\begin{matrix}C^{1111} & C^{1211} & C^{1311}\\C^{1211} & C^{2211} & C^{2311}\\C^{1311} & C^{2311} & C^{3311}\end{matrix}\right] & \left[\begin{matrix}C^{1211} & C^{1212} & C^{1312}\\C^{1212} & C^{2212} & C^{2312}\\C^{1312} & C^{2312} & C^{3312}\end{matrix}\right] & \left[\begin{matrix}C^{1311} & C^{1312} & C^{1313}\\C^{1312} & C^{2213} & C^{2313}\\C^{1313} & C^{2313} & C^{3313}\end{matrix}\right]\\\left[\begin{matrix}C^{1211} & C^{1212} & C^{1312}\\C^{1212} & C^{2212} & C^{2312}\\C^{1312} & C^{2312} & C^{3312}\end{matrix}\right] & \left[\begin{matrix}C^{2211} & C^{2212} & C^{2213}\\C^{2212} & C^{2222} & C^{2322}\\C^{2213} & C^{2322} & C^{3322}\end{matrix}\right] & \left[\begin{matrix}C^{2311} & C^{2312} & C^{2313}\\C^{2312} & C^{2322} & C^{2323}\\C^{2313} & C^{2323} & C^{3323}\end{matrix}\right]\\\left[\begin{matrix}C^{1311} & C^{1312} & C^{1313}\\C^{1312} & C^{2213} & C^{2313}\\C^{1313} & C^{2313} & C^{3313}\end{matrix}\right] & \left[\begin{matrix}C^{2311} & C^{2312} & C^{2313}\\C^{2312} & C^{2322} & C^{2323}\\C^{2313} & C^{2323} & C^{3323}\end{matrix}\right] & \left[\begin{matrix}C^{3311} & C^{3312} & C^{3313}\\C^{3312} & C^{3322} & C^{3323}\\C^{3313} & C^{3323} & C^{3333}\end{matrix}\right]\end{matrix}\right]$$



In [ ]:

    
# for i in range(3):
#     for j in range(i, 3):        
#         for k in range(3):
#             for l in range(k, 3):
#                 elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
#                 el = Symbol(elem_index)
#                 C_x[i,j,k,l] = el
#                 C_x[i,j,l,k] = el
#                 C_x[j,i,k,l] = el
#                 C_x[j,i,l,k] = el
                
                
# C_x

# for i in range(3):
#     for j in range(3):        
#         for k in range(3):
#             for l in range(3):
#                 elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
#                 el = Symbol(elem_index)
#                 C_x[i,j,k,l] = el
                
#                 C_x[k,l,i,j] = el
#         if (s < 3 and t < 3):
            
#         elif (s==t):
#             i,j = getCIndecies(s)
#             k,l = getCIndecies(t)
#             elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
#             el = Symbol(elem_index)
#             C_x[i,j,k,l] = el

# def getCIndecies(index):
#     if (index == 0):
#         return 0, 0
#     elif (index == 1):
#         return 1, 1
#     elif (index == 2):
#         return 2, 2
#     elif (index == 3):
#         return 0, 1
#     elif (index == 4):
#         return 0, 2
#     elif (index == 5):
#         return 1, 2
    

# def getCalpha(C, A, q, p, s, t):
#     res = S(0)
#     for i in range(3):
#         for j in range(3):        
#             for k in range(3):
#                 for l in range(3):
#                     res += C[i,j,k,l]*A[q,i]
#     return simplify(trigsimp(res))
                    


# C_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)
# C_alpha_empty = MutableDenseNDimArray.zeros(3, 3, 3, 3)

# m11=R1.dot(N.i)
# m12=R2.dot(N.i)
# m21=R1.dot(N.j)
# m22=R2.dot(N.j)
# A=Matrix([[m11, m12, 0], [m21, m22, 0], [0,0,1]])
# A_inv=A**-1

# for s in range(6):
#     for t in range(s, 6):
#         if (s < 3 and t < 3):
#             i,j = getCIndecies(s)
#             k,l = getCIndecies(t)
#             elem_index = 'C^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
#             el = Symbol(elem_index)
#             C_x[i,j,k,l] = el
#             C_x[k,l,i,j] = el
#         elif (s==t):
#             i,j = getCIndecies(s)
#             k,l = getCIndecies(t)
#             elem_index = 'C^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
#             el = Symbol(elem_index)
#             C_x[i,j,k,l] = el

# for i in range(3):
#     for j in range(3):        
#         for k in range(3):
#             for l in range(3):
#                 c = getCalpha(C_x, A_inv, i, j, k, l)
#                 C_alpha[i,j,k,l] = c
#                 if (c != 0):
#                     elem_index = 'C_{{\alpha}}^{}{}{}{}'.format(i+1, j+1, k+1, l+1)
#                     el = Symbol(elem_index)
#                     C_alpha_empty[i,j,k,l] = el
                
# # trigsimp(C_alpha[0,2,0,2])

# C_alpha_empty



In [31]:

    
C_x









    Out[31]:





$$\left[\begin{matrix}\left[\begin{matrix}C^{1111} & 0 & 0\\0 & C^{1122} & 0\\0 & 0 & C^{1133}\end{matrix}\right] & \left[\begin{matrix}0 & C^{1212} & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & C^{1313}\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\\\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}C^{1122} & 0 & 0\\0 & C^{2222} & 0\\0 & 0 & C^{2233}\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & C^{2323}\\0 & 0 & 0\end{matrix}\right]\\\left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}C^{1133} & 0 & 0\\0 & C^{2233} & 0\\0 & 0 & C^{3333}\end{matrix}\right]\end{matrix}\right]$$

Virtual work



In [ ]:

    
def contraction3DSameRank(A,B):
    res = S(0)
    for i in range(3):
        for j in range(3):
            res += A[i,j]*B[j,i]
    return res

def contraction3D(C,e):
    res = MutableDenseNDimArray.zeros(3, 3)
    for i in range(3):
        for j in range(3):
            res[i,j] = S(0)
            for k in range(3):
                for l in range(3):
                    res[i,j] += C[i,j,k,l]*e[k,l]
    return res



In [ ]:

    
e11 = Symbol("e_{11}")
e12 = Symbol("e_{12}")
e22 = Symbol("e_{22}")
e13 = Symbol("e_{13}")
e23 = Symbol("e_{23}")
e33 = Symbol("e_{33}")
# s11 = Symbol("s_{11}")
# s12 = Symbol("s_{12}")
# s22 = Symbol("s_{22}")

e=Matrix([[e11, e12, e13], [e12, e22, e23], [e13, e23, e33]])

s=contraction3D(C_alpha, e)

E=contraction3DSameRank(s, e)



In [ ]:

    
# e_alpha=G*A_inv*e*A_inv.T*G
# #e_alpha=A*e*A.T

# s_alpha=A_inv*s*A_inv.T
# E_alpha=contraction2D(s_alpha, e_alpha)