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

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

5. 固體模型

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

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

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

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

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

圖9

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

圖10

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

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

圖11

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

6. UDF

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖12

6.1 UDF - DEFINE_GRID_MOTION

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

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

#include "udf.h"

#include "dynamesh_tools.h"

#include "sg.h"

// 定義計(jì)算梁響應(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),存儲(chǔ)了有關(guān)網(wǎng)格的信息,但不是網(wǎng)格本身

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

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

// dtime - 時(shí)間步長(zhǎng)

// 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;

// 定義計(jì)算結(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é)點(diǎn)間數(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;

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

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

#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)

{

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

}

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

{

node=F_NODE(f,t,idNodeL);

// NODE_X, NODE_Y, NODE_Z存儲(chǔ)了節(jié)點(diǎn)坐標(biāo),但是注意node不是整數(shù)類(lèi)型,實(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é)點(diǎn) ... (省略)

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

#if !RP_HOST // SERIAL OR NODE

SET_DEFORMING_THREAD_FLAG(THREAD_T0(t));

#endif

// 計(jì)算表面力(見(jiàn)6.2節(jié))

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

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

// 將數(shù)據(jù)傳輸?shù)接?jì)算節(jié)點(diǎn) ... (省略)

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

}

6.2 力的計(jì)算

梁表面所受到的激勵(lì)分為法向力和切向力。法向力由表面壓力引起;切向力就是剪力。梁的表面覆蓋著流體網(wǎng)格的面(face)。數(shù)據(jù)的計(jì)算在面中點(diǎn)(face centroid)上進(jìn)行。法向力的大小等于壓力乘以網(wǎng)格面積。剪力等于剪應(yīng)力乘以網(wǎng)格面積。剪應(yīng)力等于粘度乘以速度的導(dǎo)數(shù)。速度的導(dǎo)數(shù)可以簡(jiǎn)單近似為為cell centre velocity 除以cell centre 到 wall的距離。當(dāng)然也可以采用更高階的近似方法。需要注意的是在并行系統(tǒng)上,網(wǎng)格被分為多個(gè)區(qū),每個(gè)區(qū)之間的交界面上的face會(huì)被重復(fù)計(jì)算。為了防止這種情況發(fā)生,需要用PRINCIPAL_FACE_P判斷是否為該face實(shí)際存在于當(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時(shí)間積分。注意Fortran 函數(shù)的末尾要加上一個(gè)額外的下劃線。另外變量名前要加&,因?yàn)镕ortran函數(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é)點(diǎn)的位置。可以采用NV_V命令賦值。NODE_COORD(node)對(duì)應(yīng)的是節(jié)點(diǎn)node的坐標(biāo)。由于節(jié)點(diǎn)循環(huán)的時(shí)候會(huì)重復(fù)訪問(wèn)兩個(gè)單元共有的節(jié)點(diǎn),因此更新的時(shí)候需要用NODE_POS_NEED_UPDATE檢查當(dāng)前節(jié)點(diǎn)是否已經(jīng)被更新過(guò)。

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. 計(jì)算結(jié)果

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

圖13

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

圖14

8. 后記

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

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

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

朋友們好好學(xué),這年頭知識(shí)就是錢(qián)

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