Fluent流固耦合基礎(chǔ)教程(下)

2016-09-04  by:CAE仿真在線  來源:互聯(lián)網(wǎng)

5. 固體模型

        固體變形和運動的求解可以采用很多方法。對于簡單的問題,可以采用解析解。目前的這個彈性梁的問題,就可以用采用解析解的方法。對于復(fù)雜結(jié)構(gòu),有限元是比較常用的方法。建立有限元模型有兩種方案。一種是采用三維實體單元完整地模擬出整個梁。這種方法的好處是流體和固體的接觸面兩邊都是實際網(wǎng)格,一邊是流體網(wǎng)格,一邊是有限元網(wǎng)格。邊界上的位移,速度,受力等計算都比較簡單,缺點是固體模型計算量大。第二種方法是采用結(jié)構(gòu)體單元,如梁,板殼等。固體模型簡單,但是接觸面上的處理稍微復(fù)雜一些。結(jié)構(gòu)體單元在工程上應(yīng)用得比較普遍。


        為了說明如何利用結(jié)構(gòu)體單元做流固耦合,這里采用三節(jié)點的梁單元。梁單元是一維單元,單元的幾何形狀只是一條中性線,沒有體積。梁的變形是由中性線的位移和梁截面的轉(zhuǎn)動描述的。在流體模型中,梁的體積是存在的。梁表面上的流體網(wǎng)格節(jié)點的位置需要由梁的變形確定。因此一個關(guān)鍵的步驟是根據(jù)梁的變形計算出梁表面的流體網(wǎng)格節(jié)點的位置。


        這里的梁采用小變形條件下的歐拉梁。根據(jù)歐拉梁理論,梁截面為剛性面且保持與中性線垂直。圖9所示為一梁單元。梁表面上有一流體網(wǎng)格節(jié)點A。點A所在的截面軸向坐標(biāo)為ξ。則流體網(wǎng)格節(jié)點的位置可以表示為:

        其中uA是A點的位移向量,uO是O點的位移向量。r是O點到A點的位置向量。Φ是旋轉(zhuǎn)矩陣。對于一般情況,Φ矩陣的分量是歐拉轉(zhuǎn)角的函數(shù),在小變形小轉(zhuǎn)動條件下,可以簡化為O點處梁截面的轉(zhuǎn)角的函數(shù):

        O點處的位移和梁截面的轉(zhuǎn)角可以通過有限元的形函數(shù)求得。

圖9

        圖10給出了梁變形后其表面上流體網(wǎng)格節(jié)點的位置。兩個端點處的流體網(wǎng)格沒有加以控制,而是交給fluent自動處理,這只是對當(dāng)前這種支撐情況有效。對于一般情況,如懸臂梁,則需要人工計算。

圖10

        與網(wǎng)格位置計算類似。流體計算所得到的力需要傳遞給固體模型。梁表面的流體網(wǎng)格上的壓力和剪力也需要利用有限元的形函數(shù)離散到有限元節(jié)點上。有限元方面的基本知識,我在這里就不羅嗦了。任何一本教材上都寫得很清楚。


        流固耦合問題是瞬態(tài)問題,因此需進行時間積分。常用的時間積分算法有Newmark,Wilson-theta,Runge-Kutta等。Newmark方法在多自由度的線性問題中應(yīng)用比較普遍。這里我們也采用Newmark方法。但是在程序編寫的時候需要注意Fluent 求解流固耦合問題的流程。Fluent必須作為整個過程的主導(dǎo)程序,如圖11所示。

圖11

        這種流程模式給程序編寫帶來一些麻煩,因為Newmark算法的循環(huán)被拆解開進行,必須按照單個時間步考慮。這就要求每個時間步之間的數(shù)據(jù)傳遞不能用通常的變量存儲。解決的辦法有兩個:一個辦法是用硬盤存儲,但是這樣耽誤在文件I/O上的時間很多;另外一個辦法是利用常駐內(nèi)存的global變量存儲,但是具體操作要看所用的編程語言環(huán)境。這里我用Fortran90編寫的有限元程序,global 變量保存在獨立的module里。如果只用UDF的話,跟C語言一致。


6. UDF

        UDF是用Fluent做流固耦合的關(guān)鍵。UDF是Fluent 用來提供可擴展功能的框架。做過windows或者linux程序開發(fā)的朋友會覺得很好理解,但不熟悉編程的朋友可能會覺得費解。在這里我簡單解釋一下它的意思。Fluent 是個編譯好的程序,想要實現(xiàn)功能擴展,最方便的辦法就是利用動態(tài)鏈接庫。在程序運行的時候?qū)⒅付ǖ膭討B(tài)鏈接庫調(diào)入內(nèi)存,然后通過預(yù)先定義好的接口函數(shù)執(zhí)行用戶定義的子程序。換句話說UDF就是一些放在動態(tài)鏈接庫里的函數(shù)。這些函數(shù)的名字和格式已經(jīng)由Fluent預(yù)先定義好,但是內(nèi)容是空的,需要用戶來寫。Fluent 會在預(yù)定的情況下調(diào)用指定的函數(shù),去運行用戶定義的代碼。所以對用戶來說,就是填寫一個指定的函數(shù),編譯成動態(tài)鏈接庫,在Fluent里鏈接上,然后等著Fluent來調(diào)用。如此簡單。


        這些預(yù)先指定好的函數(shù)由被稱為定義宏(DEFINE macro)的宏命令來定義。這是個很拗口的說法,不過一般用戶不必深究,拿它當(dāng)函數(shù)來用就是。為了簡單起見,下面用“UDF函數(shù)”來代替“定義宏”。 Fluent 提供了很多UDF函數(shù),具體有哪些,都是什么功能,諸位可以看Fluent UDF手冊。這里就只說跟流固耦合有關(guān)的一個。下面的UDF函數(shù)用來定義網(wǎng)格的運動。

DEFINE_GRID_MOTION( name, d, dt, time, dtime)


        除了編譯的方法,Fluent 還支持對UDF的逐行解釋執(zhí)行。這樣做方便調(diào)試,但是動網(wǎng)格的一些UDF函數(shù)是不能這樣做的,所以我們還得用動態(tài)鏈接庫。


        UDF是C語言編寫的。Fluent 自帶一個C編譯器和一堆頭文件。UDF的編譯可以在Fluent GUI里自動進行,但是也可以手工完成。有些情況下必須要手工操作,比如有C/Fortran混合編程的時候。順便說一句,混合程序的編譯也是個頭疼的問題,需要費很多周章。這跟所用的操作系統(tǒng)和Fortran 版本有關(guān)。如果所寫的UDF比較簡單,只包括普通的C文件,則自動編譯就可以。如果需要手工操作,則要建立如下文件結(jié)構(gòu):

work folder (工作目錄,模型放這里)

  |-> libudf (UDF目錄,有個Makefile)

  |-> src (UDF源文件放這里)

  |-> ultra_64 (這個可能根據(jù)所用操作系統(tǒng)有所不同)

  |-> 3ddp (3ddp單機版)

  |-> 3ddp_host(3ddp并行版主節(jié)點)

  |-> 3ddp_node(3ddp并行版從節(jié)點)

然后修改src/makefile文件。之后回到libudf目錄執(zhí)行make。


        Fluent 的數(shù)據(jù)結(jié)構(gòu)是cell - face - node,如圖12所示。每個網(wǎng)格的數(shù)據(jù)點是中心點 (cell center)。

圖12


6.1 UDF - DEFINE_GRID_MOTION

這一節(jié)講一下DEFINE_GRID_MOTION 的大體結(jié)構(gòu)。

// 首先引用幾個頭文件。這幾個頭文件在Fluent的src 目錄下。

#include "udf.h"

#include "dynamesh_tools.h"

#include "sg.h"


// 定義計算梁響應(yīng)的有限元函數(shù)。此函數(shù)是Fortran函數(shù)。

extern void beam_response_(

  int *nfgrid,

  double (*fgrid)[ARR_LIMIT_FACE_BEAM][3],

  double (*force)[ARR_LIMIT_FACE_BEAM][3],

  double *gamm,

  double *beta,

  double *dt,

  double *t,

  int *run,

  int *ndgrid,

  double (*dgrid)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  double (*disp)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  double (*velo)[ARR_LIMIT_FACE_BEAM*NNODE_ON_FACE][3],

  int *info

);


// ... (定義其它全局變量和函數(shù))


/* ------------------------------------------------------------------------- */

/* --------------------- macros are defined below here --------------------- */

/* ------------------------------------------------------------------------- */


/* Macro for defining dynamic mesh udf for the beam */

DEFINE_GRID_MOTION(beam,domain,dt,time,dtime)

{

  // beam - UDF 的名字

  // domain - 當(dāng)前的domain數(shù)據(jù)結(jié)構(gòu),存儲了有關(guān)網(wǎng)格的信息,但不是網(wǎng)格本身

  // dt - 指向網(wǎng)格數(shù)據(jù)結(jié)構(gòu)(thead)的指針

  // time - 當(dāng)前時刻

  // dtime - 時間步長


  // Define variables: pointers and utilities

  Thread *t = DT_THREAD(dt);

  cell_t c;

  face_t f;

  Node *node;

  int idNodeL; // local index of a node inside a face

  int countF; // number of faces

  int countN; // number of nodes

  int index;

  int i,j,k;

  int run;

  int id;


  // 定義計算結(jié)構(gòu)響應(yīng)需要的變量

  double xi; // the axial coordinate of a grid node

  double area; // area

  double pressure; // pressure

  double distance; // distance between two points

  double a_by_es; // A'A/(A'e)

  double gamm=0.5;

  double beta=0.25;


  // 定義主/從節(jié)點間數(shù)據(jù)傳遞用的數(shù)組... (省略)


  // Define constants

  const double PI=3.14159265358979323846;

  const double VISCOSITY = 0.001;

  const double DENSITY = 1000.0;

  const double TOLERANCE = 1e-15;

  const double LOWERLIMIT = 1e-8; 


  // 向量初始化 ... (省略)


  // 這一段用來取得梁表面的流體網(wǎng)格節(jié)點位置

  #if !RP_HOST // SERIAL OR NODE

  countF=0;

  countN=0;

  begin_f_loop(f,t) // 掃描當(dāng)前domain的所有face

  {

    if PRINCIPAL_FACE_P(f,t)

    {

      countF++;

      if(countF>ARR_LIMIT_FACE_BEAM)

      {

        // 錯誤處理 ... (省略)

      }

      f_node_loop(f,t,idNodeL) // 掃描當(dāng)前face的所有node

      {

         node=F_NODE(f,t,idNodeL);

         // NODE_X, NODE_Y, NODE_Z存儲了節(jié)點坐標(biāo),但是注意node不是整數(shù)類型,實際上是聯(lián)合變量。

         arrNode[countN][0]=NODE_X(node);

         arrNode[countN][1]=NODE_Y(node);

         arrNode[countN][2]=NODE_Z(node);

         countN++;

       }

    }

  }

  end_f_loop(f,t);

  #endif


  // 將數(shù)據(jù)傳輸?shù)街鞴?jié)點 ... (省略)


  // 將相鄰區(qū)域設(shè)置為動網(wǎng)格

  #if !RP_HOST // SERIAL OR NODE

  SET_DEFORMING_THREAD_FLAG(THREAD_T0(t));

  #endif


  // 計算表面力(見6.2節(jié))


  // 將數(shù)據(jù)傳輸?shù)街鞴?jié)點 ... (省略)


  // 計算梁的響應(yīng)(見6.3節(jié))


  // 將數(shù)據(jù)傳輸?shù)接嬎愎?jié)點 ... (省略)


  // 更新網(wǎng)格 (見6.4節(jié))


}


6.2 力的計算

梁表面所受到的激勵分為法向力和切向力。法向力由表面壓力引起;切向力就是剪力。梁的表面覆蓋著流體網(wǎng)格的面(face)。數(shù)據(jù)的計算在面中點(face centroid)上進行。法向力的大小等于壓力乘以網(wǎng)格面積。剪力等于剪應(yīng)力乘以網(wǎng)格面積。剪應(yīng)力等于粘度乘以速度的導(dǎo)數(shù)。速度的導(dǎo)數(shù)可以簡單近似為為cell centre velocity 除以cell centre 到 wall的距離。當(dāng)然也可以采用更高階的近似方法。需要注意的是在并行系統(tǒng)上,網(wǎng)格被分為多個區(qū),每個區(qū)之間的交界面上的face會被重復(fù)計算。為了防止這種情況發(fā)生,需要用PRINCIPAL_FACE_P判斷是否為該face實際存在于當(dāng)前的區(qū)里。


// Obtain the centroid location and the force on the centroids

index=0;

begin_f_loop(f,t)

{

  if PRINCIPAL_FACE_P(f,t)

  {

  // Save the face id

  arrFaceID[index]=f;

  // Obtain the centroid coordinates and save in arrCentr

  F_CENTROID(vecXf,f,t);

  arrCentr[index][0]=vecXf[0]; // x

  arrCentr[index][1]=vecXf[1]; // y

  arrCentr[index][2]=vecXf[2]; // z

  // Obtain the area vector and the area

  F_AREA(vecArea,f,t);

  area=sqrt(pow(vecArea[0],2)+pow(vecArea[1],2)+pow(vecArea[2],2));

  // Obtain the pressure and calculate the pressure force

  pressure=F_P(f,t);

  vecLift[0]=vecArea[0]*pressure;

  vecLift[1]=vecArea[1]*pressure;

  vecLift[2]=vecArea[2]*pressure;

  arrPForce[index][0]=vecLift[0];

  arrPForce[index][1]=vecLift[1];

  arrPForce[index][2]=vecLift[2];

  // Obtain the shear stress and calculate the viscous force

  // get the face parameters

  BOUNDARY_FACE_GEOMETRY(f,t,vecArea,distance,vecEs,a_by_es,vecDr0)

  if(distance < TOLERANCE) // distance is always positive

  {

    distance=LOWERLIMIT;

  }

  // -- get the centroid velocity of the cell

  vecVelo[0]=C_U(c,t);

  vecVelo[1]=C_V(c,t);

  vecVelo[2]=C_W(c,t);

  // -- calculate the viscous force (= shear stress * area)

  vecDrag[0]=(VISCOSITY/distance)*vecVelo[0]*area;

  vecDrag[1]=(VISCOSITY/distance)*vecVelo[1]*area;

  vecDrag[2]=(VISCOSITY/distance)*vecVelo[2]*area;

  arrVForce[index][0]=vecDrag[0];

  arrVForce[index][1]=vecDrag[1];

  arrVForce[index][2]=vecDrag[2];

  // Increase index

  index=index+1;

  }

}

end_f_loop(f,t);

// Calculate the total force

for(i=1;i<=countF;i++)

{

  arrForce[1] = arrPForce[1] + arrVForce[1];

  arrForce[2] = arrPForce[2] + arrVForce[2];

  arrForce[3] = arrPForce[3] + arrVForce[3];

}


6.3結(jié)構(gòu)響應(yīng)

結(jié)構(gòu)響應(yīng)由Fortran函數(shù)求得,采用Newmark時間積分。注意Fortran 函數(shù)的末尾要加上一個額外的下劃線。另外變量名前要加&,因為Fortran函數(shù)參數(shù)都是地址傳遞而非值傳遞。

// Perform the Newmark scheme.

if (flagInitial!=1)

{

  // Calculate the beam response

  info=0;

  for (i=1;i<=countN;i++)

  {

    arrDisp[0]=0.0;

    arrDisp[1]=0.0;

    arrDisp[2]=0.0;

    arrVelo[0]=0.0;

    arrVelo[1]=0.0;

    arrVelofor(i=1;i<=countF;i++)[2]=0.0;

  }

  run=1;

  beam_response_(

  &countF, // nfgrid,

  &arrCentr, // fgrid,

  &arrForce,        // force,

  &gamm, // gamma=0.5,

  &beta, // beta=0.25,

  &dtime, // dt,

  &time, // t,

  &run,         // run, 0-preparation only, 1-run

  &countN, // ndgrid,

  &arrNode, // dgrid,

  &arrDisp, // disp,

  &arrVelo,         // velo,

  &info // info

  );

}


6.4網(wǎng)格更新

得到結(jié)構(gòu)響應(yīng)以后,需要更新梁表面流體網(wǎng)格節(jié)點的位置??梢圆捎肗V_V命令賦值。NODE_COORD(node)對應(yīng)的是節(jié)點node的坐標(biāo)。由于節(jié)點循環(huán)的時候會重復(fù)訪問兩個單元共有的節(jié)點,因此更新的時候需要用NODE_POS_NEED_UPDATE檢查當(dāng)前節(jié)點是否已經(jīng)被更新過。



if (flagInitial!=1)

{

  // Update the grid nodes

  index=0;

  begin_f_loop(f,t)

  {

    if PRINCIPAL_FACE_P(f,t)

    {

      f_node_loop(f,t,idNodeL)

      {

        node=F_NODE(f,t,idNodeL); NV_V(NODE_COORD(node)

        // Update node if the current node has not been previously

        // visited when looping through previous faces

        if (NODE_POS_N{

  arrForce[1] = arrPForce[1] + arrVForce[1];

  arrForce[2] = arrPForce[2] + arrVForce[2];

  arrForce[3] = arrPForce[3] + arrVForce[3];

}


7. 計算結(jié)果

        為了保證模型的正確性,需要查看數(shù)值模擬的結(jié)果。這里不講如何詮釋結(jié)果,只講如何檢查流固耦合過程是否正確。最直接的方法是觀察網(wǎng)格的變形。簡單的做法是在Write/Autosave中定義自動保存的頻率。一般按照時間步來保存,比如選擇每1000個時間步保存一次。然后查看每次保存的時刻上網(wǎng)格的變形程度。保存的時候一定要將case和data文件都保存下來。另外保存頻率的選擇要保證一個振動周期內(nèi)保存至少5次以上。這樣才能看出周期運動。圖13是結(jié)構(gòu)網(wǎng)格變形的例子。

圖13

        查看流體網(wǎng)格是比較麻煩的工作,但也是必須的工作。如果只查看固體變形的過程,并不能說明流場網(wǎng)格也正確變化。筆者第開始運行程序的時候就遇到這樣的情況。;流體網(wǎng)格沒有更新,但是流體力被正確地傳遞給,而固體振動也是正常的。這樣的結(jié)果得到是接耦合的流致振動解,比如湍流激勵引起的振動。但是這樣的解是不正確的,因為沒有包括附加質(zhì)量,流體阻尼,和附加剛度。所以一個很重要的方法是觀察結(jié)構(gòu)的響應(yīng),然后根據(jù)響應(yīng)分析這些附加質(zhì)量等效果是否存在。對于這個梁,沒有附加質(zhì)量的自由振動周期為0.4秒左右,帶有附加質(zhì)量以后上升到0.57秒左右。與理論估計得到的結(jié)果一致。

圖14


8. 后記

        用每天擠牙膏的方式,終于寫完了。水平有限,時間有限,只能寫到這個層次了。這篇文章寫得很粗略,只是作為大概的介紹,很多細致的內(nèi)容還要讀者自己去探索。流固耦合還是個很有意思的話題的?,F(xiàn)在ANSYS做得好,恐怕大家以后都不用這么費勁寫UDF了。耦合問題現(xiàn)在主要需要面對的是計算量問題。解決實際工程問題,沒有并行計算恐怕很難完成。


        有限體積法和有限元已經(jīng)都是很成熟的東西。關(guān)鍵是湍流問題不好解決。大渦理論是個很好的方向,只是near wall條件不好處理。如果能做出好的wall function 就可以極大降低計算量。流固耦合的難題還是在強耦合問題上,據(jù)說ADINA做得不錯。


        希望這篇文章對打算進入流固耦合領(lǐng)域的朋友有點啟發(fā),起碼算是個快速入門。這篇文章里的代碼是具有代表性的段落,但不是完整的程序,請不要不加思考地照搬。請不要找我要完整代碼,我是不會給你的。做這個東西,我是有錢掙的。坦白地說一共掙了八千。如果你非要完整代碼也可以,我給你打到極限折,只收十分之一。順便說一句,貨幣單位是美元。


        朋友們好好學(xué),這年頭知識就是錢



開放分享:優(yōu)質(zhì)有限元技術(shù)文章,助你自學(xué)成才

相關(guān)標(biāo)簽搜索:Fluent流固耦合基礎(chǔ)教程(下) Fluent培訓(xùn) Fluent流體培訓(xùn) Fluent軟件培訓(xùn) fluent技術(shù)教程 fluent在線視頻教程 fluent資料下載 fluent分析理論 fluent化學(xué)反應(yīng) fluent軟件下載 UDF編程代做 Fluent、CFX流體分析 HFSS電磁分析 

編輯
在線報名:
  • 客服在線請直接聯(lián)系我們的客服,您也可以通過下面的方式進行在線報名,我們會及時給您回復(fù)電話,謝謝!
驗證碼

全國服務(wù)熱線

1358-032-9919

廣州公司:
廣州市環(huán)市中路306號金鷹大廈3800
電話:13580329919
          135-8032-9919
培訓(xùn)QQ咨詢:點擊咨詢 點擊咨詢
項目QQ咨詢:點擊咨詢
email:kf@1cae.com