forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
10.1.3.1.py
46 lines (46 loc) · 1.75 KB
/
10.1.3.1.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import sympy as sp
from sympy import Matrix, diff, integrate, simplify
x1,x2,x3 = sp.symbols("x_1 x_2 x_3")
s = Matrix([[x1*x2,5,0],[5,x1,0],[0,0,0]])
x = Matrix([[x1,x2,x3]])
pb1,pb2,pb3 =[-sum([diff(s[i,j],x[i]) for i in range(3)]) for j in range(3)]
print("\u03C1b_1 =",pb1)
print("\u03C1b_2 =",pb2)
print("\u03C1b_3 =",pb3)
a,b,t = sp.symbols("a b t")
ustar = Matrix([a*x1 + b*x2, 0])
print("u* =",ustar)
gradustar = Matrix([[diff(i,j) for j in x[:2]]for i in ustar])
print("\u2207u* =",gradustar)
ustara = ustar.subs({x1:2})
ustarb = ustar.subs({x2:1})
ustarc = ustar.subs({x1:0})
ustard = ustar.subs({x2:0})
s = Matrix([[x1*x2,5],[5,x1]])
na = Matrix([1,0])
nb = Matrix([0,1])
nc = Matrix([-1,0])
nd = Matrix([0,-1])
ta = (s*na).subs({x1:2})
tb = (s*nb).subs({x2:1})
tc = (s*nc).subs({x1:0})
td = (s*nd).subs({x2:0})
print("t_A =",ta,"t_B =",tb,"t_C =",tc,"t_D =",td)
estar = sp.Rational("1/2")*(gradustar+gradustar.T)
print("\u03B5* =",estar)
ststrain = sum([s[i]*estar[i] for i in range(4)])
print("ststrain =",ststrain)
IVW = integrate(ststrain,(x1,0,2),(x2,0,1),(x3,0,t))
print("Internal Virtual Work =",IVW)
EVWBodyForces = integrate(pb1*ustar[0],(x1,0,2),(x2,0,1),(x3,0,t))
print("External Virtual Work|_{Body Forces} =",EVWBodyForces)
EVWa = integrate((ta.dot(ustara)).subs({x1:2}),(x2,0,1),(x3,0,t))
print("External Virtual Work|_A =",EVWa)
EVWb = integrate((tb.dot(ustarb)).subs({x2:1}),(x1,0,2),(x3,0,t))
print("External Virtual Work|_B =",EVWb)
EVWc = integrate((tc.dot(ustarc)).subs({x1:0}),(x2,0,1),(x3,0,t))
print("External Virtual Work|_C =",EVWc)
EVWd = integrate((td.dot(ustard)).subs({x2:0}),(x1,0,2),(x3,0,t))
print("External Virtual Work|_D =",EVWd)
EVW = simplify(EVWBodyForces+EVWa+EVWb+EVWc+EVWd)
print("External Virtual Work = Internal Virtual Work =",EVW)