Python: 単位ベクトルを使う（断面力図への応用） (2023.12.29)

はじめに

下に示すような図を作った。

この図は２次元骨組構造解析で得られた断面力を図にしたものである。基本的な書き方は、要素軸線に対し直行する方向に解析で得られた断面力の値をプロットし繋いでいくものである。ここでのキーは、「要素軸に直行する方向に値をプロットする」ことである。
以前は角度を計算してゴチョゴチョとやっていたのだが、何かもっと美しい方法はないかと考えていたところ、単位ベクトルというものがあったことを思い出し、活用してみることにした。

理論

まずは、ベクトルの回転から。
あるベクトル ${\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()

以　上