Plasticity extension elements statically indeterminate#
Click –> Live Code on the top right corner of this screen to investigate some plasticity!
import micropip
await micropip.install("ipympl")
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import numpy as np
import sympy as sym
import ipywidgets as widgets
#from IPython.display import display
%matplotlib widget
wM = sym.symbols('wM',real=True,positive=True)
deltawAM = sym.symbols('deltawAM',real=True)
E = sym.nsimplify(200e9)
A1 = sym.nsimplify(10e-6)
fy = sym.nsimplify(500e6)
bm = sym.S(2)
L = sym.S(6)
kA = E * A1 / 6
kB = E * A1 / 10
kC = E * A1 / 8 #(op x = 3)
wAz = wM-deltawAM/3*(6-bm)
wBz = wM+deltawAM/3*bm
wBx = sym.cos(deltawAM/(L/(bm/6)))*6
wCz = wM-deltawAM/3*(3-bm)
wCx = sym.cos(deltawAM/(L/(bm/6)))*3
FA = kA * wAz
FC = kC * wCz
FB = kB * wBz
F = sym.symbols('F',real=True,positive=True)
eq1 = sym.Eq(FA + FB + FC- F,0)
eq2 = sym.Eq(F * 4 - FC * 3 - FB * 6,0)
F_sol = sym.solve(eq1,F)[0]
deltawAM_sol_1 = sym.solve(eq2.subs(F,F_sol),deltawAM)[0]
wyA = sym.solve(sym.Eq(FA.subs(deltawAM,deltawAM_sol_1),fy*A1),wM)[0]
wyB = sym.solve(sym.Eq(FB.subs(deltawAM,deltawAM_sol_1),fy*A1),wM)[0]
wyC = sym.solve(sym.Eq(FC.subs(deltawAM,deltawAM_sol_1),fy*A1),wM)[0]
FB = fy*A1
eq1 = sym.Eq(FA + FB + FC- F,0)
eq2 = sym.Eq(F * 4 - FC * 3 - FB * 6,0)
F_sol = sym.solve(eq1,F)[0]
deltawAM_sol_2 = sym.solve(eq2.subs(F,F_sol),deltawAM)[0]
wyA = sym.solve(sym.Eq(FA.subs(deltawAM,deltawAM_sol_2),-fy*A1),wM)[0]
wyC = sym.solve(sym.Eq(FC.subs(deltawAM,deltawAM_sol_2),fy*A1),wM)[0]
FC = fy*A1
eq1 = sym.Eq(FA + FB + FC- F,0)
eq2 = sym.Eq(F * 4 - FC * 3 - FB * 6,0)
F_sol = sym.solve(eq1,F)[0]
deltawAM_sol_3 = sym.solve(eq2.subs(F,F_sol),deltawAM)[0]
deltawAM_sol = sym.Piecewise((deltawAM_sol_1, wM < wyB), (deltawAM_sol_2, wM < wyC), (deltawAM_sol_3, True))
FB = sym.Piecewise((kB * wBz,wM < wyB), (fy*A1, True))
FC = sym.Piecewise((kC * wCz,wM < wyC), (fy*A1, True))
eq1 = sym.Eq(FA + FB + FC- F,0)
F_sol = sym.solve(eq1,F)[0]
eps = sym.symbols('eps',real=True,positive=True)
sigmaA = FA.subs(deltawAM,deltawAM_sol)/A1
epsA = wAz.subs(deltawAM,deltawAM_sol)/6
sigmaB = FB.subs(deltawAM,deltawAM_sol)/A1
epsB = wBz.subs(deltawAM,deltawAM_sol)/10
sigmaC = FC.subs(deltawAM,deltawAM_sol)/A1
epsC = wCz.subs(deltawAM,deltawAM_sol)/8
wAz_nump = sym.lambdify(wM,wAz.subs(deltawAM,deltawAM_sol))
wBz_nump = sym.lambdify(wM,wBz.subs(deltawAM,deltawAM_sol))
wBx_nump = sym.lambdify(wM,wBx.subs(deltawAM,deltawAM_sol))
wCz_nump = sym.lambdify(wM,wCz.subs(deltawAM,deltawAM_sol))
wCx_nump = sym.lambdify(wM,wCx.subs(deltawAM,deltawAM_sol))
sigmaA_nump = sym.lambdify(wM,sigmaA)
sigmaB_nump = sym.lambdify(wM,sigmaB)
sigmaC_nump = sym.lambdify(wM,sigmaC)
epsA_nump = sym.lambdify(wM,epsA)
epsB_nump = sym.lambdify(wM,epsB)
epsC_nump = sym.lambdify(wM,epsC)
F_nump = sym.lambdify(wM,F_sol.subs(deltawAM,deltawAM_sol))
fig = plt.figure(figsize=(10,5))
ax1 = fig.add_subplot(1,3,1)
ax1.set_aspect('equal')
ax2 = fig.add_subplot(2,3,2)
ax3 = fig.add_subplot(2,3,3)
ax4 = fig.add_subplot(2,3,5)
ax5 = fig.add_subplot(2,3,6)
ax2.spines['right'].set_color('none')
ax2.spines['top'].set_color('none')
ax3.spines['right'].set_color('none')
ax3.spines['top'].set_color('none')
ax4.spines['right'].set_color('none')
ax4.spines['top'].set_color('none')
ax5.spines['right'].set_color('none')
ax5.spines['top'].set_color('none')
fig.tight_layout()
scale = 30
w_slider = widgets.FloatSlider(value=0, min=0, max=0.2, step=0.0001, description='w:')
wM_list = np.array([0,wyB,wyC,0.2])
def update_plot(w):
wM_num = w
ax1.clear() # Clear the existing plot
ax2.clear()
ax3.clear()
ax4.clear()
ax5.clear()
ax1.set_xlim([-1,7])
ax1.set_ylim([-7,11])
ax1.plot([0,0],[-wAz_nump(wM_num)*scale,6])
# Get the color of the previous plot
color1 = ax1.get_lines()[-1].get_color()
if wM_num < wyB:
ax1.plot([wBx_nump(wM_num),wBx_nump(wM_num)],[-wBz_nump(wM_num)*scale,10])
else:
ax1.plot([wBx_nump(wM_num),wBx_nump(wM_num)],[-wBz_nump(wM_num)*scale,10],linewidth=3)
color2 = ax1.get_lines()[-1].get_color()
ax1.plot([0,wBx_nump(wM_num)],[-wAz_nump(wM_num)*scale,-wBz_nump(wM_num)*scale],color='black')
if wM_num < wyC:
ax1.plot([wCx_nump(wM_num),wCx_nump(wM_num)],[-wCz_nump(wM_num)*scale,8])
else:
ax1.plot([wCx_nump(wM_num),wCx_nump(wM_num)],[-wCz_nump(wM_num)*scale,8],linewidth=3)
color3 = ax1.get_lines()[-1].get_color()
ax1.annotate('', xy=(wBx_nump(wM_num)*((6-bm)/6), -w*scale), xytext=(wBx_nump(wM_num)*((6-bm)/6), 0) , arrowprops=dict(arrowstyle='fancy'))
ax1.text(3,9.5,'Vertical deformation\nscaled by 30',fontsize=12,ha='center',va='center')
ax1.axis('off')
y1 = np.array([[0-1/3,6+1/2], [0+0,6+0], [0+1/3,6+1/2]])
p1 = Polygon(y1, facecolor = '#00CC99',edgecolor='black')
ax1.add_patch(p1)
ax1.plot([0-1/2,0+1/2],[6+1/2,6+1/2],color='black')
y2 = np.array([[wBx_nump(wM_num)-1/3,10+1/2], [wBx_nump(wM_num)+0,10+0], [wBx_nump(wM_num)+1/3,10+1/2]])
p2 = Polygon(y2, facecolor = '#00CC99',edgecolor='black')
ax1.add_patch(p2)
ax1.plot([wBx_nump(wM_num)-1/2,wBx_nump(wM_num)+1/2],[10+1/2,10+1/2],color='black')
ax1.plot([wBx_nump(wM_num)-1/2,wBx_nump(wM_num)+1/2],[10+1/2+1/4,10+1/2+1/4],color='black')
y3 = np.array([[wCx_nump(wM_num)-1/3,8+1/2], [wCx_nump(wM_num)+0,8+0], [wCx_nump(wM_num)+1/3,8+1/2]])
p3 = Polygon(y3, facecolor = '#00CC99',edgecolor='black')
ax1.add_patch(p3)
ax1.plot([wCx_nump(wM_num)-1/2,wCx_nump(wM_num)+1/2],[8+1/2,8+1/2],color='black')
ax1.plot([wCx_nump(wM_num)-1/2,wCx_nump(wM_num)+1/2],[8+1/2+1/4,8+1/2+1/4],color='black')
ax3.plot(epsB_nump(wM_list),sigmaB_nump(wM_list)/1e6,color=color2)
ax3.plot(epsB_nump(wM_num),sigmaB_nump(wM_num)/1e6,'o',color=color2)
ax3.set_ylabel('$\sigma$ [MPa]')
ax3.set_xlabel('$\epsilon$ [MPa]')
ax3.set_xlim([-0.001,max(epsB_nump(0.1)+0.001,epsB_nump(wM_num)+0.001)])
ax2.plot(epsC_nump(wM_list),sigmaC_nump(wM_list)/1e6,color=color3)
ax2.plot(epsC_nump(wM_num),sigmaC_nump(wM_num)/1e6,'o',color=color3)
ax2.set_ylabel('$\sigma$ [MPa]')
ax2.set_xlabel('$\epsilon$ [MPa]')
ax2.set_xlim([-0.001,max(epsC_nump(0.1)+0.001,epsC_nump(wM_num)+0.001)])
ax4.plot(epsA_nump(wM_list),sigmaA_nump(wM_list)/1e6,color=color1)
ax4.plot(epsA_nump(wM_num),sigmaA_nump(wM_num)/1e6,'o',color=color1)
ax4.set_ylabel('$\sigma$ [MPa]')
ax4.set_xlabel('$\epsilon$ [MPa]')
ax5.plot(wM_list,F_nump(wM_list),color='black')
ax5.plot(wM_num,F_nump(wM_num),'o',color='black')
ax5.set_ylabel('$F$ [N]')
ax5.set_xlabel('$w$ [m]')
ax5.set_xlim([0-0.002,max(0.1+0.002,wM_num+0.002)])
ax2.set_yticklabels([])
ax2.set_xticklabels([])
ax3.set_yticklabels([])
ax3.set_xticklabels([])
ax4.set_yticklabels([])
ax4.set_xticklabels([])
ax5.set_yticklabels([])
ax5.set_xticklabels([])
plt.draw()
widgets.interact(update_plot, w = w_slider);