Skip to content

Commit

Permalink
-
Browse files Browse the repository at this point in the history
  • Loading branch information
BB committed Jan 9, 2023
1 parent 42d4019 commit 9a3f7e6
Show file tree
Hide file tree
Showing 87 changed files with 689 additions and 689 deletions.
30 changes: 15 additions & 15 deletions Intro_Solid_Mechanics_Adeeb/10.1.3.1.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,14 @@
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)]
display("\u03C1b_1 =",pb1)
display("\u03C1b_2 =",pb2)
display("\u03C1b_3 =",pb3)
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])
display("u* =",ustar)
print("u* =",ustar)
gradustar = Matrix([[diff(i,j) for j in x[:2]]for i in ustar])
display("\u2207u* =",gradustar)
print("\u2207u* =",gradustar)
ustara = ustar.subs({x1:2})
ustarb = ustar.subs({x2:1})
ustarc = ustar.subs({x1:0})
Expand All @@ -25,22 +25,22 @@
tb = (s*nb).subs({x2:1})
tc = (s*nc).subs({x1:0})
td = (s*nd).subs({x2:0})
display("t_A =",ta,"t_B =",tb,"t_C =",tc,"t_D =",td)
print("t_A =",ta,"t_B =",tb,"t_C =",tc,"t_D =",td)
estar = sp.Rational("1/2")*(gradustar+gradustar.T)
display("\u03B5* =",estar)
print("\u03B5* =",estar)
ststrain = sum([s[i]*estar[i] for i in range(4)])
display("ststrain =",ststrain)
print("ststrain =",ststrain)
IVW = integrate(ststrain,(x1,0,2),(x2,0,1),(x3,0,t))
display("Internal Virtual Work =",IVW)
print("Internal Virtual Work =",IVW)
EVWBodyForces = integrate(pb1*ustar[0],(x1,0,2),(x2,0,1),(x3,0,t))
display("External Virtual Work|_{Body Forces} =",EVWBodyForces)
print("External Virtual Work|_{Body Forces} =",EVWBodyForces)
EVWa = integrate((ta.dot(ustara)).subs({x1:2}),(x2,0,1),(x3,0,t))
display("External Virtual Work|_A =",EVWa)
print("External Virtual Work|_A =",EVWa)
EVWb = integrate((tb.dot(ustarb)).subs({x2:1}),(x1,0,2),(x3,0,t))
display("External Virtual Work|_B =",EVWb)
print("External Virtual Work|_B =",EVWb)
EVWc = integrate((tc.dot(ustarc)).subs({x1:0}),(x2,0,1),(x3,0,t))
display("External Virtual Work|_C =",EVWc)
print("External Virtual Work|_C =",EVWc)
EVWd = integrate((td.dot(ustard)).subs({x2:0}),(x1,0,2),(x3,0,t))
display("External Virtual Work|_D =",EVWd)
print("External Virtual Work|_D =",EVWd)
EVW = simplify(EVWBodyForces+EVWa+EVWb+EVWc+EVWd)
display("External Virtual Work = Internal Virtual Work =",EVW)
print("External Virtual Work = Internal Virtual Work =",EVW)
16 changes: 8 additions & 8 deletions Intro_Solid_Mechanics_Adeeb/10.1.3.2.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,29 +10,29 @@
q = -5*x
s = dsolve(EI*y(x).diff(x,4)-q,y(x),ics={y1:0,y2:0,yp1:0,yp2:0})
y = s.rhs
display("y(x) =",y)
print("y(x) =",y)
th = simplify(diff(y,x))
M = simplify(EI*diff(y,x,x))
V = simplify(EI*diff(y,x,x,x))
display("th =",th,"M =",M,"V =",V)
print("th =",th,"M =",M,"V =",V)
M1 = M.subs({x:0})
M2 = M.subs({x:L})
V1 = V.subs({x:0})
V2 = V.subs({x:L})
display("M@x=0",M1,"M@x=L",M2,"V@x=0",V1,"V@x=L",V2)
print("M@x=0",M1,"M@x=L",M2,"V@x=0",V1,"V@x=L",V2)
ystar = a*x**2
thstar = diff(ystar,x)
ystar1 = ystar.subs({x:0})
ystar2 = ystar.subs({x:L})
thstar1 = thstar.subs({x:0})
thstar2 = thstar.subs({x:L})
display("y*@x=0",ystar1,"y*@x=L",ystar2,"\u03B8*@x=0",thstar1,"\u03B8*@x=L",thstar2)
print("y*@x=0",ystar1,"y*@x=L",ystar2,"\u03B8*@x=0",thstar1,"\u03B8*@x=L",thstar2)
IVW = integrate(M*diff(ystar,x,x),(x,0,L))
display("Internal Virtual Work =",IVW)
print("Internal Virtual Work =",IVW)
EVWq = integrate(q*ystar,(x,0,L))
display("EVWq =",EVWq)
print("EVWq =",EVWq)
EVW1 = V1*ystar1-M1*thstar1
EVW2 = M2*thstar2-V2*ystar2
display("EVW1 =",EVW1,"EVW2 =",EVW2)
print("EVW1 =",EVW1,"EVW2 =",EVW2)
EVW = EVWq+EVW1+EVW2
display("External Virtual Work = Internal Virtual Work =",EVW)
print("External Virtual Work = Internal Virtual Work =",EVW)
26 changes: 13 additions & 13 deletions Intro_Solid_Mechanics_Adeeb/10.2.3.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,19 @@
s = solve(Equil,th)
s=[(i.expand(complex=True)).evalf() for i in s]
degs = [(i/sp.pi*180).evalf() for i in s]
display("Equilibrium =",Equil)
display("there are 4 real roots of the above expression")
display(s)
display("\u03B8_1 =",degs[1])
display("\u03B8_2 =",degs[2])
display("\u03B8_3 =",degs[3])
display("\u03B8_4 =",degs[0]+360)
print("Equilibrium =",Equil)
print("there are 4 real roots of the above expression")
print(s)
print("\u03B8_1 =",degs[1])
print("\u03B8_2 =",degs[2])
print("\u03B8_3 =",degs[3])
print("\u03B8_4 =",degs[0]+360)
PE = 750/2*(2*sin(th))**2-140*th-900*sin(th)
DPE = diff(PE,th)
D2PE = diff(DPE,th)
display("Potential Energy =",PE)
display("Derivative of PE",DPE,"Second Derivative of Pe",D2PE)
display("D2PE@\u03B8=\u03B8_1 =",D2PE.subs({th:s[1]}))
display("D2PE@\u03B8=\u03B8_2 =",D2PE.subs({th:s[2]}))
display("D2PE@\u03B8=\u03B8_3 =",D2PE.subs({th:s[3]}))
display("D2PE@\u03B8=\u03B8_4 =",D2PE.subs({th:s[0]+2*sp.pi}))
print("Potential Energy =",PE)
print("Derivative of PE",DPE,"Second Derivative of Pe",D2PE)
print("D2PE@\u03B8=\u03B8_1 =",D2PE.subs({th:s[1]}))
print("D2PE@\u03B8=\u03B8_2 =",D2PE.subs({th:s[2]}))
print("D2PE@\u03B8=\u03B8_3 =",D2PE.subs({th:s[3]}))
print("D2PE@\u03B8=\u03B8_4 =",D2PE.subs({th:s[0]+2*sp.pi}))
10 changes: 5 additions & 5 deletions Intro_Solid_Mechanics_Adeeb/11.1.1.1.b.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,13 @@
u = a2*x**2+a1*x
uL = u.subs(x,L)
PE = integrate((1/2)*E*A*(u.diff(x)**2), (x,0,L)) - P*uL - integrate(c*x*u, (x,0,L))
display("Potential Energy: ", PE)
print("Potential Energy: ", PE)
Eq1 = PE.diff(a2)
Eq2 = PE.diff(a1)
display("Minimize PE: ", Eq1, Eq2)
print("Minimize PE: ", Eq1, Eq2)
s = solve((Eq1, Eq2), (a2, a1))
display("Solve: ", s)
print("Solve: ", s)
u = u.subs({a1:s[a1], a2:s[a2]})
display("Best second degree Polynomial(Rayleigh Ritz method): ", u)
print("Best second degree Polynomial(Rayleigh Ritz method): ", u)
stress = E * u.diff(x)
display("stress: ", simplify(stress))
print("stress: ", simplify(stress))
10 changes: 5 additions & 5 deletions Intro_Solid_Mechanics_Adeeb/11.1.1.1.c.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,14 @@
u = a3*x**3+a2*x**2+a1*x
uL = u.subs(x,L)
PE = integrate((1/2)*E*A*(u.diff(x)**2), (x,0,L)) - P*uL - integrate(c*x*u, (x,0,L))
display("Potential Energy: ", PE)
print("Potential Energy: ", PE)
Eq1 = PE.diff(a2)
Eq2 = PE.diff(a1)
Eq3 = PE.diff(a3)
display("Minimize PE: ", Eq1, Eq2, Eq3)
print("Minimize PE: ", Eq1, Eq2, Eq3)
s = solve((Eq1, Eq2, Eq3), (a2, a1, a3))
display("Solve: ", s)
print("Solve: ", s)
u = u.subs({a1:s[a1], a2:s[a2], a3:s[a3]})
display("Best third degree Polynomial(Rayleigh Ritz method): ", u)
print("Best third degree Polynomial(Rayleigh Ritz method): ", u)
stress = E * u.diff(x)
display("stress: ", simplify(stress))
print("stress: ", simplify(stress))
4 changes: 2 additions & 2 deletions Intro_Solid_Mechanics_Adeeb/11.1.1.1.d.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,8 @@
u2 = u(x).diff(x).subs(x,L)
Ax = 25/100*(5/10-25/200*x)
Ax1 = Ax.subs(x,L)
display("Area: ", Ax)
print("Area: ", Ax)
s = dsolve(E*u(x).diff(x,2)*Ax+E*u(x).diff(x)*Ax.diff(x), u(x), ics = {u1:0, u2:200/E/Ax1})
u = s.rhs.subs({E:100000, L:2, P:200})
stress = (E * u.diff(x)).subs({E:100000, L:2, P:200})
display("displacement and stress: ", u, stress)
print("displacement and stress: ", u, stress)
44 changes: 22 additions & 22 deletions Intro_Solid_Mechanics_Adeeb/11.1.1.1.e.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,59 +6,59 @@
sp.init_printing(use_latex = "mathjax")
x, a1, a2, a3, E, P = symbols("x a_1 a_2 a_3 E P")
Ax = 25/100*(5/10-25/200*x)
display("Exact")
print("Exact")
u = Function("u")
u1_1 = u(x).subs(x,0)
u2_1 = u(x).diff(x).subs(x,L)
Ax1 = Ax.subs(x,L)
display("Area: ", Ax)
print("Area: ", Ax)
s = dsolve(E*u(x).diff(x,2)*Ax+E*u(x).diff(x)*Ax.diff(x), u(x), ics = {u1_1:0, u2_1:200/E/Ax1})
u = s.rhs.subs({E:100000, L:2, P:200})
stress0 = (E * u.diff(x)).subs({E:100000, L:2, P:200})
display("displacement and stress: ", u, stress0)
display("First Degree")
print("displacement and stress: ", u, stress0)
print("First Degree")
u1 = a1*x
uL = u1.subs(x,L)
PE = integrate((1/2)*E*((u1.diff(x))**2)*Ax, (x,0,L)) - P*uL
PE = PE.subs({E:100000, L:2, P:200})
display("Potential Energy: ", PE)
print("Potential Energy: ", PE)
Eq1 = PE.diff(a1)
display("Minimize:", Eq1)
print("Minimize:", Eq1)
s = solve(Eq1,a1)
display("Solve: ", s)
print("Solve: ", s)
u1 = u1.subs(a1,s[0])
display("Best First degree Polynomial(Rayleigh Ritz method): ", u1)
print("Best First degree Polynomial(Rayleigh Ritz method): ", u1)
stress1 = (E * u1.diff(x)).subs(E,100000)
display("stress: ", stress1)
display("Second Degree")
print("stress: ", stress1)
print("Second Degree")
u2 = a1*x+a2*x**2
PE2 = integrate((1/2)*E*((u2.diff(x))**2)*Ax, (x,0,L)) - P*u2.subs(x,L)
PE2 = PE2.subs({E:100000, L:2, P:200})
display("Potential Energy: ", PE2)
print("Potential Energy: ", PE2)
Eq1_2 = PE2.diff(a1)
Eq2_2 = PE2.diff(a2)
s = solve((Eq1_2, Eq2_2),a1, a2)
display("Minimize:", Eq1_2, Eq2_2)
display("Solve: ", s)
print("Minimize:", Eq1_2, Eq2_2)
print("Solve: ", s)
u2 = u2.subs({a1:s[a1], a2:s[a2]})
display("Best Second degree Polynomial(Rayleigh Ritz method): ", u2)
print("Best Second degree Polynomial(Rayleigh Ritz method): ", u2)
stress2 = (E * u2.diff(x)).subs(E,100000)
display("stress: ", stress2)
display("Third Degree")
print("stress: ", stress2)
print("Third Degree")
u3 = a1*x+a2*x**2+a3*x**3
PE3 = integrate((1/2)*E*((u3.diff(x))**2)*Ax, (x,0,L)) - P*u3.subs(x,L)
PE3 = PE3.subs({E:100000, L:2, P:200})
display("Potential Energy: ", PE3)
print("Potential Energy: ", PE3)
Eq1_3 = PE3.diff(a1)
Eq2_3 = PE3.diff(a2)
Eq3_3 = PE3.diff(a3)
s = solve((Eq1_3, Eq2_3, Eq3_3),a1, a2, a3)
display("Minimize:", Eq1_3, Eq2_3, Eq3_3)
display("Solve: ", s)
print("Minimize:", Eq1_3, Eq2_3, Eq3_3)
print("Solve: ", s)
u3 = u3.subs({a1:s[a1], a2:s[a2], a3:s[a3]})
display("Best Third degree Polynomial(Rayleigh Ritz method): ", u3)
print("Best Third degree Polynomial(Rayleigh Ritz method): ", u3)
stress3 = (E * u3.diff(x)).subs(E,100000)
display("stress: ", stress3)
print("stress: ", stress3)
fig, ax = plt.subplots(1,2, figsize = (15,6))
plt.setp(ax[0], xlabel = "X1 m ", ylabel = "u(m)")
plt.setp(ax[1], xlabel = "X1 m", ylabel = "Stress(Pa)")
Expand All @@ -71,7 +71,7 @@
stress1_list = [stress1.subs({x:i}) for i in x1]
stress2_list = [stress2.subs({x:i}) for i in x1]
stress3_list = [stress3.subs({x:i}) for i in x1]
display("Comparison")
print("Comparison")
ax[0].plot(x1, u_list, 'blue', label = "exact")
ax[0].plot(x1, u1_list, 'orange', label = "u1")
ax[0].plot(x1, u2_list, 'green', label = "u2")
Expand Down
48 changes: 24 additions & 24 deletions Intro_Solid_Mechanics_Adeeb/11.1.1.1.f.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,65 +6,65 @@
sp.init_printing(use_latex = "mathjax")
x, a1, a2, a3, a4, E, p, A, L= symbols("x a_1 a_2 a_3 a_4 E p A L")
p = 5*x**2
display("Exact")
print("Exact")
u = Function("u")
u1_1 = u(x).subs(x,0)
u2_1 = u(x).subs(x,L)
s = dsolve(E*u(x).diff(x,2)+ p/A, u(x), ics = {u1_1:0, u2_1:0})
u = s.rhs.subs({E:100000, L:2, A:(25/100)**2})
stress0 = (E * u.diff(x)).subs({E:100000, L:2, A:(25/100)**2})
display("displacement and stress: ", u, stress0)
display("Second Degree")
print("displacement and stress: ", u, stress0)
print("Second Degree")
u2 = a1*x+a2*x**2
sol = solve(u2.subs(x,L), a2)
u2 = u2.subs(a2,sol[0])
display("a2:", sol)
display(u2)
print("a2:", sol)
print(u2)
PE2 = (1/2)*(integrate(E*A*((u2.diff(x))**2), (x,0,L))) - integrate(p*u2, (x,0,L))
PE2 = PE2.subs({E:100000, L:2, A:(25/100)**2})
display("Potential Energy: ", PE2)
print("Potential Energy: ", PE2)
Eq1_2 = PE2.diff(a1)
sol1 = solve(Eq1_2,a1)
display("a1:", sol1)
print("a1:", sol1)
u2 = u2.subs(a1,sol1[0]).subs(L,2)
display("Best Second degree Polynomial(Rayleigh Ritz method): ", u2)
print("Best Second degree Polynomial(Rayleigh Ritz method): ", u2)
stress2 = (E * u2.diff(x)).subs(E,100000)
display("stress: ", stress2)
display("Third Degree")
print("stress: ", stress2)
print("Third Degree")
u3 = a1*x+a2*x**2+a3*x**3
sol = solve(u3.subs(x,L), a3)
display("a3: ", sol)
print("a3: ", sol)
u3 = u3.subs(a3,sol[0])
PE3 = (1/2)*(integrate(E*A*((u3.diff(x))**2), (x,0,L))) - integrate(p*u3, (x,0,L))
PE3 = PE3.subs({E:100000, L:2, A:(25/100)**2})
display("Potential Energy: ", PE3)
print("Potential Energy: ", PE3)
Eq1_3 = PE3.diff(a1)
Eq2_3 = PE3.diff(a2)
sol1 = solve((Eq1_3, Eq2_3),a1, a2)
display("Minimize:", Eq1_3, Eq2_3)
display("Solve: ", sol1)
print("Minimize:", Eq1_3, Eq2_3)
print("Solve: ", sol1)
u3 = u3.subs({a1:sol1[a1], a2:sol1[a2], L:2})
display("Best Third degree Polynomial(Rayleigh Ritz method): ", u3)
print("Best Third degree Polynomial(Rayleigh Ritz method): ", u3)
stress3 = (E * u3.diff(x)).subs(E,100000)
display("stress: ", stress3)
display("Fourth Degree")
print("stress: ", stress3)
print("Fourth Degree")
u4 = a1*x+a2*x**2+a3*x**3+a4*x**4
sol = solve(u4.subs(x,L), a4)
display("a4: ", sol)
print("a4: ", sol)
u4 = u4.subs(a4,sol[0])
PE4 = (1/2)*(integrate(E*A*((u4.diff(x))**2), (x,0,L))) - integrate(p*u4, (x,0,L))
PE4 = PE4.subs({E:100000, L:2, A:(25/100)**2})
display("Potential Energy: ", PE4)
print("Potential Energy: ", PE4)
Eq1_4 = PE4.diff(a1)
Eq2_4 = PE4.diff(a2)
Eq3_4 = PE4.diff(a3)
sol1 = solve((Eq1_4, Eq2_4, Eq3_4),(a1, a2, a3))
display("Minimize:", Eq1_4, Eq2_4, Eq3_4)
display("Solve: ", sol1)
print("Minimize:", Eq1_4, Eq2_4, Eq3_4)
print("Solve: ", sol1)
u4 = u4.subs({a1:sol1[a1], a2:sol1[a2], a3:sol1[a3]}).subs(L,2)
display("Best Fourth degree Polynomial(Rayleigh Ritz method): ", u3)
print("Best Fourth degree Polynomial(Rayleigh Ritz method): ", u3)
stress4 = (E * u4.diff(x)).subs(E,100000)
display("stress: ", stress4)
print("stress: ", stress4)
fig, ax = plt.subplots(1,2, figsize = (15,6))
plt.setp(ax[0], xlabel = "X1 m ", ylabel = "u(m)")
plt.setp(ax[1], xlabel = "X1 m", ylabel = "Stress(Pa)")
Expand All @@ -77,7 +77,7 @@
stress2_list = [stress2.subs({x:i}) for i in x1]
stress3_list = [stress3.subs({x:i}) for i in x1]
stress4_list = [stress4.subs({x:i}) for i in x1]
display("Comparison")
print("Comparison")
ax[0].plot(x1, u_list, 'blue', label = "exact")
ax[0].plot(x1, u2_list, 'orange', label = "u2")
ax[0].plot(x1, u3_list, 'green', label = "u3")
Expand Down
2 changes: 1 addition & 1 deletion Intro_Solid_Mechanics_Adeeb/11.1.1.1.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,4 +8,4 @@
a = dsolve(u(x).diff(x,2) + C*x/(E*A), u(x), ics = {u1:0, u2:P/(E*A)})
u = a.rhs
sigma = E * u.diff(x)
display("displacement and stress: ", u, simplify(sigma))
print("displacement and stress: ", u, simplify(sigma))
Loading

0 comments on commit 9a3f7e6

Please sign in to comment.