はじめに

格子桁構造解析プログラムを紹介します。このプログラムは、平面骨組構造解析プログラムの軸力・軸方向変位をねじりモーメント・ねじり角に書き換えたものです。構造は平面形状とし、x - y 平面上で定義します。

格子桁解析は、立体骨組構造解析プログラムで行うことができるのですが、床あるいは壁のような単純な板状構造物の解析には、入力データの単純さを考えると格子桁解析も有用であることから、作ったものです。

なお、このプログラムでは温度変化による断面力は考慮できません。それは、考慮できる断面力は、ねじりモーメント、曲げモーメント、せん断力であり、軸力項が無いためです。

（注意事項）

充実矩形断面のねじり定数 $J$ は、断面の短辺を $a$ 、長辺を $b$ として、以下の式で計算できます。

$\begin{equation} J=\cfrac{1}{3} b a^3 \left(1-\cfrac{2}{\pi} \cfrac{a}{b}\right) \end{equation}$

ここで、コンクリート床のような版構造を格子桁でモデル化する場合、ねじり定数 $J$ はゼロとみなせるくらいに小さな値を入力します。これは、ある要素のねじり角は、それに交差する要素の曲げ角に影響するため、曲げモーメントがねじりの影響で不連続になってしまうためです。

本来の格子桁をモデル化する場合は、もちろん、計算により求めた正規のねじり定数を入力する必要があります。

Grid Girder Analysis

Outline of Grid Girder Analysis Program

This is a program for Grid Girder Analysis. Elastic and small displacement problems can be treated. This program has been made using 'Python 3.'
Beam element with 2 nodes is used. 1 node has 3 degrees of freedom which are torsional rotation, bending rotation and deflection.
Simultaneous linear equations are solved using numpy.linalg.solve.
Input / Output file name can be defined arbitrarily as text files with the separator of blank, and those are inputted from command line of a terminal.

Workable condition
Item	Description
Element	Beam element which has 3 degrees of freedom such as torsional rotation, bending rotation and deflection.
Load	Specify the node number and load values on the node.
Displacement	Specify the node number and displacements on the node. Any values can be applied including zero

Global Coordinate System	Local Coordinate Syatem

The Stiffness matrix of an element is shown below:

$\begin{equation} \begin{Bmatrix} T_i \\ M_i \\ Q_i \\ T_j \\ M_j \\ Q_j \end{Bmatrix} =\begin{bmatrix} GJ/L & 0 & 0 & -GJ/L & 0 & 0 \\ 0 & 4 EI/L & -6 EI/L^2 & 0 & 2 EI/L & 6 EI/L^2 \\ 0 & -6 EI/L^2 & 12 EI/L^3 & 0 & -6 EI/L^2 & -12 EI/L^3 \\ -GJ/L & 0 & 0 & GJ/L & 0 & 0 \\ 0 & 2 EI/L & -6 EI/L^2 & 0 & 4 EI/L & 6 EI/L^2 \\ 0 & 6 EI/L^2 & -12 EI/L^3 & 0 & 6 EI/L^2 & 12 EI/L^3 \end{bmatrix} \begin{Bmatrix} \phi_i \\ \theta_i \\ w_i \\ \phi_j \\ \theta_j \\ w_j \end{Bmatrix} \end{equation}$

$\begin{align} &GJ & &\text{: Torsional rigidity} & &T & &\text{: Torsional moment} & &\phi & &\text{: Tortional rotation around x-axis}\\ &EI & &\text{: Bending rigidity} & &M & &\text{: Bending moment} & &\theta & &\text{: Bending rotation around y-axis}\\ &L & &\text{: Length of element} & &Q && \text{: Shearing force} & &w & &\text{: Displacement in z-direction} \end{align}$

A shear modulus $G$ is calculated from an elastic modulus of $E$ and Poisson's ration $\nu$ using following equation in the program.

$\begin{equation} G=\cfrac{E}{2 (1+\nu)} \end{equation}$

Since the coordinate transformation is carried out on the x-y plane only, the transfoemation matrix is the same type as it for 2D frame analysis.

Format of input data file

npoin nele nsec npfix nlod  # Basic values for analysis
E po AA AI AJ gamma         # material properties
    ....(1 to nsec)....     #
node-1,node-2,matno         # Element connectivity, material set number, uniformly distributed load per unit length
    ....(1 to nele)....     #
x,y                         # Node coordinates, temperature change of node
    ....(1 to npoin)....    #
node nokx noky nokz rdisx rdisy rdisz # node number, restrict conditions and forced displacement
    ....(1 to npfix)....    #
node,fx,fy,fz               # node number, Load values of torsional moment, bending moment, vertical force
    ....(1 to nlod)....     # (Omit data input if nlod=0)

npoin : Number of nodes                   E     : Elastic modulus of element
nele  : Number of elements                po    : Poison's ratio
nsec  : Number of material sets           AA    : Section area of element
npfix : Number of restricted nodes        AI    : Moment of Inertia of element
nlod  : Number of loaded nodes            AJ    : Tosion constant
                                          gamma : unit weight of material
                                          matno : Material set number

The structure shall be defined on the x-y plane in the global coordinate system.
x-direction in the local coordinate system means the member axis.
The displacement in local x-direction means the rotation around the x-axis, which is same as torsional rotation.
The displacement in local y-direction means the rotation around the y-axis, which is same as bending rotation.
The displacement in local z-direction means the deflection of the beam.
Restricted node means the node which has known (given) displacement. As a known (given) value of nodal displacement, any value can be given including zero for a restricted node.
Since stress resultants of element are defined as equivalent nodal forces in local coordinate system, it is necessary to note that signs are different from it on general structural mechanics. Positive directions are always right-direction, upward-direction for each node in local coordinate system.

Format of output file

npoin  nele  nsec npfix  nlod
    (Each value for above item)
  sec   E   po   A   I   J   gamma
    (material caracteristics)
    .....(1 to nsec).....
 node   x   y   fx   fy   fz   kox   koy   koz
    node: node number
    x   : x-coordinate
    y   : y-coordinate
    fx  : torsional load
    fy  : bending moment load
    fz  : vertical force load
    kox : restriction around x-axis
    koy : restriction around y-axis
    koz : restriction in z-direction
    .....(1 to npoin).....
 node   kox   koy   koz   rdis_x   rdis_y   rdis_z
    node   : node number
    kox    : restriction around x-axis
    koy    : restriction around y-axis
    koz    : restriction in z-direction
    rdis_x : given rotation around x-axis
    rdis_y : given rotation around y-axis
    rdis_z : given displacement in z-direction
    .....(1 to npfix).....
 elem   i   j   sec
    elem : element number
    i    : node-1
    j    : node-2
    sec  : material number
    .....(1 to nele).....
 node   rot-x   rot-y   dis-z
    node  : node number
    rot-x : rotation arouns x-axis
    rot-y : rotation around y-axis
    dis-z : displacement in z-direction
    .....(1 to npoin).....
 elem   T_i   M_i   Q_i   T_j   M_j   Q_j
    elem : element number
    T_i  : torsional moment at node-i
    M_i  : bending moment at node-i
    Q_i  : shearing force at node-i
    T_j  : torsional moment at node-j
    M_j  : bending moment at node-j
    Q_j  : shearing force at node-j
    .....(1 to nele).....
n=693  time=0.078 sec   # reference information

Grid Girder Analysis Program

###################################
# Grid Girder Analysis
###################################
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_grd(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    # array declaration
    x     =np.zeros([2,npoin],dtype=np.float64)         # Coordinates of nodes
    ae    =np.zeros([6,nsec],dtype=np.float64)          # Section characteristics
    node  =np.zeros([nod+1,nele],dtype=np.int)          # Node-element relationship
    fp    =np.zeros(nfree*npoin,dtype=np.float64)       # External force vector
    mpfix =np.zeros([nfree,npoin],dtype=np.int)         # Ristrict conditions
    rdis  =np.zeros([nfree,npoin],dtype=np.float64)     # Ristricted displacement
    # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) # E  : Elastic modulus
        ae[1,i]=float(text[1]) # po : Poisson's ratio
        ae[2,i]=float(text[2]) # AA : Section area
        ae[3,i]=float(text[3]) # AI : Moment of inertia
        ae[4,i]=float(text[4]) # AJ : Tortional constant
        ae[5,i]=float(text[5]) # gamma : Unit weight of beam material
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]=float(text[0])    # x-coordinate
        x[1,i]=float(text[1])    # y-coordinate
    # boundary conditions (0:free, 1:restricted)
    for i in range(0,npfix):
        text=f.readline()
        text=text.strip()
        text=text.split()
        lp=int(text[0])             #fixed node
        mpfix[0,lp-1]=int(text[1])  #index of restriction for tortional rotation
        mpfix[1,lp-1]=int(text[2])  #index of restriction for bending rotation
        mpfix[2,lp-1]=int(text[3])  #index of restriction for vertical displacement
        rdis[0,lp-1]=float(text[4]) #given tortional rotation
        rdis[1,lp-1]=float(text[5]) #given bending rotation
        rdis[2,lp-1]=float(text[6]) #given vertical displacement
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[3*lp-3]=float(text[1]) #tortional moment load
            fp[3*lp-2]=float(text[2]) #bending moment load
            fp[3*lp-1]=float(text[3]) #vertical load
    f.close()
    return npoin, nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp


def prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('sec','E','po','A','I','J','gamma'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s} {8:>5s}'
    .format('node','x','y','fx','fy','fz','kox','koy','koz'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d} {8:5d}'
        .format(lp,x[0,i],x[1,i],fp[3*lp-3],fp[3*lp-2],fp[3*lp-1],mpfix[0,i],mpfix[1,i],mpfix[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>15s} {5:>15s} {6:>15s}'
    .format('node','kox','koy','koz','rdis_x','rdis_y','rdis_z'),file=fout)
    for i in range(0,npoin):
        if 0<mpfix[0,i]+mpfix[1,i]+mpfix[2,i]:
            lp=i+1
            print('{0:5d} {1:5d} {2:5d} {3:5d} {4:15.7e} {5:15.7e} {6:15.7e}'
            .format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],rdis[0,i],rdis[1,i],rdis[2,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s}'.format('elem','i','j','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
    fout.close()

    
def prout_grd(fnameW,npoin,nele,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'.format('node','rot-x','rot-y','dis-z'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'.format(lp,disg[3*lp-3],disg[3*lp-2],disg[3*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s}'
    .format('elem','T_i','M_i','Q_i','T_j','M_j','Q_j'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
        .format(ne+1,fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
    fout.close()


def sm_grd(GJ,EI,x1,y1,x2,y2):
    ek=np.zeros([6,6],dtype=np.float64) # local stiffness matrix
    tt=np.zeros([6,6],dtype=np.float64) # transformation matrix
    xx=x2-x1
    yy=y2-y1
    el=np.sqrt(xx**2+yy**2)
    ek[0,0]= GJ/el; ek[0,1]=        0.0; ek[0,2]=         0.0; ek[0,3]=-GJ/el; ek[0,4]=        0.0; ek[0,5]=         0.0
    ek[1,0]=   0.0; ek[1,1]= 4*EI/el   ; ek[1,2]= -6*EI/el**2; ek[1,3]=   0.0; ek[1,4]= 2*EI/el   ; ek[1,5]=  6*EI/el**2
    ek[2,0]=   0.0; ek[2,1]=-6*EI/el**2; ek[2,2]= 12*EI/el**3; ek[2,3]=   0.0; ek[2,4]=-6*EI/el**2; ek[2,5]=-12*EI/el**3
    ek[3,0]=-GJ/el; ek[3,1]=        0.0; ek[3,2]=         0.0; ek[3,3]= GJ/el; ek[3,4]=        0.0; ek[3,5]=         0.0
    ek[4,0]=   0.0; ek[4,1]= 2*EI/el   ; ek[4,2]= -6*EI/el**2; ek[4,3]=   0.0; ek[4,4]= 4*EI/el   ; ek[4,5]=  6*EI/el**2
    ek[5,0]=   0.0; ek[5,1]= 6*EI/el**2; ek[5,2]=-12*EI/el**3; ek[5,3]=   0.0; ek[5,4]= 6*EI/el**2; ek[5,5]= 12*EI/el**3
    s=yy/el
    c=xx/el
    tt[0,0]=  c; tt[0,1]=  s; tt[0,2]=0.0; tt[0,3]=0.0; tt[0,4]=0.0; tt[0,5]=0.0
    tt[1,0]= -s; tt[1,1]=  c; tt[1,2]=0.0; tt[1,3]=0.0; tt[1,4]=0.0; tt[1,5]=0.0
    tt[2,0]=0.0; tt[2,1]=0.0; tt[2,2]=1.0; tt[2,3]=0.0; tt[2,4]=0.0; tt[2,5]=0.0
    tt[3,0]=0.0; tt[3,1]=0.0; tt[3,2]=0.0; tt[3,3]=  c; tt[3,4]=  s; tt[3,5]=0.0
    tt[4,0]=0.0; tt[4,1]=0.0; tt[4,2]=0.0; tt[4,3]= -s; tt[4,4]=  c; tt[4,5]=0.0
    tt[5,0]=0.0; tt[5,1]=0.0; tt[5,2]=0.0; tt[5,3]=0.0; tt[5,4]=0.0; tt[5,5]=1.0
    return ek,tt


def bfvec_grd(gamma,A,x1,y1,x2,y2):
    # Distributed vertical load vector
    bfv=np.zeros(6,dtype=np.float64)
    el=np.sqrt((x2-x1)**2+(y2-y1)**2)
    gkv=-1.0 # downward direction
    bfv[0]=0.0
    bfv[1]=0.0
    bfv[2]=0.5*gamma*A*el*gkv
    bfv[3]=0.0
    bfv[4]=0.0
    bfv[5]=0.5*gamma*A*el*gkv
    return bfv


def main():
    # Main routine
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=2              # Number of nodes per element
    nfree=3            # Degree of freedom per node
    # data input and input data print
    npoin,nele,nsec,npfix,nlod,ae,node,x,mpfix,rdis,fp=inpdata_grd(fnameR,nod,nfree)
    prinp_grd(fnameW,npoin,nele,nsec,npfix,nlod,ae,x,fp,mpfix,rdis,node)
    # assembly of stiffness matrix & load vectors
    ir    =np.zeros(nod*nfree,dtype=np.int)             # Work vector for matrix assembly
    gk    =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64) # Global stiffness matrix
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        aa=ae[2,m]
        ai=ae[3,m]
        aj=ae[4,m]
        gamma=ae[5,m]
        A=aa
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)      # Stiffness matrix in local coordinate
        ck   =np.dot(np.dot(tt.T,ek),tt)     # Stiffness matrix in global coordinate
        bfe  =bfvec_grd(gamma,A,x1,y1,x2,y2) # Body force vector downward direction
        ir[5]=3*j+2; ir[4]=ir[5]-1; ir[3]=ir[4]-1
        ir[2]=3*i+2; ir[1]=ir[2]-1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+ck[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of section force
    fsec=np.zeros([nod*nfree,nele],dtype=np.float64) # Section force vector
    dis =np.zeros(nod*nfree,dtype=np.float64) # work vector for section force calculation
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        m=node[2,ne]-1
        x1=x[0,i]
        y1=x[1,i]
        x2=x[0,j]
        y2=x[1,j]
        ee=ae[0,m]
        po=ae[1,m]
        ai=ae[3,m]
        aj=ae[4,m]
        GJ=0.5*ee/(1.0+po)*aj
        EI=ee*ai
        ek,tt=sm_grd(GJ,EI,x1,y1,x2,y2)
        dis[0]=disg[3*i]; dis[1]=disg[3*i+1]; dis[2]=disg[3*i+2]
        dis[3]=disg[3*j]; dis[4]=disg[3*j+1]; dis[5]=disg[3*j+2]
        fsec[:,ne]=np.dot(ek,np.dot(tt,dis))
    # print out of result
    prout_grd(fnameW,npoin,nele,disg,fsec)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()


#==============
# Execution
#==============
if __name__ == '__main__': main()