2節點梁單元的彈性分析(python,有限元)

import numpy as np

import A

nels=4

np_types=2

loaded_nodes=4

fixed_freedoms=2

nod=2

nprops=1

nn=nels+1

nr=2

nodof=2

ndof=nod*nodof

ell=np.array([[2.5],[2.5],[3.0],[2.0]])

# 初始化定義陣列

nf=np.ones((nodof,nn),dtype=np.int64)

g=np.ones((ndof,1),dtype=np.int64)

num=np.ones((nod,1),dtype=np.int64)

etype=np.ones((nels,1),dtype=np.int64)

eld=np.ones((ndof,1))

km=np.ones((ndof,ndof))

mm=np.zeros((ndof,ndof))

action=np.ones((ndof,1))

g_g=np.ones((ndof,nels))

prop=np.ones((nprops,np_types))

if np_types==1:

   etype[:,0]=1

else:

   etype[:,0]=(1,1,2,2)

if np_types==1:

   prop[:nprops,np_types-1]=5e5

else:

   prop[0,:]=(4.0e4,2.0e4)

# 讀取 nr

if nr!=0:

   Dim_1=[1,4]

   nf_value=np.array([[0,0],[1,1]])

   m=0

   for i in Dim_1:

       for j in range(1,nodof+1):

         nf[j-1,i-1]=nf_value[j-1,m]

       m=m+1  

#form nf

A.formnf(nf)

neq=int(max(nf.reshape(nf.shape[0]*nf.shape[1],1)))

kdiag=np.zeros((neq,1),dtype=np.int64)

loads=np.zeros((neq+1,1))

## 累加單元去發現全域性尺寸

for iel in range(1,nels+1):

   num[:nod,0]=(iel,iel+1)

   A.num_to_g(num,nf,g)

   g_g[:,iel-1]=g[:,0]

#call fkidiag

   A.fkdiag(kdiag,g)

print(neq)

for i in range(1,neq):

   kdiag[i]=kdiag[i]+kdiag[i-1]

kv=np.zeros((kdiag[neq-1,0],1),dtype=float)

print(' 一共有 '+str(neq)外匯跟單gendan5.com+' 等式和天際線儲存個數為 '+str(kdiag[neq-1]))

# 全域性剛度矩陣安裝

for iel in range(1,nels+1):

#call pin_jointed

   A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])

   g[:,0]=g_g[:,iel-1]

   if nprops>1:

     A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])

#call fsparv

   A.fsparv(kv,km+mm,g,kdiag)

### 讀荷載和位移

if loaded_nodes!=0:

   Dim_2 = [2,3,4,5]

   load_value=np.array([[-20.0,-6.0,-8.8,-1.2],[0.0,-3.0,2.2,0.5333]])

   m=0

   for i in Dim_2:

       for j in range(1,nodof+1):

         loads[nf[j-1,i-1]-1]=load_value[j-1,m]

       m=m+1

if fixed_freedoms!=0:

   node=np.ones((fixed_freedoms,1),dtype=np.int64)

   sense=np.ones((fixed_freedoms,1),dtype=np.int64)

   value=np.ones((fixed_freedoms,1))

   no=np.ones((fixed_freedoms,1),dtype=np.int64)

   node[:,0]=(1,3)

   sense[:,0]=(2,1)

   value[:,0]=(-0.001,-0.005)

#### 位移賦值

   for i in range(1,fixed_freedoms+1):

     no[i-1]=nf[sense[i-1]-1,node[i-1]-1]

   kv[kdiag[no-1]-1]=kv[kdiag[no-1]-1]+1e20

   loads[no-1,0]=kv[kdiag[no-1,0]-1,0]*value

##### 方程求解

#call sparin  

A.sparin(kv,kdiag)

##call spabac

A.spabac(kv,loads,kdiag)

loads[neq]=0

print(' 節點 位移    轉角 ')

for k in range(1,nn+1):

   print(k,end=' ')

   for m in range(1,nodof+1):

     print('{:9.4e}'.format(loads[nf[m-1,k-1]-1,0]),end=' ')

   print( )

## 重新取得單元和節點力

print(' 單元  力       力矩       力       力矩 ')

axial=[0]

for iel in range(1,nels+1):

#call pin_jointed

     #call pin_jointed

   A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])

   g[:,0]=g_g[:,iel-1]

   if nprops>1:

     A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])

##########

   eld=loads[g-1,0]

   action=np.dot(km+mm,eld)

   print(iel,end=' ')

   for i in range(1,ndof+1):

     print('{:9.4e}'.format(action[i-1,0]),end=' ')

   print( )