はじめに
下に示すような図を作った。
この図は２次元骨組構造解析で得られた断面力を図にしたものである。基本的な書き方は、要素軸線に対し直行する方向に解析で得られた断面力の値をプロットし繋いでいくものである。ここでのキーは、「要素軸に直行する方向に値をプロットする」ことである。
以前は角度を計算してゴチョゴチョとやっていたのだが、何かもっと美しい方法はないかと考えていたところ、単位ベクトルというものがあったことを思い出し、活用してみることにした。
理論
まずは、ベクトルの回転から。
あるベクトル ${\boldsymbol{a}}=(x_a,y_a)$ を角度 $\theta$ で回転させた後のベクトル ${\boldsymbol{b}}=(x_b,y_b)$ は以下のようにして得られる。
\begin{equation}
\begin{Bmatrix}
x_b \\
y_b
\end{Bmatrix}
=
\begin{bmatrix}
\cos\theta & -\sin\theta \\
\sin\theta & \cos\theta \\
\end{bmatrix}
\begin{Bmatrix}
x_a \\
y_a
\end{Bmatrix}
\end{equation}
ここで、回転角を $\theta=90^o$ とすれば
\begin{equation}
\begin{Bmatrix}
x_b \\
y_b
\end{Bmatrix}
=
\begin{bmatrix}
0 & -1 \\
1 & 0 \\
\end{bmatrix}
\begin{Bmatrix}
x_a \\
y_a
\end{Bmatrix}
=
\begin{Bmatrix}
-y_a \\
x_a
\end{Bmatrix}
\end{equation}
と簡単な形となる。
この関係を用いて、要素軸に直行する方向に与えられた距離だけ離れた点をプロットすることを考える。
計算のステップは以下の通り。
- ２つの点 Node-1、Node-2 およびこれらを結ぶ線分を定義する。
Node-1の座標を、$(x_1,y_1)$,、Node-2 の座標を、$(x_2,y_2)$ とする。
＋線分1-2 上にこれに沿う方向の単位ベクトル ${\boldsymbol{e_1}}$ を計算する。 - 単位ベクトル ${\boldsymbol{e_1}}$ に直行する単位ベクトル ${\boldsymbol{e_2}}$ を求める。
- 求めた単位ベクトル ${\boldsymbol{e_2}}$ を用い、Node-1 および Node-2 から与えられた距離 $D_1$お よび $D_2$ を持つ点 $(x_4,y_4)$ および $(x_3,y_3)$ を求める。
上記を式で表せば、以下の通りとなる。（上図参照）
\begin{equation}
\ell=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}
\end{equation}
\begin{align}
&\{\boldsymbol{e_1}\}=(x_s,y_s)=\left(\cfrac{x_2-x_1}{\ell},\cfrac{y_2-y_1}{\ell}\right) \\
&\{\boldsymbol{e_1}\}=(x_t,y_t)=(-y_s,x_s)
\end{align}
\begin{align}
&(x_4,y_4)=(x_1+D_1\cdot x_t, \quad y_1+D_1\cdot y_t) \\
&(x_3,y_3)=(x_2+D_2\cdot x_t, \quad y_2+D_2\cdot y_t)
\end{align}
これらを活用することにより、冒頭で示した断面力ずを描画することができるわけである。
プログラム例
一応、この投稿で用いた説明図を描くコードを載せておく。
import numpy as np
import matplotlib.pyplot as plt
def sarrow(x1,y1,x2,y2):
col='#0000ff'
arst='->,head_width=3,head_length=10'
aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=3)
plt.annotate('',
xy=(x1,y1), xycoords='data',
xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)
def dtext(xs,ys,ss,col,sha,sva,fszl,ang):
plt.text(xs,ys,ss,color=col,ha=sha,va=sva,fontsize=fszl,rotation=ang)
def main():
# node coordinates and section force
d1,d2=1.5,3 # section force
xx1,yy1=2.5,1.0 # node-1
xx2,yy2=6.5,4.0 # node-2
# calculation of unit vector
al=np.sqrt((xx2-xx1)**2+(yy2-yy1)**2)
xxs,yys=(xx2-xx1)/al,(yy2-yy1)/al
xxt,yyt=-yys,xxs
# calculation of plotting point
xx3,yy3=xx2+d2*xxt,yy2+d2*yyt
xx4,yy4=xx1+d1*xxt,yy1+d1*yyt
fsz=14
xmin,xmax,dx=0,8,1
ymin,ymax,dy=0,8,1
iw=6
ih=iw/(xmax-xmin)*(ymax-ymin)
plt.figure(figsize=(iw*2,ih),facecolor='w')
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
for iii in [1,2]:
if iii==1: tstr='Unit vector'
if iii==2: tstr='Plot of section force'
nplt=120+iii
plt.subplot(nplt)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xlabel('x_axis')
plt.ylabel('y_axis')
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.gca().set_aspect('equal',adjustable='box')
plt.title(tstr,loc='center',fontsize=fsz)
# element and node
plt.plot([xx1,xx2],[yy1,yy2],'o-',ms=8,color='#aaaaaa',lw=3)
dsx,dsy=0.7,0.5
col,sha,sva,fszl,ang='#000000','center','center',fsz,0
for i in [1,2]:
if i==1: xs,ys,ss=xx1,yy1-dsy,'1\n(x1,y1)'
if i==2: xs,ys,ss=xx2+dsx,yy2,'2\n(x2,y2)'
dtext(xs,ys,ss,col,sha,sva,fszl,ang)
if iii==1:
# plot of unit vector
dv=0.2
col,sha,sva,fszl,ang='#0000ff','center','center',fsz,0
for i in [1,2]:
if i==1:
x1,y1,x2,y2=xx1+xxs,yy1+yys,xx1,yy1
xs,ys,ss=0.5*(x1+x2)+dv,0.5*(y1+y2)-dv,'e1=(xs,ys)'
sha,sva='left','center'
if i==2:
x1,y1,x2,y2=xx1+xxt,yy1+yyt,xx1,yy1
xs,ys,ss=0.5*(x1+x2)-dv,0.5*(y1+y2-dv),'e2=(xt,yt)'
sha,sva='right','center'
sarrow(x1,y1,x2,y2)
dtext(xs,ys,ss,col,sha,sva,fszl,ang)
if iii==2:
# plot of diagram
xx=np.array([xx1,xx2,xx3,xx4])
yy=np.array([yy1,yy2,yy3,yy4])
plt.fill(xx,yy,color='#bbbbbb', alpha=0.5)
col,sha,sva,fszl,ang='#000000','center','center',fsz,0
for i in [1,2]:
if i==1: xs,ys,ss=xx3,yy3+dsy,'3\n(x3,y3)'
if i==2: xs,ys,ss=xx4-dsx,yy4,'4\n(x4,y4)'
dtext(xs,ys,ss,col,sha,sva,fszl,ang)
ang=np.degrees(np.arctan2(yyt,xxt))-180
ds=0.2
dsx=ds*np.cos(np.radians(ang+90))
dsy=ds*np.sin(np.radians(ang+90))
col,sha,sva,fszl,ang='#ff0000','center','center',fsz+2,ang
for i in [1,2]:
if i==1: xs,ys,ss=0.5*(xx1+xx4)+dsx,0.5*(yy1+yy4)+dsy,'D1'
if i==2: xs,ys,ss=0.5*(xx2+xx3)+dsx,0.5*(yy2+yy3)+dsy,'D2'
dtext(xs,ys,ss,col,sha,sva,fszl,ang)
plt.plot([xx1,xx4],[yy1,yy4],'-',color=col,lw=1)
plt.plot([xx2,xx3],[yy2,yy3],'-',color=col,lw=1)
plt.tight_layout()
fnameF='fig_ex.jpg'
plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
#plt.show()
#---------------
# Execute
#---------------
if __name__ == '__main__': main()
以 上