如何创建以特定方式定位的绘图功能（门户框架 - 力矩图）？
How to create a plot functions orientated in a specific way (portal frame - moment diagram)?
我是一名机械工程专业的学生，我正在尝试为 Python 中的力矩图创建一个脚本。我的代码中缺少的是如何定位力矩函数以便像门户框架一样对齐。
Mpilar1 是第一列（从左到右）的矩函数。
Masna1 是第一个梁的力矩函数（从左到右）。
Masna2 是第二个梁（从左到右）的力矩函数。
Mpilar2 是第二列（从左到右）的矩函数。
代码：
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
alfa = (mt.atan((int(ht)-int(h))/(int(v)/2)))*180/((mt.pi))
print("Ângulo da vertente:", round(alfa, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
def draw_line():
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.show()
if __name__ == '__main__':
draw_line()
Mpilar1 = 1500 * h1 ** 2 + 350 * h1
Masna1 = 300 * ha1 ** 2 + 15 * ha1
Masna2 = 200 * ha1 ** 2 + 15 * ha1
Mpilar2 = 1400 * h1 ** 2 + 10 * h1
plt.plot(h1, Mpilar1)
plt.plot(ha1, Masna1)
plt.plot(ha1, Masna2)
plt.plot(h1, Mpilar2)
您必须在代表曲线的点上使用变换矩阵。特别是，您必须使用 roto-translation 矩阵将曲线旋转并平移到正确的位置和方向，并且您可能必须应用镜像矩阵以使力矩根据您的约定对齐。
请注意，我在结构工程领域的日子已经一去不复返了，所以我不记得正确定位时刻的惯例。留给你作为练习。
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
alfa_rad = mt.atan((int(ht)-int(h))/(int(v)/2))
alfa_deg = alfa_rad*180/mt.pi
print("Ângulo da vertente:", round(alfa_deg, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
# Roto-translation matrix:
# Rotates the points by an angle theta and translates
# them by x in the horizontal direction, and y in the
# vertical direction
R = lambda x, y, theta: np.array([
[np.cos(theta), np.sin(theta), x],
[-np.sin(theta), np.cos(theta), y],
[0, 0, 1],
])
# mirror matrix about the x-axis
Mx = np.array([
[1, 0, 0], [0, -1, 0], [0, 0, 1]
])
# mirror matrix about the y-axis
My = np.array([
[-1, 0, 0], [0, 1, 0], [0, 0, 1]
])
Mpilar1 = 1500 * h1 ** 2 + 350 * h1
Masna1 = 300 * ha1 ** 2 + 15 * ha1
Masna2 = 200 * ha1 ** 2 + 15 * ha1
Mpilar2 = 1400 * h1 ** 2 + 10 * h1
def draw_line():
plt.figure()
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
# left column
points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])
points1 = np.matmul(R(0, 0, -np.pi/2), points1)
plt.plot(points1[0, :], points1[1, :], label="Mpilar1")
# right column
points2 = np.stack([h1, Mpilar2 / max(Mpilar2), np.ones_like(h1)])
points2 = np.matmul(R(20, 0, -np.pi/2), points2)
plt.plot(points2[0, :], points2[1, :], label="Mpilar2")
# left asna
points3 = np.stack([ha1, Masna1 / max(Masna1), np.ones_like(ha1)])
points3 = np.matmul(R(0, 6, -alfa_rad), points3)
plt.plot(points3[0, :], points3[1, :], label="Masna1")
# right asna
points4 = np.stack([ha1, Masna2 / max(Masna2), np.ones_like(ha1)])
points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)
plt.plot(points4[0, :], points4[1, :], label="Masna2")
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.legend()
plt.show()
draw_line()
上面的代码有几点需要注意：
让我们考虑一下points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])。它创建一个 3xn 坐标矩阵。第一行是x-coordinates、h1。第二行是当前的 y-coordinates，Mpilar1 / max(Mpilar1)（请注意，为了适合图表，我已经对其进行了维度化）。第三行是1，能够应用平移矩阵是一个技巧。在绘图命令中，我们将只使用第一行和第二行（x 和 y 坐标）。
points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)这里我先镜像了关于y-axis的点，然后我应用了旋转和平移。你必须玩才能正确定位时刻！
解决方案有效，但在这种情况下，我没有正确的力矩图。
代码：
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
nm = 7 #("Número de madres por asna:")
npo = 5 #("Número de pórticos:")
dp = 5 #("Distância entre pórticos em metros:")
alfa_rad = mt.atan((int(ht)-int(h))/(int(v)/2))
alfa_deg = alfa_rad*180/mt.pi
alfa = (mt.atan((int(ht)-int(h))/(int(v)/2)))*180/((mt.pi))
print("Ângulo da vertente:", round(alfa, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
dm = lasna / nm
print(("Distância entre madres da cobertura em metros:"), round(dm, 2), "m")
lp = npo * dp
print("Comprimento total do pavilhão:", round(lp, 1), "m")
def draw_line():
# x axis value list.
x_number_list = [0, 0, (v/2), v, v]
# y axis value list.
y_number_list = [0, h, ht, h, 0]
# Plot the number in the list and set the line thickness.
plt.plot(x_number_list, y_number_list, linewidth=3)
# Set the line chart title and the text font size.
plt.title("Pórtico", fontsize=15)
# Set x axis label.
plt.xlabel("Vão (m)", fontsize=10)
# Set y axis label.
plt.ylabel("Altura (m)", fontsize=10)
# Set the x, y axis tick marks text size.
plt.tick_params(axis='both', labelsize=9)
# Display the plot in the matplotlib's viewer.
plt.show()
if __name__ == '__main__':
draw_line()
#EXEMPLO - IPE 100
pppilar = 8.1 #variavel do peso próprio do pilar
arpilar = 10.32 * 10 ** -4 #variavél da área de secção do pilar - m^2
inypilar = 171 * 10 ** -8 #variavél da inércia pilar segundo y - m^4
ppasna = 8.1 #variavel do peso próprio da asna
arasna = 10.32 * 10 ** -4 #variavél da área de secção do pilar - m^2
inyasna = 171 * 10 ** -8 #variavél da inércia pilar - m^4
pesomadre = 4.84 #(Peso linear das madres de cobertura em kg/m)
pesopainel = 11.2 #(Peso painel revestimento kg/m2)
rotulado = 1
base = rotulado
#Definir ações
ppcobertura = (pesopainel * dp + ((pesomadre*dp)/dm) + ppasna) * 9.81 / 1000 #kN/m
#print("Peso próprio da cobertura:", round(ppcobertura, 1), "kN/m")
#Sobrecarga
qk = 0.4 #Sobrecarga em kN/m2
sb = qk * dp #Sobrecarga em kN/m
#Neve
Sk = 0.5 #Neve em kN/m2
ne = Sk * dp #Neve em kN/m
#Vento
Qp1 = 0.5 #Vento pilar 1 em kN/m2
vnp1 = Qp1 * dp #Vento pilar 1 em kN/m
Qp2 = 0.5 #Vento pilar 2 em kN/m2
vnp2 = Qp2 * dp #Vento pilar 2 em kN/m
Qa1 = 0.5 #Vento asna 1 em kN/m2
vna1 = Qa1 * dp #Vento asna 1 em kN/m
Qa2 = 0.5 #Vento asna 2 em kN/m2
vna2 = Qa2 * dp #Vento asna 2 em kN/m
#Decompor as ações em normal e tagencial
ppcoberturan = ppcobertura * mt.cos(alfa*mt.pi/180)
ppcoberturat = ppcobertura * mt.sin(alfa*mt.pi/180)
sbn = sb * mt.cos(alfa*mt.pi/180) * mt.cos(alfa*mt.pi/180)
sbt = sb * mt.cos(alfa*mt.pi/180) * mt.sin(alfa*mt.pi/180)
nen = ne * mt.cos(alfa*mt.pi/180) * mt.cos(alfa*mt.pi/180)
net = ne * mt.cos(alfa*mt.pi/180) * mt.sin(alfa*mt.pi/180)
#Coeficientes de majoração
psipp = 1.35
#print("\u03A8 peso próprio:", psipp)
psi0sb = 0
#print("\u03A8 0 sobrecarga:", psi0sb)
psi1sb = 0
#print("\u03A8 1 sobrecarga:", psi1sb)
psi2sb = 0
#print("\u03A8 2 sobrecarga:", psi2sb)
ne1 = 1
ne2 = 2
nete = ne1
if nete == ne1:
psi0ne = 0.5
#print("\u03A8 0 neve:", psi0ne)
psi1ne = 0.2
#print("\u03A8 1 neve:", psi1ne)
psi2ne = 0
#print("\u03A8 2 neve:", psi2ne)
psi0vn = 0.6
#print("\u03A8 0 vento:", psi0vn)
psi1vn = 0.2
#print("\u03A8 1 vento:", psi1vn)
psi2vn = 0
#print("\u03A8 2 vento:", psi2vn)
#Combinação das ações para a cobertura - ELU - asna 1 - normal
comb_sbn = (psipp * ppcoberturan + sbn * 1.5 + (1.5 * psi0ne * nen + 1.5 * psi0vn * vna1)) * 1000 #N/m
comb_vnn = (psipp * ppcoberturan + vna1 * 1.5 + (1.5 * psi0ne * nen + 1.5 * sbn * psi0sb)) * 1000 #N/m
comb_nen =( psipp * ppcoberturan + nen * 1.5 + (1.5 * psi0vn * vna1 + 1.5 * sbn * psi0sb)) * 1000 #N/m
if (comb_sbn >= comb_vnn) and (comb_sbn >= comb_nen):
comb_a1n = comb_sbn
elif (comb_vnn >= comb_sbn) and (comb_vnn >= comb_nen):
comb_a1n = comb_vnn
else:
comb_a1n = comb_nen
#Combinação das ações para a cobertura - ELU - asna 1 - tangencial
comb_sbt = (psipp * ppcoberturat + sbt * 1.5 + (1.5 * psi0ne * net + 1.5 * psi0vn * 0)) * 1000 #N/m
comb_vnt = (psipp * ppcoberturat + 0 * 1.5 + (1.5 * psi0ne * net + 1.5 * sbt * psi0sb)) * 1000 #N/m
comb_net = (psipp * ppcoberturat + net * 1.5 + (1.5 * psi0vn * 0 + 1.5 * sbt * psi0sb)) * 1000 #N/m
if (comb_sbn >= comb_vnn) and (comb_sbn >= comb_nen):
comb_a1t = comb_sbt
elif (comb_vnn >= comb_sbn) and (comb_vnn >= comb_nen):
comb_a1t = comb_vnt
else:
comb_a1t = comb_net
#Combinação das ações para a cobertura - ELU - asna 2 - normal
comb_sb2n = (psipp * ppcoberturan + sbn * 1.5 + (1.5 * psi0ne * nen + 1.5 * psi0vn * vna2)) * 1000 #N/m
comb_vn2n = (psipp * ppcoberturan + vna2 * 1.5 + (1.5 * psi0ne * nen + 1.5 * sbn * psi0sb)) * 1000 #N/m
comb_ne2n = (psipp * ppcoberturan + nen * 1.5 + (1.5 * psi0vn * vna2 + 1.5 * sbn * psi0sb)) * 1000 #N/m
if (comb_sb2n >= comb_vn2n) and (comb_sb2n >= comb_ne2n):
comb_a2n = comb_sb2n
elif (comb_vn2n >= comb_sb2n) and (comb_vn2n >= comb_ne2n):
comb_a2n = comb_vn2n
else:
comb_a2n = comb_ne2n
#Combinação das ações para a cobertura - ELU - asna 2 - tangencial
comb_sbt2 = (psipp * ppcoberturat + sbt * 1.5 + (1.5 * psi0ne * net + 1.5 * psi0vn * 0)) * 1000 #N/m
comb_vnt2 = (psipp * ppcoberturat + 0 * 1.5 + (1.5 * psi0ne * net + 1.5 * sbt * psi0sb)) * 1000 #N/m
comb_net2 = (psipp * ppcoberturat + net * 1.5 + (1.5 * psi0vn * 0 + 1.5 * sbt * psi0sb)) * 1000 #N/m
if (comb_sb2n >= comb_vn2n) and (comb_sb2n >= comb_ne2n):
comb_a2t = comb_sbt2
elif (comb_vn2n >= comb_sb2n) and (comb_vn2n >= comb_ne2n):
comb_a2t = comb_vnt2
else:
comb_a2t = comb_net2
#Elementos finitos - Reações e deslocamentos
E = 210 * 10 ** 9 #módulo de elasticidade do aço em Pa
#Elemento 1 - asna1
a1 = E * arpilar / h
b1 = 12 * E * inypilar / h**3
c1 = 6 * E * inypilar / h**2
d1 = 4 * E * inypilar / h
e1 = 2 * E * inypilar / h
alfa1 = 90 * mt.pi / 180
l1 = mt.cos(alfa1)
m1 = mt.sin(alfa1)
t1 = np.matrix([[l1, m1, 0, 0, 0, 0],
[-m1, l1, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l1, m1, 0],
[0, 0, 0, -m1, l1, 0],
[0 , 0, 0, 0, 0, 1]])
k1local = np.matrix([[a1, 0, 0, -a1, 0, 0],
[0, b1, c1, 0, -b1, c1],
[0, c1, d1, 0, -c1, e1],
[-a1, 0, 0, a1, 0, 0],
[0, -b1, -c1, 0, b1, -c1],
[0, c1, e1, 0, -c1, d1]])
invt1 = np.matrix.transpose(t1)
k1global = np.linalg.multi_dot([invt1, k1local, t1])
#Elmento 2 - asna 1
a2 = E * arasna / lasna
b2 = 12 * E *inyasna / lasna**3
c2 = 6 * E * inyasna / lasna**2
d2 = 4 * E * inyasna / lasna
e2 = 2 * E *inyasna / lasna
alfa2 = ((alfa) * mt.pi) / 180
l2 = mt.cos(alfa2)
m2 = mt.sin(alfa2)
t2 = np.matrix([[l2, m2, 0, 0, 0, 0],
[-m2, l2, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l2, m2, 0],
[0, 0, 0, -m2, l2, 0],
[0 , 0, 0, 0, 0, 1]])
k2local = np.matrix([[a2, 0, 0, -a2, 0, 0],
[0, b2, c2, 0, -b2, c2],
[0, c2, d2, 0, -c2, e2],
[-a2, 0, 0, a2, 0, 0],
[0, -b2, -c2, 0, b2, -c2],
[0, c2, e2, 0, -c2, d2]])
invt2 = np.matrix.transpose(t2)
k2global = np.linalg.multi_dot([invt2, k2local, t2])
#Elmento 3 - asna 2
a3 = E * arasna / lasna
b3 = 12 * E *inyasna / lasna**3
c3 = 6 * E * inyasna / lasna**2
d3 = 4 * E * inyasna / lasna
e3 = 2 * E *inyasna / lasna
alfa3 = -alfa2
l3 = mt.cos(alfa3)
m3 = mt.sin(alfa3)
t3 = np.matrix([[l3, m3, 0, 0, 0, 0],
[-m3, l3, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l3, m3, 0],
[0, 0, 0, -m3, l3, 0],
[0 , 0, 0, 0, 0, 1]])
k3local = np.matrix([[a3, 0, 0, -a3, 0, 0],
[0, b3, c3, 0, -b3, c3],
[0, c3, d3, 0, -c3, e3],
[-a3, 0, 0, a3, 0, 0],
[0, -b3, -c3, 0, b3, -c3],
[0, c3, e3, 0, -c3, d3]])
invt3 = np.matrix.transpose(t3)
k3global = np.linalg.multi_dot([invt3, k3local, t3])
#Elmento 4 - pilar 2
a4 = E * arpilar / h
b4 = 12 * E *inypilar / h**3
c4 = 6 * E * inypilar / h**2
d4 = 4 * E * inypilar / h
e4 = 2 * E *inypilar / h
alfa4 = -(90 * mt.pi/180)
l4 = mt.cos(alfa4)
m4 = mt.sin(alfa4)
t4 = np.matrix([[l4, m4, 0, 0, 0, 0],
[-m4, l4, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l4, m4, 0],
[0, 0, 0, -m4, l4, 0],
[0 , 0, 0, 0, 0, 1]])
k4local = np.matrix([[a4, 0, 0, -a4, 0, 0],
[0, b4, c4, 0, -b4, c4],
[0, c4, d4, 0, -c4, e4],
[-a4, 0, 0, a4, 0, 0],
[0, -b4, -c4, 0, b4, -c4],
[0, c4, e4, 0, -c4, d4]])
invt4 = np.matrix.transpose(t4)
k4global = np.linalg.multi_dot([invt4, k4local, t4])
k = [k1global, k2global, k3global, k4global]
kportico = np.zeros([15,15])
for i,m in enumerate(k):
kportico[i*3:i*3+6,i*3:i*3+6] += m
#[K] * {U} = {F} - ELU
F12x = (comb_a1n * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F12y = (comb_a1n * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F22x = (comb_a1t * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F22y = (comb_a1t* lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F13x = F12x
F13y = F12y
F23x = F22x
F23y = F22y
F33x = (comb_a2n * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F33y = (comb_a2n * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F43x = (comb_a2t * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F43y = (comb_a2t * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F14x = F33x
F14y = F33y
F24x = F43x
F24y = F43y
F1x = (vnp1 * 1000 * h) / 2
F1y = 0
M1 = -(vnp1 * 1000 * h ** 2) / 12
F2x = int((vnp1 * 1000 * h / 2) + F12x - F22x)
F2y = - F12y - F22y
M2 = ((vnp1 * 1000 * h ** 2) / 12) - ((comb_a1n * lasna ** 2) / 12)
F3x = F13x - F23x - F33x + F43x
F3y = -F13y - F23y - F33y - F43y
M3 = ((comb_a1n * lasna ** 2) / 12) - ((comb_a2n * lasna ** 2) / 12)
F4x = (vnp2 * h * 1000 / 2) + F24x - F14x
F4y = - F14y - F24y
M4 = - ((vnp2 * 1000 * h ** 2) / 12) + ((comb_a2n * lasna ** 2) / 12)
F5x = (vnp2 * 1000 * h) / 2
F5y = 0
M5 = (vnp2 * 1000 * h ** 2) / 12
f = np.array([[F1x],
[F1y],
[M1],
[F2x],
[F2y],
[M2],
[F3x],
[F3y],
[M3],
[F4x],
[F4y],
[M4],
[F5x],
[F5y],
[M5]])
fel1 = np.array([[(vnp1 * 1000 * h) / 2],
[0],
[-(vnp1 * 1000 * h ** 2) / 12],
[(vnp1 * 1000 * h) / 2],
[0],
[(vnp1 * 1000 * h ** 2) / 12]])
fel2 = np.array([[(comb_a1n * lasna / 2) * mt.sin(alfa2) - (comb_a1t * lasna / 2) * mt.cos(alfa2)],
[- (comb_a1n * lasna / 2) * mt.cos(alfa2) - (comb_a1t * lasna / 2) * mt.sin(alfa2)],
[-(comb_a1n * lasna ** 2) / 12],
[(comb_a1n * lasna / 2) * mt.sin(alfa2) - (comb_a1t * lasna / 2) * mt.cos(alfa2)],
[- (comb_a1n * lasna / 2) * mt.cos(alfa2) - (comb_a1t * lasna / 2) * mt.sin(alfa2)],
[(comb_a1n * lasna ** 2) / 12]])
fel3 = np.array([[( - comb_a2n * lasna / 2) * mt.sin(-alfa3) + (comb_a2t * lasna / 2) * mt.cos(-alfa3)],
[( - comb_a2n * lasna / 2) * mt.cos(-alfa3) - (comb_a2t * lasna / 2) * mt.sin(-alfa3)],
[-(comb_a2n * lasna ** 2) / 12],
[( - comb_a2n * lasna / 2) * mt.sin(-alfa3) + (comb_a2t * lasna / 2) * mt.cos(-alfa3)],
[( - comb_a2n * lasna / 2) * mt.cos(-alfa3) - (comb_a2t * lasna / 2) * mt.sin(-alfa3)],
[(comb_a2n * lasna ** 2) / 12]])
fel4 = np.array([[(vnp2 * 1000 * h) / 2],
[0],
[-(vnp2 * 1000 * h ** 2) / 12],
[(vnp2 * 1000 * h) / 2],
[0],
[(vnp2 * 1000 * h ** 2) / 12]])
if base == rotulado:
u = np.dot(np.linalg.pinv(kportico[np.ix_([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14], [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14])]), f[np.ix_([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14])])
utotal = np.array([[0],
[0],
[u[0]],
[u[1]],
[u[2]],
[u[3]],
[u[4]],
[u[5]],
[u[6]],
[u[7]],
[u[8]],
[u[9]],
[0],
[0],
[u[10]]])
r = np.dot(kportico, utotal) - f
rp = np.array([[r[0]],
[r[1]],
[r[12]],
[r[13]]])
rv = np.array([["R1x ="],
["R1y ="],
["R5x ="],
["R5y ="]])
uni = np.array([["N"],
["N"],
["N"],
["N"]])
size = len(rv)
print(" ")
print(" ")
if len(rv) == len(rp) and len(rv) == len(uni):
for x in range(size):
print(rv[x],rp[x],uni[x])
Fpilar1 = np.dot(k1global, utotal[0:6]) - fel1 #ESFORÇOS NÓS PILAR 1
Fasna1 = np.dot(k2global, utotal[3:9]) - fel2 #ESFORÇOS NÓS ASNA 1
Fasna2 = np.dot(k3global, utotal[6:12]) - fel3 #ESFORÇOS NÓS ASNA 2
Fpilar2 = np.dot(k4global, utotal[9:15]) - fel4 #ESFORÇOS NÓS PILAR 2
#Diagrama de esforço transverso e momentos fletores
fig, ax = plt.subplots()
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
ha2 = np.arange(0, lasna, 0.1)
hp2 = np.arange(0, h+1, 1)
Mpilar1 = ((-vnp1 * 1000 * h1 ** 2 / 2) - float(Fpilar1[0]) * h1 - float(Fpilar1[2])) / 1000
Masna1 = (((- comb_a1n * ha1 ** 2 ) / 2) + (float(Fasna1[1]) * mt.cos(alfa2) - float(Fasna1[0]) * mt.sin(alfa2)) * ha1 - float(Fasna1[2])) / 1000
Masna2 = (((- comb_a2n * ha2 ** 2 ) / 2) + (float(Fasna2[1]) * mt.cos(alfa2) + float(Fasna2[0]) * mt.sin(alfa2)) * ha2 - float(Fasna2[2])) / 1000
Mpilar2 = ((vnp2 * 1000 * hp2 ** 2 / 2) + float(Fpilar2[0]) * hp2 - float(Fpilar2[2])) / 1000
# Roto-translation matrix:
# Rotates the points by an angle theta and translates
# them by x in the horizontal direction, and y in the
# vertical direction
R = lambda x, y, theta: np.array([
[np.cos(theta), np.sin(theta), x],
[-np.sin(theta), np.cos(theta), y],
[0, 0, 1],
])
# mirror matrix about the x-axis
Mx = np.array([
[1, 0, 0], [0, -1, 0], [0, 0, 1]
])
# mirror matrix about the y-axis
My = np.array([
[-1, 0, 0], [0, 1, 0], [0, 0, 1]
])
def draw_line():
plt.figure()
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
# left column
points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])
points1 = np.matmul(R(0, 0, -np.pi/2), points1)
plt.plot(points1[0, :], points1[1, :], label="Mpilar1")
# right column
points2 = np.stack([h1, Mpilar2 / max(Mpilar2), np.ones_like(h1)])
points2 = np.matmul(R(20, 0, -np.pi/2), points2)
plt.plot(points2[0, :], points2[1, :], label="Mpilar2")
# left asna
points3 = np.stack([ha1, Masna1 / max(Masna1), np.ones_like(ha1)])
points3 = np.matmul(R(0, 6, -alfa_rad), points3)
plt.plot(points3[0, :], points3[1, :], label="Masna1")
# right asna
points4 = np.stack([ha1, Masna2 / max(Masna2), np.ones_like(ha1)])
points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)
plt.plot(points4[0, :], points4[1, :], label="Masna2")
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.legend()
plt.show()
draw_line()
我是一名机械工程专业的学生，我正在尝试为 Python 中的力矩图创建一个脚本。我的代码中缺少的是如何定位力矩函数以便像门户框架一样对齐。
Mpilar1 是第一列（从左到右）的矩函数。
Masna1 是第一个梁的力矩函数（从左到右）。
Masna2 是第二个梁（从左到右）的力矩函数。
Mpilar2 是第二列（从左到右）的矩函数。
代码：
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
alfa = (mt.atan((int(ht)-int(h))/(int(v)/2)))*180/((mt.pi))
print("Ângulo da vertente:", round(alfa, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
def draw_line():
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.show()
if __name__ == '__main__':
draw_line()
Mpilar1 = 1500 * h1 ** 2 + 350 * h1
Masna1 = 300 * ha1 ** 2 + 15 * ha1
Masna2 = 200 * ha1 ** 2 + 15 * ha1
Mpilar2 = 1400 * h1 ** 2 + 10 * h1
plt.plot(h1, Mpilar1)
plt.plot(ha1, Masna1)
plt.plot(ha1, Masna2)
plt.plot(h1, Mpilar2)
您必须在代表曲线的点上使用变换矩阵。特别是，您必须使用 roto-translation 矩阵将曲线旋转并平移到正确的位置和方向，并且您可能必须应用镜像矩阵以使力矩根据您的约定对齐。
请注意，我在结构工程领域的日子已经一去不复返了，所以我不记得正确定位时刻的惯例。留给你作为练习。
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
alfa_rad = mt.atan((int(ht)-int(h))/(int(v)/2))
alfa_deg = alfa_rad*180/mt.pi
print("Ângulo da vertente:", round(alfa_deg, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
# Roto-translation matrix:
# Rotates the points by an angle theta and translates
# them by x in the horizontal direction, and y in the
# vertical direction
R = lambda x, y, theta: np.array([
[np.cos(theta), np.sin(theta), x],
[-np.sin(theta), np.cos(theta), y],
[0, 0, 1],
])
# mirror matrix about the x-axis
Mx = np.array([
[1, 0, 0], [0, -1, 0], [0, 0, 1]
])
# mirror matrix about the y-axis
My = np.array([
[-1, 0, 0], [0, 1, 0], [0, 0, 1]
])
Mpilar1 = 1500 * h1 ** 2 + 350 * h1
Masna1 = 300 * ha1 ** 2 + 15 * ha1
Masna2 = 200 * ha1 ** 2 + 15 * ha1
Mpilar2 = 1400 * h1 ** 2 + 10 * h1
def draw_line():
plt.figure()
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
# left column
points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])
points1 = np.matmul(R(0, 0, -np.pi/2), points1)
plt.plot(points1[0, :], points1[1, :], label="Mpilar1")
# right column
points2 = np.stack([h1, Mpilar2 / max(Mpilar2), np.ones_like(h1)])
points2 = np.matmul(R(20, 0, -np.pi/2), points2)
plt.plot(points2[0, :], points2[1, :], label="Mpilar2")
# left asna
points3 = np.stack([ha1, Masna1 / max(Masna1), np.ones_like(ha1)])
points3 = np.matmul(R(0, 6, -alfa_rad), points3)
plt.plot(points3[0, :], points3[1, :], label="Masna1")
# right asna
points4 = np.stack([ha1, Masna2 / max(Masna2), np.ones_like(ha1)])
points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)
plt.plot(points4[0, :], points4[1, :], label="Masna2")
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.legend()
plt.show()
draw_line()
上面的代码有几点需要注意：
让我们考虑一下
points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])。它创建一个 3xn 坐标矩阵。第一行是x-coordinates、h1。第二行是当前的 y-coordinates，Mpilar1 / max(Mpilar1)（请注意，为了适合图表，我已经对其进行了维度化）。第三行是1，能够应用平移矩阵是一个技巧。在绘图命令中，我们将只使用第一行和第二行（x 和 y 坐标）。points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)这里我先镜像了关于y-axis的点，然后我应用了旋转和平移。你必须玩才能正确定位时刻！
解决方案有效，但在这种情况下，我没有正确的力矩图。
代码：
import math as mt
import numpy as np
import warnings
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
#Definir parâmetros do pavilhão
v = 20 #(Vão em metros)
h = 6 #("Altura do pilar em metros:")
ht = 8 #("Altura total metros:")
nm = 7 #("Número de madres por asna:")
npo = 5 #("Número de pórticos:")
dp = 5 #("Distância entre pórticos em metros:")
alfa_rad = mt.atan((int(ht)-int(h))/(int(v)/2))
alfa_deg = alfa_rad*180/mt.pi
alfa = (mt.atan((int(ht)-int(h))/(int(v)/2)))*180/((mt.pi))
print("Ângulo da vertente:", round(alfa, 1), "º")
lasna = ((v/2) ** 2 + (ht-h) ** 2) ** 0.5
print("Comprimento de cada asna: ", round(lasna, 2), "m")
dm = lasna / nm
print(("Distância entre madres da cobertura em metros:"), round(dm, 2), "m")
lp = npo * dp
print("Comprimento total do pavilhão:", round(lp, 1), "m")
def draw_line():
# x axis value list.
x_number_list = [0, 0, (v/2), v, v]
# y axis value list.
y_number_list = [0, h, ht, h, 0]
# Plot the number in the list and set the line thickness.
plt.plot(x_number_list, y_number_list, linewidth=3)
# Set the line chart title and the text font size.
plt.title("Pórtico", fontsize=15)
# Set x axis label.
plt.xlabel("Vão (m)", fontsize=10)
# Set y axis label.
plt.ylabel("Altura (m)", fontsize=10)
# Set the x, y axis tick marks text size.
plt.tick_params(axis='both', labelsize=9)
# Display the plot in the matplotlib's viewer.
plt.show()
if __name__ == '__main__':
draw_line()
#EXEMPLO - IPE 100
pppilar = 8.1 #variavel do peso próprio do pilar
arpilar = 10.32 * 10 ** -4 #variavél da área de secção do pilar - m^2
inypilar = 171 * 10 ** -8 #variavél da inércia pilar segundo y - m^4
ppasna = 8.1 #variavel do peso próprio da asna
arasna = 10.32 * 10 ** -4 #variavél da área de secção do pilar - m^2
inyasna = 171 * 10 ** -8 #variavél da inércia pilar - m^4
pesomadre = 4.84 #(Peso linear das madres de cobertura em kg/m)
pesopainel = 11.2 #(Peso painel revestimento kg/m2)
rotulado = 1
base = rotulado
#Definir ações
ppcobertura = (pesopainel * dp + ((pesomadre*dp)/dm) + ppasna) * 9.81 / 1000 #kN/m
#print("Peso próprio da cobertura:", round(ppcobertura, 1), "kN/m")
#Sobrecarga
qk = 0.4 #Sobrecarga em kN/m2
sb = qk * dp #Sobrecarga em kN/m
#Neve
Sk = 0.5 #Neve em kN/m2
ne = Sk * dp #Neve em kN/m
#Vento
Qp1 = 0.5 #Vento pilar 1 em kN/m2
vnp1 = Qp1 * dp #Vento pilar 1 em kN/m
Qp2 = 0.5 #Vento pilar 2 em kN/m2
vnp2 = Qp2 * dp #Vento pilar 2 em kN/m
Qa1 = 0.5 #Vento asna 1 em kN/m2
vna1 = Qa1 * dp #Vento asna 1 em kN/m
Qa2 = 0.5 #Vento asna 2 em kN/m2
vna2 = Qa2 * dp #Vento asna 2 em kN/m
#Decompor as ações em normal e tagencial
ppcoberturan = ppcobertura * mt.cos(alfa*mt.pi/180)
ppcoberturat = ppcobertura * mt.sin(alfa*mt.pi/180)
sbn = sb * mt.cos(alfa*mt.pi/180) * mt.cos(alfa*mt.pi/180)
sbt = sb * mt.cos(alfa*mt.pi/180) * mt.sin(alfa*mt.pi/180)
nen = ne * mt.cos(alfa*mt.pi/180) * mt.cos(alfa*mt.pi/180)
net = ne * mt.cos(alfa*mt.pi/180) * mt.sin(alfa*mt.pi/180)
#Coeficientes de majoração
psipp = 1.35
#print("\u03A8 peso próprio:", psipp)
psi0sb = 0
#print("\u03A8 0 sobrecarga:", psi0sb)
psi1sb = 0
#print("\u03A8 1 sobrecarga:", psi1sb)
psi2sb = 0
#print("\u03A8 2 sobrecarga:", psi2sb)
ne1 = 1
ne2 = 2
nete = ne1
if nete == ne1:
psi0ne = 0.5
#print("\u03A8 0 neve:", psi0ne)
psi1ne = 0.2
#print("\u03A8 1 neve:", psi1ne)
psi2ne = 0
#print("\u03A8 2 neve:", psi2ne)
psi0vn = 0.6
#print("\u03A8 0 vento:", psi0vn)
psi1vn = 0.2
#print("\u03A8 1 vento:", psi1vn)
psi2vn = 0
#print("\u03A8 2 vento:", psi2vn)
#Combinação das ações para a cobertura - ELU - asna 1 - normal
comb_sbn = (psipp * ppcoberturan + sbn * 1.5 + (1.5 * psi0ne * nen + 1.5 * psi0vn * vna1)) * 1000 #N/m
comb_vnn = (psipp * ppcoberturan + vna1 * 1.5 + (1.5 * psi0ne * nen + 1.5 * sbn * psi0sb)) * 1000 #N/m
comb_nen =( psipp * ppcoberturan + nen * 1.5 + (1.5 * psi0vn * vna1 + 1.5 * sbn * psi0sb)) * 1000 #N/m
if (comb_sbn >= comb_vnn) and (comb_sbn >= comb_nen):
comb_a1n = comb_sbn
elif (comb_vnn >= comb_sbn) and (comb_vnn >= comb_nen):
comb_a1n = comb_vnn
else:
comb_a1n = comb_nen
#Combinação das ações para a cobertura - ELU - asna 1 - tangencial
comb_sbt = (psipp * ppcoberturat + sbt * 1.5 + (1.5 * psi0ne * net + 1.5 * psi0vn * 0)) * 1000 #N/m
comb_vnt = (psipp * ppcoberturat + 0 * 1.5 + (1.5 * psi0ne * net + 1.5 * sbt * psi0sb)) * 1000 #N/m
comb_net = (psipp * ppcoberturat + net * 1.5 + (1.5 * psi0vn * 0 + 1.5 * sbt * psi0sb)) * 1000 #N/m
if (comb_sbn >= comb_vnn) and (comb_sbn >= comb_nen):
comb_a1t = comb_sbt
elif (comb_vnn >= comb_sbn) and (comb_vnn >= comb_nen):
comb_a1t = comb_vnt
else:
comb_a1t = comb_net
#Combinação das ações para a cobertura - ELU - asna 2 - normal
comb_sb2n = (psipp * ppcoberturan + sbn * 1.5 + (1.5 * psi0ne * nen + 1.5 * psi0vn * vna2)) * 1000 #N/m
comb_vn2n = (psipp * ppcoberturan + vna2 * 1.5 + (1.5 * psi0ne * nen + 1.5 * sbn * psi0sb)) * 1000 #N/m
comb_ne2n = (psipp * ppcoberturan + nen * 1.5 + (1.5 * psi0vn * vna2 + 1.5 * sbn * psi0sb)) * 1000 #N/m
if (comb_sb2n >= comb_vn2n) and (comb_sb2n >= comb_ne2n):
comb_a2n = comb_sb2n
elif (comb_vn2n >= comb_sb2n) and (comb_vn2n >= comb_ne2n):
comb_a2n = comb_vn2n
else:
comb_a2n = comb_ne2n
#Combinação das ações para a cobertura - ELU - asna 2 - tangencial
comb_sbt2 = (psipp * ppcoberturat + sbt * 1.5 + (1.5 * psi0ne * net + 1.5 * psi0vn * 0)) * 1000 #N/m
comb_vnt2 = (psipp * ppcoberturat + 0 * 1.5 + (1.5 * psi0ne * net + 1.5 * sbt * psi0sb)) * 1000 #N/m
comb_net2 = (psipp * ppcoberturat + net * 1.5 + (1.5 * psi0vn * 0 + 1.5 * sbt * psi0sb)) * 1000 #N/m
if (comb_sb2n >= comb_vn2n) and (comb_sb2n >= comb_ne2n):
comb_a2t = comb_sbt2
elif (comb_vn2n >= comb_sb2n) and (comb_vn2n >= comb_ne2n):
comb_a2t = comb_vnt2
else:
comb_a2t = comb_net2
#Elementos finitos - Reações e deslocamentos
E = 210 * 10 ** 9 #módulo de elasticidade do aço em Pa
#Elemento 1 - asna1
a1 = E * arpilar / h
b1 = 12 * E * inypilar / h**3
c1 = 6 * E * inypilar / h**2
d1 = 4 * E * inypilar / h
e1 = 2 * E * inypilar / h
alfa1 = 90 * mt.pi / 180
l1 = mt.cos(alfa1)
m1 = mt.sin(alfa1)
t1 = np.matrix([[l1, m1, 0, 0, 0, 0],
[-m1, l1, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l1, m1, 0],
[0, 0, 0, -m1, l1, 0],
[0 , 0, 0, 0, 0, 1]])
k1local = np.matrix([[a1, 0, 0, -a1, 0, 0],
[0, b1, c1, 0, -b1, c1],
[0, c1, d1, 0, -c1, e1],
[-a1, 0, 0, a1, 0, 0],
[0, -b1, -c1, 0, b1, -c1],
[0, c1, e1, 0, -c1, d1]])
invt1 = np.matrix.transpose(t1)
k1global = np.linalg.multi_dot([invt1, k1local, t1])
#Elmento 2 - asna 1
a2 = E * arasna / lasna
b2 = 12 * E *inyasna / lasna**3
c2 = 6 * E * inyasna / lasna**2
d2 = 4 * E * inyasna / lasna
e2 = 2 * E *inyasna / lasna
alfa2 = ((alfa) * mt.pi) / 180
l2 = mt.cos(alfa2)
m2 = mt.sin(alfa2)
t2 = np.matrix([[l2, m2, 0, 0, 0, 0],
[-m2, l2, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l2, m2, 0],
[0, 0, 0, -m2, l2, 0],
[0 , 0, 0, 0, 0, 1]])
k2local = np.matrix([[a2, 0, 0, -a2, 0, 0],
[0, b2, c2, 0, -b2, c2],
[0, c2, d2, 0, -c2, e2],
[-a2, 0, 0, a2, 0, 0],
[0, -b2, -c2, 0, b2, -c2],
[0, c2, e2, 0, -c2, d2]])
invt2 = np.matrix.transpose(t2)
k2global = np.linalg.multi_dot([invt2, k2local, t2])
#Elmento 3 - asna 2
a3 = E * arasna / lasna
b3 = 12 * E *inyasna / lasna**3
c3 = 6 * E * inyasna / lasna**2
d3 = 4 * E * inyasna / lasna
e3 = 2 * E *inyasna / lasna
alfa3 = -alfa2
l3 = mt.cos(alfa3)
m3 = mt.sin(alfa3)
t3 = np.matrix([[l3, m3, 0, 0, 0, 0],
[-m3, l3, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l3, m3, 0],
[0, 0, 0, -m3, l3, 0],
[0 , 0, 0, 0, 0, 1]])
k3local = np.matrix([[a3, 0, 0, -a3, 0, 0],
[0, b3, c3, 0, -b3, c3],
[0, c3, d3, 0, -c3, e3],
[-a3, 0, 0, a3, 0, 0],
[0, -b3, -c3, 0, b3, -c3],
[0, c3, e3, 0, -c3, d3]])
invt3 = np.matrix.transpose(t3)
k3global = np.linalg.multi_dot([invt3, k3local, t3])
#Elmento 4 - pilar 2
a4 = E * arpilar / h
b4 = 12 * E *inypilar / h**3
c4 = 6 * E * inypilar / h**2
d4 = 4 * E * inypilar / h
e4 = 2 * E *inypilar / h
alfa4 = -(90 * mt.pi/180)
l4 = mt.cos(alfa4)
m4 = mt.sin(alfa4)
t4 = np.matrix([[l4, m4, 0, 0, 0, 0],
[-m4, l4, 0, 0, 0, 0],
[0, 0 , 1, 0, 0, 0],
[0, 0, 0, l4, m4, 0],
[0, 0, 0, -m4, l4, 0],
[0 , 0, 0, 0, 0, 1]])
k4local = np.matrix([[a4, 0, 0, -a4, 0, 0],
[0, b4, c4, 0, -b4, c4],
[0, c4, d4, 0, -c4, e4],
[-a4, 0, 0, a4, 0, 0],
[0, -b4, -c4, 0, b4, -c4],
[0, c4, e4, 0, -c4, d4]])
invt4 = np.matrix.transpose(t4)
k4global = np.linalg.multi_dot([invt4, k4local, t4])
k = [k1global, k2global, k3global, k4global]
kportico = np.zeros([15,15])
for i,m in enumerate(k):
kportico[i*3:i*3+6,i*3:i*3+6] += m
#[K] * {U} = {F} - ELU
F12x = (comb_a1n * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F12y = (comb_a1n * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F22x = (comb_a1t * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F22y = (comb_a1t* lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F13x = F12x
F13y = F12y
F23x = F22x
F23y = F22y
F33x = (comb_a2n * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F33y = (comb_a2n * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F43x = (comb_a2t * lasna / 2) * mt.cos((alfa*mt.pi) / 180)
F43y = (comb_a2t * lasna / 2) * mt.sin((alfa*mt.pi) / 180)
F14x = F33x
F14y = F33y
F24x = F43x
F24y = F43y
F1x = (vnp1 * 1000 * h) / 2
F1y = 0
M1 = -(vnp1 * 1000 * h ** 2) / 12
F2x = int((vnp1 * 1000 * h / 2) + F12x - F22x)
F2y = - F12y - F22y
M2 = ((vnp1 * 1000 * h ** 2) / 12) - ((comb_a1n * lasna ** 2) / 12)
F3x = F13x - F23x - F33x + F43x
F3y = -F13y - F23y - F33y - F43y
M3 = ((comb_a1n * lasna ** 2) / 12) - ((comb_a2n * lasna ** 2) / 12)
F4x = (vnp2 * h * 1000 / 2) + F24x - F14x
F4y = - F14y - F24y
M4 = - ((vnp2 * 1000 * h ** 2) / 12) + ((comb_a2n * lasna ** 2) / 12)
F5x = (vnp2 * 1000 * h) / 2
F5y = 0
M5 = (vnp2 * 1000 * h ** 2) / 12
f = np.array([[F1x],
[F1y],
[M1],
[F2x],
[F2y],
[M2],
[F3x],
[F3y],
[M3],
[F4x],
[F4y],
[M4],
[F5x],
[F5y],
[M5]])
fel1 = np.array([[(vnp1 * 1000 * h) / 2],
[0],
[-(vnp1 * 1000 * h ** 2) / 12],
[(vnp1 * 1000 * h) / 2],
[0],
[(vnp1 * 1000 * h ** 2) / 12]])
fel2 = np.array([[(comb_a1n * lasna / 2) * mt.sin(alfa2) - (comb_a1t * lasna / 2) * mt.cos(alfa2)],
[- (comb_a1n * lasna / 2) * mt.cos(alfa2) - (comb_a1t * lasna / 2) * mt.sin(alfa2)],
[-(comb_a1n * lasna ** 2) / 12],
[(comb_a1n * lasna / 2) * mt.sin(alfa2) - (comb_a1t * lasna / 2) * mt.cos(alfa2)],
[- (comb_a1n * lasna / 2) * mt.cos(alfa2) - (comb_a1t * lasna / 2) * mt.sin(alfa2)],
[(comb_a1n * lasna ** 2) / 12]])
fel3 = np.array([[( - comb_a2n * lasna / 2) * mt.sin(-alfa3) + (comb_a2t * lasna / 2) * mt.cos(-alfa3)],
[( - comb_a2n * lasna / 2) * mt.cos(-alfa3) - (comb_a2t * lasna / 2) * mt.sin(-alfa3)],
[-(comb_a2n * lasna ** 2) / 12],
[( - comb_a2n * lasna / 2) * mt.sin(-alfa3) + (comb_a2t * lasna / 2) * mt.cos(-alfa3)],
[( - comb_a2n * lasna / 2) * mt.cos(-alfa3) - (comb_a2t * lasna / 2) * mt.sin(-alfa3)],
[(comb_a2n * lasna ** 2) / 12]])
fel4 = np.array([[(vnp2 * 1000 * h) / 2],
[0],
[-(vnp2 * 1000 * h ** 2) / 12],
[(vnp2 * 1000 * h) / 2],
[0],
[(vnp2 * 1000 * h ** 2) / 12]])
if base == rotulado:
u = np.dot(np.linalg.pinv(kportico[np.ix_([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14], [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14])]), f[np.ix_([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14])])
utotal = np.array([[0],
[0],
[u[0]],
[u[1]],
[u[2]],
[u[3]],
[u[4]],
[u[5]],
[u[6]],
[u[7]],
[u[8]],
[u[9]],
[0],
[0],
[u[10]]])
r = np.dot(kportico, utotal) - f
rp = np.array([[r[0]],
[r[1]],
[r[12]],
[r[13]]])
rv = np.array([["R1x ="],
["R1y ="],
["R5x ="],
["R5y ="]])
uni = np.array([["N"],
["N"],
["N"],
["N"]])
size = len(rv)
print(" ")
print(" ")
if len(rv) == len(rp) and len(rv) == len(uni):
for x in range(size):
print(rv[x],rp[x],uni[x])
Fpilar1 = np.dot(k1global, utotal[0:6]) - fel1 #ESFORÇOS NÓS PILAR 1
Fasna1 = np.dot(k2global, utotal[3:9]) - fel2 #ESFORÇOS NÓS ASNA 1
Fasna2 = np.dot(k3global, utotal[6:12]) - fel3 #ESFORÇOS NÓS ASNA 2
Fpilar2 = np.dot(k4global, utotal[9:15]) - fel4 #ESFORÇOS NÓS PILAR 2
#Diagrama de esforço transverso e momentos fletores
fig, ax = plt.subplots()
h1 = np.arange(0, h+1, 1)
ha1 = np.arange(0, lasna, 0.1)
ha2 = np.arange(0, lasna, 0.1)
hp2 = np.arange(0, h+1, 1)
Mpilar1 = ((-vnp1 * 1000 * h1 ** 2 / 2) - float(Fpilar1[0]) * h1 - float(Fpilar1[2])) / 1000
Masna1 = (((- comb_a1n * ha1 ** 2 ) / 2) + (float(Fasna1[1]) * mt.cos(alfa2) - float(Fasna1[0]) * mt.sin(alfa2)) * ha1 - float(Fasna1[2])) / 1000
Masna2 = (((- comb_a2n * ha2 ** 2 ) / 2) + (float(Fasna2[1]) * mt.cos(alfa2) + float(Fasna2[0]) * mt.sin(alfa2)) * ha2 - float(Fasna2[2])) / 1000
Mpilar2 = ((vnp2 * 1000 * hp2 ** 2 / 2) + float(Fpilar2[0]) * hp2 - float(Fpilar2[2])) / 1000
# Roto-translation matrix:
# Rotates the points by an angle theta and translates
# them by x in the horizontal direction, and y in the
# vertical direction
R = lambda x, y, theta: np.array([
[np.cos(theta), np.sin(theta), x],
[-np.sin(theta), np.cos(theta), y],
[0, 0, 1],
])
# mirror matrix about the x-axis
Mx = np.array([
[1, 0, 0], [0, -1, 0], [0, 0, 1]
])
# mirror matrix about the y-axis
My = np.array([
[-1, 0, 0], [0, 1, 0], [0, 0, 1]
])
def draw_line():
plt.figure()
x_number_list = [0, 0, (v/2), v, v]
y_number_list = [0, h, ht, h, 0]
plt.plot(x_number_list, y_number_list, linewidth=3)
# left column
points1 = np.stack([h1, Mpilar1 / max(Mpilar1), np.ones_like(h1)])
points1 = np.matmul(R(0, 0, -np.pi/2), points1)
plt.plot(points1[0, :], points1[1, :], label="Mpilar1")
# right column
points2 = np.stack([h1, Mpilar2 / max(Mpilar2), np.ones_like(h1)])
points2 = np.matmul(R(20, 0, -np.pi/2), points2)
plt.plot(points2[0, :], points2[1, :], label="Mpilar2")
# left asna
points3 = np.stack([ha1, Masna1 / max(Masna1), np.ones_like(ha1)])
points3 = np.matmul(R(0, 6, -alfa_rad), points3)
plt.plot(points3[0, :], points3[1, :], label="Masna1")
# right asna
points4 = np.stack([ha1, Masna2 / max(Masna2), np.ones_like(ha1)])
points4 = np.matmul(np.matmul(R(20, 6, alfa_rad), My), points4)
plt.plot(points4[0, :], points4[1, :], label="Masna2")
plt.title("Pórtico", fontsize=15)
plt.xlabel("Vão (m)", fontsize=10)
plt.ylabel("Altura (m)", fontsize=10)
plt.tick_params(axis='both', labelsize=9)
plt.legend()
plt.show()
draw_line()