UEL單元開發教程 | 三節點梁單元

易公子

2022年11月26日 16:25

瀏覽：2975 評論：2

知乎、B站：[易木木響叮當]
關注可了解更多的有限元數值仿真技巧。問題或建議，請公眾號留言;
如果你覺得木木同學對你有幫助，歡迎贊賞。

本次給大家帶來的主要內容是：如何使用UEL子程序中開發三節點平面梁單元？

關于梁單元的介紹，我在之前的推文中介紹過如何用Matlab編制梁單元有限元程序，及兩節點的梁單元UEL子程序編制，感興趣的朋友可以點擊跳轉閱讀瀏覽。

要點

為防止軸力過大出現“過度的剛性行為”，常常在梁單元內部增加一個節點，發揮“緩和”作用，增加的這個節點只有切線自由度，即切線位移，單元剛度矩陣的求解需要用到高斯數值積分。

如上圖所示，三節點梁單元共有 7 個自由度，需要兩個高斯積分點保證數值結果精度，節點順序按照上圖表示的來，中間節點為第三個節點。

有限元格式

為應變關系矩陣，為軸向應變，為曲率，軸向力，力矩，單元長度。應變矩陣：

局部坐標系下的B矩陣為：

坐標轉換矩陣：

通過坐標轉換矩陣將局部坐標系下的轉換為整體坐標系下的矩陣。本構關系：

單元剛度矩陣：

單元內力矢量：

代碼解讀

dx=coords(1,2)-coords(1,1)
dy=coords(2,2)-coords(2,1)
dl2=dx**2+dy**2
dl=sqrt(dl2)fortran
hdl=dl/two
acos=dx/dl
asin=dy/dl

計算單元長度，單元轉角和；
coords(1,2)表示第 2 個節點的坐標，相應的coords(1,1)表示第 1 個節點的坐標。

b(1,1)=(-three+four*g)*acos/dl
b(1,2)=(-three+four*g)*asin/dl
b(1,3)=zero
b(1,4)=(-one+four*g)*acos/dl
b(1,5)=(-one+four*g)*asin/dl
b(1,6)=zero
b(1,7)=(four-eight*g)/dl
b(2,1)=(-six+twelve*g)*(-asin)/dl2
b(2,2)=(-six+twelve*g)*acos/dl2
b(2,3)=(-four+six*g)/dl
b(2,4)=(six-twelve*g)*(-asin)/dl2
b(2,5)=(six-twelve*g)*acos/dl2
b(2,6)=(-two+six*g)/dl
b(2,7)=zero

以上代碼用于計算應變關系矩陣，由于Fortran語言不方便矩陣的乘積運算，可直接運算好寫入程序。

如上面代碼片所示，還可以使用我之前做的Fortran矩陣運算相關推文里面的矩陣相乘子程序，效果是一樣的。

eps=zero ! 當前應變
deps=zero ! 應變增量
cap=zero ! 當前曲率
dcap=zero ! 曲率增量
do k=1,7
    eps=eps+b(1,k)*u(k) 
    deps=deps+b(1,k)*du(k,1)
    cap=cap+b(2,k)*u(k)
    dcap=dcap+b(2,k)*du(k,1)
end do

上式是對應變矩陣的編譯；

對 7 個自由度進行遍歷，u(k)表示單元第個自由度的值;

du(k,1)表示第個自由度的增量值，有關du的解釋，官方培訓手冊給出如下解讀，一般情況下KRHS取1即可

    isvint=1+(kintk-1)*nsvint
    bn=zero ! 軸力
    bm=zero ! 彎矩
    daxial=zero ! 切線軸向剛度
    dcoupl=zero ! 切線剛度耦合項
    call ugenb(bn,bm,daxial,dbend,dcoupl,eps,deps,cap,dcap,
1              svars(isvint),nsvint,props,nprops)

isvint：在當前積分點獲取SDV；
使用ugenb本構關系子程序獲得軸向應變、曲率、軸向力和彎矩，存儲在高斯積分點上。

         do k1=1,7
           rhs(k1,1)=rhs(k1,1)-hdl*(bn*b(1,k1)+bm*b(2,k1))
           bd1=hdl*(daxial*b(1,k1)+dcoupl*b(2,k1))
           bd2=hdl*(dcoupl*b(1,k1)+dbend*b(2,k1))
           do k2=1,7
             amatrx(k1,k2)=amatrx(k1,k2)+bd1*b(1,k2)+bd2*b(2,k2)
           end do
         end do
       end do

計算單元內力矢量：，目前階段的UEL分析中不涉及等效節點荷載處理的情況，均為0；
離散化：

，為高斯積分系數。

 subroutine ugenb(bn,bm,daxial,dbend,dcoupl,eps,deps,cap,dcap,
1                 svint,nsvint,props,nprops)
 include 'aba_param.inc'
 parameter(zero=0.d0,twelve=12.d0)
 dimension svint(*),props(*)
 h=props(1)! 截面高度
 w=props(2)! 截面寬度
 E=props(3)! 彈性模量
 daxial=E*h*w
 dbend=E*w*h**3/twelve
 dcoupl=zero
 bn=svint(1)+daxial*deps
 bm=svint(2)+dbend*dcap
 svint(1)=bn
 svint(2)=bm
 svint(3)=eps
 svint(4)=cap
 return
 end