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

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

5. 固體模型

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

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

這里的梁采用小變形條件下的歐拉梁。根據(jù)歐拉梁理論,梁截面為剛性面且保持與中性線垂直。圖9所示為一梁單元。梁表面上有一流體網(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)矩陣。對于一般情況,Φ矩陣的分量是歐拉轉(zhuǎn)角的函數(shù),在小變形小轉(zhuǎn)動條件下,可以簡化為O點(diǎn)處梁截面的轉(zhuǎn)角的函數(shù):

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

圖9

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

圖10

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

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

圖11

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

6. UDF

UDF是用Fluent做流固耦合的關(guān)鍵。UDF是Fluent 用來提供可擴(kuò)展功能的框架。做過windows或者linux程序開發(fā)的朋友會覺得很好理解,但不熟悉編程的朋友可能會覺得費(fèi)解。在這里我簡單解釋一下它的意思。Fluent 是個編譯好的程序,想要實(shí)現(xiàn)功能擴(kuò)展,最方便的辦法就是利用動態(tài)鏈接庫。在程序運(yùn)行的時候?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ù),去運(yùn)行用戶定義的代碼。所以對用戶來說,就是填寫一個指定的函數(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)格的運(yùn)動。

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

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

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

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

|-> ultra_64 (這個可能根據(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所示。每個網(wǎng)格的數(shù)據(jù)點(diǎn)是中心點(diǎn) (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é)點(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;

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

// 這一段用來取得梁表面的流體網(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)

{

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

}

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

{

node=F_NODE(f,t,idNodeL);

// NODE_X, NODE_Y, NODE_Z存儲了節(jié)點(diǎn)坐標(biāo),但是注意node不是整數(shù)類型,實(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è)置為動網(wǎng)格

#if !RP_HOST // SERIAL OR NODE

SET_DEFORMING_THREAD_FLAG(THREAD_T0(t));

#endif

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

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

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

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

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

}

6.2 力的計算

梁表面所受到的激勵分為法向力和切向力。法向力由表面壓力引起;切向力就是剪力。梁的表面覆蓋著流體網(wǎng)格的面(face)。數(shù)據(jù)的計算在面中點(diǎn)(face centroid)上進(jìn)行。法向力的大小等于壓力乘以網(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實(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時間積分。注意Fortran 函數(shù)的末尾要加上一個額外的下劃線。另外變量名前要加&,因?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)的位置?？梢圆捎肗V_V命令賦值。NODE_COORD(node)對應(yīng)的是節(jié)點(diǎn)node的坐標(biāo)。由于節(jié)點(diǎn)循環(huán)的時候會重復(fù)訪問兩個單元共有的節(jié)點(diǎn),因此更新的時候需要用NODE_POS_NEED_UPDATE檢查當(dāng)前節(jié)點(diǎn)是否已經(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次以上。這樣才能看出周期運(yùn)動。圖13是結(jié)構(gòu)網(wǎng)格變形的例子。

圖13

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

圖14

8. 后記

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

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

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

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

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