9.2.2.c.py

import sympy as sp
import numpy as np
from sympy import Matrix, simplify, diff, eye, det
print("Extension")
theta = sp.symbols("theta")
X1,X2,X3 = sp.symbols("X_1 X_2 X_3")
X = Matrix([X1,X2,X3])
a,b,c = sp.symbols("a b c")
ea,eb,ec,t = sp.symbols("epsilon_a epsilon_b epsilon_c t")
f0 = sp.symbols("f_0")
xe1 = (1+t*ea)*X1
xe2 = (1-t*eb)*X2
xe3 = (1-t*ec)*X3
xe = Matrix([xe1, xe2, xe3])
F = Matrix([[diff(i,j) for j in X] for i in xe])
print("F =",F)
Elagrange = simplify(1/2*(F.T*F-eye(3)))
print("Lagrange Strain =",Elagrange)
Elagrangedot = diff(Elagrange,t)
print("Lagrange Strain rate =",Elagrangedot)
Fdot = simplify(diff(F,t))
print("Fdot =",Fdot)
L = simplify(Fdot*F.inv())
print("Velocity gradient (L) =",L)
Dd = simplify(1/2*(L+L.T))
print("Symmetric part of L", Dd)
Cauchy = Matrix([[t*f0/(b*c*(1-t*eb)*(1-t*ec)),0,0],[0,0,0],[0,0,0]])
print("Cauchy =",Cauchy)
P = simplify(det(F)*Cauchy*F.T.inv())
print("First Piola =",P)
S = simplify(F.inv()*P)
print("Second Piola =",S)
e1 = sum([Cauchy[i]*Dd[i] for i in range(9)])
print("Cauchy",Cauchy,"Dd",Dd)
print("D and Cauchy are energy conjugates: dU/dt", simplify(e1))
e2 = sum([P[i]*Fdot[i] for i in range(9)])
print("P",P,"Fdot",Fdot)
print("P and Fdot are energy conjugates: dW/dt", simplify(e2))
e2 = sum([S[i]*Elagrangedot[i] for i in range(9)])
print("S",S,"E_Green",Elagrangedot)
print("S and E_Green are energy conjugates: dW/dt", simplify(e2))
print("Rotation")
xe = xe.subs({t:1})
Q = Matrix([[sp.cos(theta*(t-1)), sp.sin(theta*(t-1)),0],
            [-sp.sin(theta*(t-1)),sp.cos(theta*(t-1)), 0],
            [0,0,1]])
Q = Q.subs({theta:sp.pi/2})
print("Q =", Q)
x = Q*xe
F = Matrix([[diff(i,j) for j in X] for i in x])
print("F =",F)
Elagrange = simplify(1/2*(F.T*F-eye(3)))
print("Lagrange Strain =",Elagrange)
Elagrangedot = diff(Elagrange,t)
print("Lagrange Strain rate =",Elagrangedot)
Fdot = simplify(diff(F,t))
print("Fdot =",Fdot)
L = simplify(Fdot*F.inv())
print("Velocity gradient (L) =",L)
Dd = simplify(1/2*(L+L.T))
print("Symmetric part of L", Dd)
sigmalocal = Cauchy.subs({t:1})
Cauchy = simplify(Q*sigmalocal*Q.T)
print("Cauchy =",Cauchy)
P = simplify(det(F)*Cauchy*F.T.inv())
print("First Piola =",P)
S = simplify(F.inv()*P)
print("Second Piola =",S)
e1 = sum([Cauchy[i]*Dd[i] for i in range(9)])
print("Cauchy",Cauchy,"Dd",Dd)
print("D and Cauchy are energy conjugates: dU/dt", simplify(e1))
e2 = sum([P[i]*Fdot[i] for i in range(9)])
print("P",P,"Fdot",Fdot)
print("P and Fdot are energy conjugates: dW/dt", simplify(e2))
e2 = sum([S[i]*Elagrangedot[i] for i in range(9)])
print("S",S,"E_Green",Elagrangedot)
print("S and E_Green are energy conjugates: dW/dt", simplify(e2))