網(wǎng)格實(shí)現(xiàn)翼型俯仰+尾緣變形完整攻略)
做風(fēng)力機(jī)葉片或者機(jī)翼的氣動(dòng)彈性分析時(shí)我經(jīng)常要面對(duì)一個(gè)不算特別復(fù)雜、但也非常容易翻車的需求翼型本身在繞某一點(diǎn)做俯仰振蕩與此同時(shí)尾緣還要疊加一定幅度的柔性變形。前者是典型的剛體運(yùn)動(dòng)對(duì)應(yīng)Fluent動(dòng)網(wǎng)格里的剛體區(qū)域加CG運(yùn)動(dòng)后者是典型的邊界變形對(duì)應(yīng)Deforming區(qū)域加網(wǎng)格光順。兩個(gè)拆開做教程一抓一大把都不難一旦要求它們同時(shí)作用在同一個(gè)翼型上很多人就開始犯嘀咕了——?jiǎng)泳W(wǎng)格區(qū)域到底怎么設(shè)UDF怎么寫同一個(gè)邊界能不能又剛體又變形網(wǎng)格會(huì)不會(huì)在某個(gè)時(shí)刻直接被擠爆……這篇文章就圍繞“Fluent動(dòng)網(wǎng)格實(shí)現(xiàn)翼型俯仰振蕩同時(shí)尾緣變形”這個(gè)組合把從方案選型、網(wǎng)格準(zhǔn)備、UDF編寫到求解器設(shè)置、常見問題排錯(cuò)的完整過程梳理一遍。案例本身選的是最簡(jiǎn)單的NACA0012翼型低速來(lái)流剛體俯仰加尾緣二次型變形但里面的方法論可以直接遷移到風(fēng)機(jī)葉片、渦輪葉片、直升機(jī)旋翼這類工程問題上。適合已經(jīng)會(huì)跑基本Fluent仿真、想進(jìn)階動(dòng)網(wǎng)格的工程師也適合正在做氣動(dòng)彈性課題、被“剛體柔性變形疊加”卡住的學(xué)生。1. 同樣都是動(dòng)網(wǎng)格為什么這個(gè)案例要單獨(dú)拿出來(lái)講1.1 先想清楚翼型在“動(dòng)”的到底是什么很多新手上來(lái)就寫UDF結(jié)果連自己要模擬的運(yùn)動(dòng)都沒拆清楚。翼型俯仰振蕩加尾緣變形聽起來(lái)是“一個(gè)運(yùn)動(dòng)”數(shù)學(xué)上其實(shí)是兩個(gè)獨(dú)立位移場(chǎng)的疊加。第一部分是剛體俯仰。整個(gè)翼型繞固定的彈性軸做正弦轉(zhuǎn)動(dòng)典型形式是θ(t) θ? A·sin(2πft)其中θ?是平均攻角A是俯仰振幅f是俯仰頻率。這個(gè)運(yùn)動(dòng)的特點(diǎn)是翼型表面上所有網(wǎng)格節(jié)點(diǎn)的相對(duì)位置不改變整體繞旋轉(zhuǎn)中心轉(zhuǎn)一個(gè)角度。在Fluent里最經(jīng)典的實(shí)現(xiàn)方式是DEFINE_CG_MOTION直接給定瞬時(shí)角速度即可。第二部分是尾緣變形。這個(gè)更微妙——它不是整體轉(zhuǎn)動(dòng)而是翼型尾緣附近的一段邊界相對(duì)它自身的初始位置做一個(gè)連續(xù)、光滑的偏移。比如尾緣點(diǎn)在某個(gè)時(shí)刻下移0.03倍弦長(zhǎng)然后離尾緣越近變形越大越往上游變形越小到某個(gè)起點(diǎn)位置變形恰好為零。這種運(yùn)動(dòng)不能用CG_MOTION表達(dá)因?yàn)樗茐牧恕皠傂浴薄獏^(qū)域內(nèi)節(jié)點(diǎn)之間的相對(duì)距離變了。所以本質(zhì)上是剛體位移場(chǎng) 柔性變形場(chǎng)按時(shí)間同步疊加。搞清楚這一點(diǎn)后面所有技術(shù)選擇都順理成章。1.2 剛體俯仰和尾緣變形的“疊加”難點(diǎn)難點(diǎn)恰恰出在Fluent動(dòng)網(wǎng)格的框架上。在一個(gè)動(dòng)網(wǎng)格計(jì)算里一個(gè)邊界區(qū)域通常只允許一種運(yùn)動(dòng)屬性設(shè)成Rigid Body整個(gè)區(qū)域做剛性平移/旋轉(zhuǎn)你無(wú)法單獨(dú)挑出尾緣那一小塊讓它再多動(dòng)一點(diǎn)設(shè)成Deforming所有節(jié)點(diǎn)的位置都由用戶自己控制剛體運(yùn)動(dòng)也得自己寫到UDF里。也就是說你不能把同一個(gè)壁面既設(shè)成剛體又設(shè)成Deforming。那怎么辦直覺的解決辦法是把翼型一分為二前段到尾緣上游設(shè)為Rigid Body做俯仰尾緣單獨(dú)設(shè)為Deforming做變形。中間切一個(gè)interface。這個(gè)思路我最早也試過結(jié)果是交界面兩側(cè)的網(wǎng)格密度必須高度匹配否則插值誤差會(huì)直接污染壁面壓力場(chǎng)而且邊界層在interface處被硬生生切斷動(dòng)網(wǎng)格光順在這種位置經(jīng)常出現(xiàn)負(fù)體積。真實(shí)算下來(lái)穩(wěn)定性很差。后來(lái)我換了一種思路既然Deforming模式允許我直接控制節(jié)點(diǎn)位置那為什么不在UDF里同時(shí)寫上“剛性旋轉(zhuǎn)位移 尾緣變形位移”讓Fluent的光順?biāo)惴ㄈヌ幚砭W(wǎng)格內(nèi)部更新這個(gè)思路最終穩(wěn)定跑通了也是下面這篇文章所有內(nèi)容的核心。整個(gè)方案的代價(jià)是UDF要自己寫但換來(lái)的是對(duì)流場(chǎng)物理的更可控、對(duì)運(yùn)動(dòng)疊加的完全接管。2. 方案選型放棄CG Motion把整個(gè)翼型交給GRID_MOTION2.1 Fluent動(dòng)網(wǎng)格三種方法哪些能用哪些不能用Fluent動(dòng)網(wǎng)格大體有三類手段Smoothing、Remeshing、Overset。選型前建議先把適用范圍圈清楚。方法適用場(chǎng)景本案例可行性Smoothing光順邊界位移中小幅度網(wǎng)格拓?fù)洳蛔兒诵氖侄闻浜螪iffusion光順Remeshing局部重構(gòu)大位移、大轉(zhuǎn)動(dòng)局部網(wǎng)格拓?fù)渲亟ㄝo助手段大俯仰角時(shí)需要Overset重疊網(wǎng)格多體大位移、相互穿越可用但殺雞用牛刀交界面插值額外耗成本本案例的運(yùn)動(dòng)量級(jí)如果控制在工程界很常見的范圍內(nèi)——俯仰振幅±5°、尾緣變形3%~5%弦長(zhǎng)——那么Smoothing加局部Remeshing足夠。Diffusion-based光順比彈簧光順的魯棒性好得多尤其適合邊界做旋轉(zhuǎn)運(yùn)動(dòng)的情況因?yàn)樗鼤?huì)把變形量“均勻擴(kuò)散”到全場(chǎng)而不是像彈簧一樣在局部積累。真正需要糾結(jié)的并不是這三種方法的取舍而是邊界區(qū)域運(yùn)動(dòng)屬性怎么設(shè)。我最終使用的是Deforming區(qū)域加DEFINE_GRID_MOTION讓一個(gè)UDF同時(shí)控制剛體旋轉(zhuǎn)和尾緣變形內(nèi)部網(wǎng)格交給Diffusion光順吸收變形。2.2 本案例的動(dòng)網(wǎng)格區(qū)域設(shè)置在Fluent里按照下面的方式設(shè)置動(dòng)網(wǎng)格區(qū)域思路很清晰打開Dynamic Mesh啟用Smoothing和Remeshing在Dynamic Mesh Zones中將翼型壁面設(shè)置為DeformingMotion UDF指定為后面要寫的airfoil_pitch_deform遠(yuǎn)場(chǎng)邊界設(shè)置為Stationary保持固定計(jì)算域內(nèi)部不需要額外指定動(dòng)態(tài)區(qū)域光順?biāo)惴〞?huì)自動(dòng)把壁面的位移擴(kuò)散到整個(gè)網(wǎng)格。這里有個(gè)容易誤解的點(diǎn)Deforming區(qū)域指定的雖然是翼型壁面但Fluent在調(diào)用UDF得到壁面節(jié)點(diǎn)位移后會(huì)通過Smoothing算法把位移逐步傳遞到內(nèi)部網(wǎng)格節(jié)點(diǎn)。所以UDF里只需要管壁面節(jié)點(diǎn)的位置內(nèi)部網(wǎng)格怎么動(dòng)是求解器的事。Smoothing參數(shù)里我習(xí)慣把Diffusion參數(shù)調(diào)到1.5~2.0。Diffusion參數(shù)越大壁面附近的位移衰減越快網(wǎng)格更傾向于在遠(yuǎn)場(chǎng)吸收變形這對(duì)保護(hù)邊界層質(zhì)量非常關(guān)鍵。如果案例的俯仰角較大再把Remeshing里Minimum Length Scale和Maximum Length Scale設(shè)成當(dāng)?shù)鼐W(wǎng)格尺寸的0.5倍和2倍目標(biāo)偏斜度0.7左右局部網(wǎng)格壞了就讓Fluent自動(dòng)重構(gòu)。2.3 計(jì)算域、網(wǎng)格與邊界層準(zhǔn)備要點(diǎn)網(wǎng)格是整個(gè)動(dòng)網(wǎng)格算例的地基。翼型幾何本身可以用NACA0012標(biāo)準(zhǔn)型值點(diǎn)生成弦長(zhǎng)取1m。計(jì)算域推薦C型或O型拓?fù)溥h(yuǎn)場(chǎng)半徑取20~30倍弦長(zhǎng)太小會(huì)污染氣動(dòng)力系數(shù)太大浪費(fèi)網(wǎng)格量。近壁網(wǎng)格按y≈1來(lái)準(zhǔn)備。粗略估算第一層網(wǎng)格高度可以直接用平板邊界層公式uτ U·√(Cf/2)Cf ≈ 0.0576·Re_x^(-1/5)y? y·μ/(ρ·uτ)以Re3×10?、來(lái)流50m/s估算第一層高度大約在10??m量級(jí)。邊界層內(nèi)網(wǎng)格增長(zhǎng)率建議1.1~1.2尾緣變形區(qū)域x/c從0.7到1.0的弦向網(wǎng)格間距建議不大于0.01c否則變形后網(wǎng)格會(huì)被拉伸得很難看。如果你用的是Fluent Meshing而非ICEM有個(gè)操作細(xì)節(jié)容易踩新建的Group在網(wǎng)格顯示里不出現(xiàn)多半不是網(wǎng)格沒建上而是顯示對(duì)象的勾選沒對(duì)上。在Graphics面板里新建Scene把要顯示的Zone或Group加進(jìn)去并勾選刷新一下就會(huì)顯示。直接在Mesh Display的默認(rèn)設(shè)置里找新Group經(jīng)常找不到。網(wǎng)格生成后先檢查質(zhì)量skewness最好小于0.7minimum orthogonal quality大于0.2。動(dòng)網(wǎng)格案例的網(wǎng)格質(zhì)量余量要比定常算例留得更足因?yàn)樽冃芜^程會(huì)讓網(wǎng)格質(zhì)量持續(xù)下降。3. UDF逐行拆解俯仰加尾緣變形的核心邏輯3.1 為什么剛體旋轉(zhuǎn)要用角度增量而不是絕對(duì)角度這是整個(gè)UDF里最容易寫錯(cuò)的地方也最需要理解清楚。Fluent的DEFINE_GRID_MOTION在每一個(gè)時(shí)間步開始前被調(diào)用它執(zhí)行的是在當(dāng)前網(wǎng)格位置基礎(chǔ)上施加一個(gè)位移增量而不是直接把節(jié)點(diǎn)挪到某個(gè)絕對(duì)位置。如果你寫成“把節(jié)點(diǎn)坐標(biāo)設(shè)置為旋轉(zhuǎn)后的絕對(duì)坐標(biāo)”那么第一個(gè)時(shí)間步網(wǎng)格挪到位第二個(gè)時(shí)間步又基于被挪過的位置再設(shè)置一次位移會(huì)反復(fù)累積幾個(gè)步之后網(wǎng)格就廢了。正確做法是計(jì)算當(dāng)前步的角度增量dθ θ(tΔt) - θ(t)然后把這個(gè)增量對(duì)應(yīng)的節(jié)點(diǎn)位移加到當(dāng)前坐標(biāo)上。寫成代碼就是dx xr·(cos(dθ) - 1) - yr·sin(dθ)dy xr·sin(dθ) yr·(cos(dθ) - 1)其中(xr, yr)是節(jié)點(diǎn)相對(duì)旋轉(zhuǎn)中心的相對(duì)坐標(biāo)。這個(gè)方式無(wú)論dθ多大都能保持剛性旋轉(zhuǎn)的精確性不會(huì)產(chǎn)生小角度近似誤差。還有一個(gè)好處Fluent在一個(gè)物理時(shí)間步內(nèi)可能因?yàn)榫W(wǎng)格重構(gòu)問題多次調(diào)用動(dòng)網(wǎng)格函數(shù)但傳入的time和dtime是同一個(gè)時(shí)間步的這樣dθ每次算出來(lái)都是0不會(huì)重復(fù)疊加位移。用絕對(duì)角度寫法就會(huì)出大問題。3.2 DEFINE_GRID_MOTION完整代碼與注釋下面給出我實(shí)際使用的完整UDF代碼層面做了參數(shù)化處理方便調(diào)到自己的工況。#include udf.h #include dynamesh_tools.h #include math.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif /* 翼型參數(shù) */ #define CHORD 1.0 /* 弦長(zhǎng) */ #define PITCH_CX 0.25 /* 俯仰旋轉(zhuǎn)中心x坐標(biāo)取1/4弦點(diǎn) */ #define PITCH_CY 0.0 /* 旋轉(zhuǎn)中心y坐標(biāo) */ /* 俯仰運(yùn)動(dòng)參數(shù) */ #define PITCH_AMP 5.0 /* 俯仰振幅單位度 */ #define PITCH_FREQ 2.0 /* 俯仰頻率單位Hz */ /* 尾緣變形參數(shù) */ #define TAIL_X0 0.7 /* 尾緣變形起始位置x/c */ #define TAIL_AMP 0.03 /* 尾緣最大變形量單位m0.03倍弦長(zhǎng) */ #define TAIL_FREQ 4.0 /* 尾緣變形頻率單位Hz */ DEFINE_GRID_MOTION(airfoil_pitch_deform, domain, dt, time, dtime) { Thread *tf DT_THREAD(dt); face_t f; Node *v; int n; real theta, theta_pdt, dtheta; real cos_d, sin_d; real x, y, xr, yr, dxr, dyr; real dx, dy; real xr_norm, shape, dy_deform; /* 標(biāo)記當(dāng)前線程為網(wǎng)格變形線程這一步不能少 */ SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf)); /* 計(jì)算本時(shí)間步的角度增量 */ theta PITCH_AMP * M_PI / 180.0 * sin(2.0 * M_PI * PITCH_FREQ * time); theta_pdt PITCH_AMP * M_PI / 180.0 * sin(2.0 * M_PI * PITCH_FREQ * (time dtime)); dtheta theta_pdt - theta; cos_d cos(dtheta); sin_d sin(dtheta); begin_f_loop(f, tf) { f_node_loop(f, tf, n) { v F_NODE(f, tf, n); if (NODE_POS_NEED_UPDATE(v)) { NODE_POS_UPDATED(v); x NODE_X(v); y NODE_Y(v); /* 1. 剛體俯仰旋轉(zhuǎn)繞旋轉(zhuǎn)中心的位移增量 */ xr x - PITCH_CX; yr y - PITCH_CY; dxr xr * (cos_d - 1.0) - yr * sin_d; dyr xr * sin_d yr * (cos_d - 1.0); dx dxr; dy dyr; /* 2. 尾緣變形在尾緣局部疊加y向變形 */ if (x TAIL_X0) { xr_norm (x - TAIL_X0) / (CHORD - TAIL_X0); /* 二次形狀函數(shù)起始位置變形為0尾緣處變形最大 */ shape xr_norm * xr_norm; dy_deform TAIL_AMP * shape * sin(2.0 * M_PI * TAIL_FREQ * time); dy dy_deform; } NODE_X(v) dx; NODE_Y(v) dy; } } } end_f_loop(f, tf) }代碼本身并不長(zhǎng)但有幾個(gè)關(guān)鍵點(diǎn)需要重點(diǎn)解釋。首先SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf))是必須的它告訴Fluent當(dāng)前線程的網(wǎng)格節(jié)點(diǎn)需要更新位置。如果不設(shè)置UDF雖然會(huì)被調(diào)用但節(jié)點(diǎn)位置可能完全不變化問題是“函數(shù)執(zhí)行了網(wǎng)格紋絲不動(dòng)”很多新手在這個(gè)坑里耗很久。其次NODE_POS_NEED_UPDATE和NODE_POS_UPDATED是配套使用的防重復(fù)更新機(jī)制。在一個(gè)時(shí)間步里同一個(gè)節(jié)點(diǎn)可能被多個(gè)面共享如果不做這個(gè)判斷節(jié)點(diǎn)會(huì)在這個(gè)循環(huán)里被反復(fù)更新位移被疊加多次。這是寫好動(dòng)網(wǎng)格UDF的基本功。第三f_node_loop(f, tf, n)里的n是節(jié)點(diǎn)在當(dāng)前面上的局部編號(hào)每次循環(huán)拿到的v是一個(gè)指向節(jié)點(diǎn)的指針。Fluent允許一個(gè)節(jié)點(diǎn)被多個(gè)面共享但因?yàn)橛辛薔ODE_POS_NEED_UPDATE的判斷共享節(jié)點(diǎn)只被更新一次。編譯時(shí)選擇Compiled UDF不能用Interpreted模式因?yàn)榇a里用了dynamesh_tools.h。編譯成功后在Dynamic Mesh Zones的Deforming區(qū)域的Motion UDF下拉列表里選擇airfoil_pitch_deform。3.3 尾緣變形的形狀函數(shù)與變形范圍設(shè)定代碼里的shape xr_norm * xr_norm也就是從變形起始點(diǎn)x/c0.7到尾緣x/c1.0采用二次函數(shù)過渡。這樣保證了在起始位置變形量及其斜率都為0避免在x0.7處出現(xiàn)幾何突變——如果變形函數(shù)在起始點(diǎn)不光滑那個(gè)位置附近會(huì)產(chǎn)生很大的網(wǎng)格畸變很容易直接負(fù)體積。如果想要更光滑的過渡可以用Hermite型形狀函數(shù)shape xr_norm3 · (6·xr_norm2 - 15·xr_norm 10)這個(gè)函數(shù)在起始點(diǎn)和終點(diǎn)的一階導(dǎo)都是0變形輪廓更接近結(jié)構(gòu)模態(tài)里的懸臂梁一階彎曲振型。我實(shí)際對(duì)比過這個(gè)函數(shù)對(duì)網(wǎng)格質(zhì)量的保護(hù)明顯好于簡(jiǎn)單二次型代價(jià)只是多一行代碼。變形方向這里用了全局y方向是因?yàn)镹ACA0012上下表面本身關(guān)于x軸對(duì)稱尾緣垂直方向變形可以近似用y向表達(dá)。如果要做有彎度的翼型或變形方向沿局部表面法向需要更精細(xì)的處理在face循環(huán)里通過F_AREA(f, tf)取面的面積矢量歸一化得到法向再把變形位移沿法向施加。這會(huì)讓UDF更復(fù)雜但物理意義更準(zhǔn)確。還有一個(gè)經(jīng)驗(yàn)變形量和頻率不要一開始就拉滿。建議先用俯仰UDF單獨(dú)跑通再疊加尾緣變形。疊加時(shí)先給一半振幅確認(rèn)網(wǎng)格沒問題再逐步加大。4. 求解設(shè)置與收斂控制先穩(wěn)后動(dòng)是鐵律4.1 定常初場(chǎng)為什么要先“凍住”翼型算穩(wěn)定動(dòng)網(wǎng)格計(jì)算最忌諱的就是從均勻流場(chǎng)直接啟動(dòng)。如果初始化后就直接開瞬態(tài)、動(dòng)網(wǎng)格第一個(gè)時(shí)間步翼型開始轉(zhuǎn)動(dòng)尾緣開始變形流場(chǎng)會(huì)感受到一個(gè)劇烈的“沖擊”壁面附近必然產(chǎn)生非物理的壓力波輕則前幾個(gè)周期升力系數(shù)亂跳重則直接發(fā)散。我的標(biāo)準(zhǔn)流程是先用定常求解器在平均攻角位置把翼型固定住算一個(gè)穩(wěn)態(tài)流場(chǎng)待殘差降到1×10??以下升阻力系數(shù)不再明顯變化再切換為瞬態(tài)在瞬態(tài)計(jì)算開始的同時(shí)打開動(dòng)網(wǎng)格讓網(wǎng)格運(yùn)動(dòng)在一個(gè)相對(duì)真實(shí)的流場(chǎng)上逐步啟動(dòng)。如果定常計(jì)算本身就很難收斂先解決網(wǎng)格和湍流模型的問題不要指望動(dòng)網(wǎng)格能幫你“兜底”。動(dòng)網(wǎng)格只會(huì)放大初場(chǎng)的不穩(wěn)定不會(huì)修正它。實(shí)際操作中很多老工程師還會(huì)在切瞬態(tài)后先固定翼型再跑幾百步讓定常流場(chǎng)在瞬態(tài)格式下進(jìn)一步穩(wěn)定然后才真正啟動(dòng)動(dòng)網(wǎng)格。這個(gè)方法尤其適合雷諾數(shù)較高、邊界層敏感的算例。4.2 時(shí)間步長(zhǎng)、約化頻率與每步迭代次數(shù)本案例的參數(shù)我建議這樣取項(xiàng)目設(shè)定值說明來(lái)流速度50 m/s低速不可壓縮弦長(zhǎng)1 m參考長(zhǎng)度俯仰頻率2 Hz周期T0.5s約化頻率kπfc/U≈0.126較低約化頻率對(duì)時(shí)間步長(zhǎng)要求不算苛刻時(shí)間步長(zhǎng)0.0005~0.001 sT/500到T/1000每步內(nèi)迭代20~40次視殘差和升力系數(shù)穩(wěn)定情況網(wǎng)格最大位移小于最小網(wǎng)格尺寸的1/3重要的穩(wěn)定性判據(jù)約化頻率k是無(wú)量綱的振蕩頻率公式是kπfc/U。它直接決定了流動(dòng)非定常性的強(qiáng)弱。k0.126屬于低頻大幅振蕩流場(chǎng)能較快響應(yīng)翼型運(yùn)動(dòng)如果k接近0.5甚至更高時(shí)間步長(zhǎng)必須成倍縮小。時(shí)間步長(zhǎng)的選擇除了滿足每個(gè)周期的采樣點(diǎn)數(shù)還要考慮網(wǎng)格位移。一個(gè)時(shí)間步內(nèi)翼型表面節(jié)點(diǎn)移動(dòng)的距離不能超過當(dāng)?shù)刈钚【W(wǎng)格尺寸的三分之一否則Diffusion光順很容易產(chǎn)生負(fù)體積。用這個(gè)判據(jù)反推往往比單純按周期取步長(zhǎng)更有效。并行計(jì)算方面Fluent Launcher啟動(dòng)時(shí)在Parallel選項(xiàng)卡里設(shè)置的Processor進(jìn)程數(shù)一般不超過機(jī)器物理核心數(shù)。動(dòng)網(wǎng)格加UDF的計(jì)算我建議先跑串行或4核以下排查UDF和網(wǎng)格問題確認(rèn)穩(wěn)定后再上大規(guī)模并行否則日志文件里全是網(wǎng)格畸變報(bào)錯(cuò)排查效率極低。4.3 湍流模型與離散格式選型依據(jù)低速翼型繞流壓力基求解器是自然選擇。湍流模型我在這個(gè)案例里推薦兩段式策略先用Spalart-Allmaras模型把計(jì)算框架跑通得到初步結(jié)果后再用SST k-ω模型做正式計(jì)算。理由很直接SA模型只有單一湍流輸運(yùn)方程數(shù)值魯棒性好收斂難度低特別適合在動(dòng)網(wǎng)格調(diào)通階段使用但SA對(duì)強(qiáng)逆壓梯度下的流動(dòng)分離預(yù)測(cè)偏粗糙。翼型俯仰振蕩往往伴隨動(dòng)態(tài)失速尾緣附近的流動(dòng)會(huì)出現(xiàn)周期性分離與再附這時(shí)候SST k-ω對(duì)分離點(diǎn)和再附的捕捉準(zhǔn)確得多代價(jià)是收斂難度上升對(duì)網(wǎng)格質(zhì)量更敏感。離散格式的設(shè)置我有明確偏好壓力用Second Order動(dòng)量用Second Order Upwind湍流量也用Second Order Upwind。瞬態(tài)格式先用First Order Implicit起跑50~100步等流場(chǎng)結(jié)構(gòu)穩(wěn)定后切換到Bounded Second Order Implicit。壓力速度耦合推薦Coupled雖然每個(gè)迭代步成本高一些但整體時(shí)間步內(nèi)收斂更快動(dòng)網(wǎng)格工況下比SIMPLE族算法更穩(wěn)。我不建議在正式計(jì)算時(shí)開Solution Steering的自動(dòng)模式它為了魯棒性會(huì)主動(dòng)降低離散格式迎風(fēng)階數(shù)掩蓋網(wǎng)格運(yùn)動(dòng)帶來(lái)的真實(shí)數(shù)值行為結(jié)果就是“算完看著沒發(fā)散但曲線一塌糊涂”。手動(dòng)控制格式每步監(jiān)控殘差和升力系數(shù)才是最可靠的。5. 實(shí)戰(zhàn)踩坑記錄負(fù)體積、發(fā)散和曲線導(dǎo)出5.1 負(fù)體積網(wǎng)格到底在哪里被“擠爆”了跑動(dòng)網(wǎng)格的人對(duì)這條報(bào)錯(cuò)一定不陌生Negative Cell Volume。第一次遇到時(shí)基本是凌晨?jī)扇c(diǎn)盯著Console里的報(bào)錯(cuò)代碼一臉茫然。根據(jù)我的經(jīng)驗(yàn)負(fù)體積最常出現(xiàn)在兩個(gè)位置一是旋轉(zhuǎn)中心附近的邊界層網(wǎng)格二是尾緣變形起始點(diǎn)x/c0.7附近。前者的機(jī)理是網(wǎng)格在旋轉(zhuǎn)過程中被“壓扁”后者則是變形函數(shù)曲率突變導(dǎo)致的拉伸過度。排查套路要系統(tǒng)看Console里報(bào)出的單元ID在Display面板用Cell ID顯示方式把這幾個(gè)網(wǎng)格高亮出來(lái)先確認(rèn)“爆”在哪回看報(bào)錯(cuò)時(shí)間是物理時(shí)間還是某一步——是所有周期都會(huì)爆還是只在某個(gè)俯仰角附近爆如果是所有周期都爆優(yōu)先減小時(shí)間步長(zhǎng)如果只在大角度時(shí)刻爆說明光順?biāo)惴▉?lái)不及吸收位移需要加大Diffusion參數(shù)或調(diào)整變形形狀函數(shù)確認(rèn)負(fù)體積網(wǎng)格集中在極薄的邊界層區(qū)域時(shí)考慮減少每個(gè)時(shí)間步內(nèi)壁面節(jié)點(diǎn)位移而不是盲目加密網(wǎng)格——網(wǎng)格越密允許的位移反而越小。Debug時(shí)最好先關(guān)閉Remeshing把問題全部歸因到Smoothing身上。如果關(guān)閉重構(gòu)后網(wǎng)格能跑通再開啟Remeshing并仔細(xì)設(shè)置最小和最大尺寸標(biāo)尺。用尾緣變形做調(diào)試時(shí)我會(huì)故意把變形量設(shè)成0先跑一個(gè)32核、幾百步的純俯仰算例確認(rèn)光順參數(shù)沒有大問題再逐步加變形量。這樣能把變量隔離快速定位問題源頭。5.2 升力曲線毛刺與網(wǎng)格更新頻率的關(guān)系動(dòng)網(wǎng)格計(jì)算完成后把升力系數(shù)時(shí)程導(dǎo)出來(lái)經(jīng)常能看到曲線底部有密密麻麻的小鋸齒。試幾個(gè)案例之后你會(huì)明白這通常不是物理現(xiàn)象而是數(shù)值噪聲。毛刺的來(lái)源往往是每個(gè)時(shí)間步開始時(shí)的網(wǎng)格位置更新導(dǎo)致壁面附近的壓力場(chǎng)瞬間被擾動(dòng)。如果每步內(nèi)迭代次數(shù)太少、壓力場(chǎng)尚未充分恢復(fù)就進(jìn)入下一步鋸齒就會(huì)逐級(jí)累積。我的對(duì)策順序是增加每步內(nèi)迭代次數(shù)到30~50次觀察毛刺是否收斂如果仍然有鋸齒把時(shí)間步長(zhǎng)減半代價(jià)是計(jì)算量翻倍調(diào)大Diffusion參數(shù)讓網(wǎng)格位移在空間上更平滑減小局部網(wǎng)格體積突變率檢查Smoothing里的Spring Constant或Diffusion參數(shù)設(shè)置不要同時(shí)開多個(gè)高剛度設(shè)置。我自己的經(jīng)驗(yàn)總結(jié)是升力曲線小幅鋸齒可以接受但如果鋸齒幅值超過升力平均值的1%說明數(shù)值噪聲已經(jīng)大到會(huì)污染高階統(tǒng)計(jì)量必須處理后再繼續(xù)計(jì)算。5.3 Report Definition曲線導(dǎo)出與后處理經(jīng)驗(yàn)很多人在Fluent里定義了Report Definition計(jì)算完了卻不知道怎么把曲線數(shù)據(jù)導(dǎo)出來(lái)。操作其實(shí)很簡(jiǎn)單在樹形菜單的Results → Report Definitions里找到你定義的升力系數(shù)或阻力系數(shù)右鍵選擇Export會(huì)彈出一個(gè)保存CSV或文本文件的對(duì)話框里面就是每個(gè)時(shí)間步的物理時(shí)間與對(duì)應(yīng)量值。如果想要計(jì)算過程中自動(dòng)保存就在創(chuàng)建Report Definition時(shí)勾選Write to File并設(shè)置輸出文件路徑。計(jì)算結(jié)束后直接拿到完整時(shí)程數(shù)據(jù)不用再手動(dòng)點(diǎn)Export。如果要用Origin或matplotlib畫圖我習(xí)慣把第一個(gè)參數(shù)列設(shè)為t/T無(wú)量綱周期數(shù)第二列設(shè)為Cl第三列為Cd第四列為Cm。對(duì)比文獻(xiàn)數(shù)據(jù)時(shí)注意力矩參考點(diǎn)很多文獻(xiàn)的Cm參考點(diǎn)是1/4弦線Fluent里默認(rèn)參考點(diǎn)要自己核對(duì)。需要導(dǎo)出壁面壓力分布時(shí)用File → Export → Solution Data在Location里選翼型壁面Variable里勾選Pressure和Wall Shear可以導(dǎo)出各個(gè)工位的壓力系數(shù)分布。后續(xù)對(duì)所有時(shí)刻的壓力分布做積分還能還原出升力系數(shù)時(shí)程和Report Definition的直出結(jié)果互相校驗(yàn)。最后再多說一句關(guān)于這個(gè)算例我前后調(diào)了很多個(gè)版本最大的體會(huì)是動(dòng)網(wǎng)格發(fā)散十有八九不是求解器不行而是運(yùn)動(dòng)寫得不干凈要么是剛體旋轉(zhuǎn)用了絕對(duì)位移導(dǎo)致累計(jì)誤差要么是變形函數(shù)在節(jié)點(diǎn)處不光滑要么是Deforming區(qū)域和Smoothing參數(shù)配合不合適。尤其是“剛體加柔性變形”這種疊加運(yùn)動(dòng)坐標(biāo)系、增量寫法、起始范圍這三點(diǎn)想清楚整個(gè)案例就成功了一半。如果后面你遇到類似問題可以先用一個(gè)只有俯仰的UDF把網(wǎng)格魯棒性測(cè)利索再漸進(jìn)式地把尾緣變形加進(jìn)來(lái)。動(dòng)網(wǎng)格這個(gè)東西慢就是快急不來(lái)。 Fluent動(dòng)網(wǎng)格實(shí)現(xiàn)翼型俯仰+尾緣變形完整攻略
做風(fēng)力機(jī)葉片或者機(jī)翼的氣動(dòng)彈性分析時(shí)我經(jīng)常要面對(duì)一個(gè)不算特別復(fù)雜、但也非常容易翻車的需求翼型本身在繞某一點(diǎn)做俯仰振蕩與此同時(shí)尾緣還要疊加一定幅度的柔性變形。前者是典型的剛體運(yùn)動(dòng)對(duì)應(yīng)Fluent動(dòng)網(wǎng)格里的剛體區(qū)域加CG運(yùn)動(dòng)后者是典型的邊界變形對(duì)應(yīng)Deforming區(qū)域加網(wǎng)格光順。兩個(gè)拆開做教程一抓一大把都不難一旦要求它們同時(shí)作用在同一個(gè)翼型上很多人就開始犯嘀咕了——?jiǎng)泳W(wǎng)格區(qū)域到底怎么設(shè)UDF怎么寫同一個(gè)邊界能不能又剛體又變形網(wǎng)格會(huì)不會(huì)在某個(gè)時(shí)刻直接被擠爆……這篇文章就圍繞“Fluent動(dòng)網(wǎng)格實(shí)現(xiàn)翼型俯仰振蕩同時(shí)尾緣變形”這個(gè)組合把從方案選型、網(wǎng)格準(zhǔn)備、UDF編寫到求解器設(shè)置、常見問題排錯(cuò)的完整過程梳理一遍。案例本身選的是最簡(jiǎn)單的NACA0012翼型低速來(lái)流剛體俯仰加尾緣二次型變形但里面的方法論可以直接遷移到風(fēng)機(jī)葉片、渦輪葉片、直升機(jī)旋翼這類工程問題上。適合已經(jīng)會(huì)跑基本Fluent仿真、想進(jìn)階動(dòng)網(wǎng)格的工程師也適合正在做氣動(dòng)彈性課題、被“剛體柔性變形疊加”卡住的學(xué)生。1. 同樣都是動(dòng)網(wǎng)格為什么這個(gè)案例要單獨(dú)拿出來(lái)講1.1 先想清楚翼型在“動(dòng)”的到底是什么很多新手上來(lái)就寫UDF結(jié)果連自己要模擬的運(yùn)動(dòng)都沒拆清楚。翼型俯仰振蕩加尾緣變形聽起來(lái)是“一個(gè)運(yùn)動(dòng)”數(shù)學(xué)上其實(shí)是兩個(gè)獨(dú)立位移場(chǎng)的疊加。第一部分是剛體俯仰。整個(gè)翼型繞固定的彈性軸做正弦轉(zhuǎn)動(dòng)典型形式是θ(t) θ? A·sin(2πft)其中θ?是平均攻角A是俯仰振幅f是俯仰頻率。這個(gè)運(yùn)動(dòng)的特點(diǎn)是翼型表面上所有網(wǎng)格節(jié)點(diǎn)的相對(duì)位置不改變整體繞旋轉(zhuǎn)中心轉(zhuǎn)一個(gè)角度。在Fluent里最經(jīng)典的實(shí)現(xiàn)方式是DEFINE_CG_MOTION直接給定瞬時(shí)角速度即可。第二部分是尾緣變形。這個(gè)更微妙——它不是整體轉(zhuǎn)動(dòng)而是翼型尾緣附近的一段邊界相對(duì)它自身的初始位置做一個(gè)連續(xù)、光滑的偏移。比如尾緣點(diǎn)在某個(gè)時(shí)刻下移0.03倍弦長(zhǎng)然后離尾緣越近變形越大越往上游變形越小到某個(gè)起點(diǎn)位置變形恰好為零。這種運(yùn)動(dòng)不能用CG_MOTION表達(dá)因?yàn)樗茐牧恕皠傂浴薄獏^(qū)域內(nèi)節(jié)點(diǎn)之間的相對(duì)距離變了。所以本質(zhì)上是剛體位移場(chǎng) 柔性變形場(chǎng)按時(shí)間同步疊加。搞清楚這一點(diǎn)后面所有技術(shù)選擇都順理成章。1.2 剛體俯仰和尾緣變形的“疊加”難點(diǎn)難點(diǎn)恰恰出在Fluent動(dòng)網(wǎng)格的框架上。在一個(gè)動(dòng)網(wǎng)格計(jì)算里一個(gè)邊界區(qū)域通常只允許一種運(yùn)動(dòng)屬性設(shè)成Rigid Body整個(gè)區(qū)域做剛性平移/旋轉(zhuǎn)你無(wú)法單獨(dú)挑出尾緣那一小塊讓它再多動(dòng)一點(diǎn)設(shè)成Deforming所有節(jié)點(diǎn)的位置都由用戶自己控制剛體運(yùn)動(dòng)也得自己寫到UDF里。也就是說你不能把同一個(gè)壁面既設(shè)成剛體又設(shè)成Deforming。那怎么辦直覺的解決辦法是把翼型一分為二前段到尾緣上游設(shè)為Rigid Body做俯仰尾緣單獨(dú)設(shè)為Deforming做變形。中間切一個(gè)interface。這個(gè)思路我最早也試過結(jié)果是交界面兩側(cè)的網(wǎng)格密度必須高度匹配否則插值誤差會(huì)直接污染壁面壓力場(chǎng)而且邊界層在interface處被硬生生切斷動(dòng)網(wǎng)格光順在這種位置經(jīng)常出現(xiàn)負(fù)體積。真實(shí)算下來(lái)穩(wěn)定性很差。后來(lái)我換了一種思路既然Deforming模式允許我直接控制節(jié)點(diǎn)位置那為什么不在UDF里同時(shí)寫上“剛性旋轉(zhuǎn)位移 尾緣變形位移”讓Fluent的光順?biāo)惴ㄈヌ幚砭W(wǎng)格內(nèi)部更新這個(gè)思路最終穩(wěn)定跑通了也是下面這篇文章所有內(nèi)容的核心。整個(gè)方案的代價(jià)是UDF要自己寫但換來(lái)的是對(duì)流場(chǎng)物理的更可控、對(duì)運(yùn)動(dòng)疊加的完全接管。2. 方案選型放棄CG Motion把整個(gè)翼型交給GRID_MOTION2.1 Fluent動(dòng)網(wǎng)格三種方法哪些能用哪些不能用Fluent動(dòng)網(wǎng)格大體有三類手段Smoothing、Remeshing、Overset。選型前建議先把適用范圍圈清楚。方法適用場(chǎng)景本案例可行性Smoothing光順邊界位移中小幅度網(wǎng)格拓?fù)洳蛔兒诵氖侄闻浜螪iffusion光順Remeshing局部重構(gòu)大位移、大轉(zhuǎn)動(dòng)局部網(wǎng)格拓?fù)渲亟ㄝo助手段大俯仰角時(shí)需要Overset重疊網(wǎng)格多體大位移、相互穿越可用但殺雞用牛刀交界面插值額外耗成本本案例的運(yùn)動(dòng)量級(jí)如果控制在工程界很常見的范圍內(nèi)——俯仰振幅±5°、尾緣變形3%~5%弦長(zhǎng)——那么Smoothing加局部Remeshing足夠。Diffusion-based光順比彈簧光順的魯棒性好得多尤其適合邊界做旋轉(zhuǎn)運(yùn)動(dòng)的情況因?yàn)樗鼤?huì)把變形量“均勻擴(kuò)散”到全場(chǎng)而不是像彈簧一樣在局部積累。真正需要糾結(jié)的并不是這三種方法的取舍而是邊界區(qū)域運(yùn)動(dòng)屬性怎么設(shè)。我最終使用的是Deforming區(qū)域加DEFINE_GRID_MOTION讓一個(gè)UDF同時(shí)控制剛體旋轉(zhuǎn)和尾緣變形內(nèi)部網(wǎng)格交給Diffusion光順吸收變形。2.2 本案例的動(dòng)網(wǎng)格區(qū)域設(shè)置在Fluent里按照下面的方式設(shè)置動(dòng)網(wǎng)格區(qū)域思路很清晰打開Dynamic Mesh啟用Smoothing和Remeshing在Dynamic Mesh Zones中將翼型壁面設(shè)置為DeformingMotion UDF指定為后面要寫的airfoil_pitch_deform遠(yuǎn)場(chǎng)邊界設(shè)置為Stationary保持固定計(jì)算域內(nèi)部不需要額外指定動(dòng)態(tài)區(qū)域光順?biāo)惴〞?huì)自動(dòng)把壁面的位移擴(kuò)散到整個(gè)網(wǎng)格。這里有個(gè)容易誤解的點(diǎn)Deforming區(qū)域指定的雖然是翼型壁面但Fluent在調(diào)用UDF得到壁面節(jié)點(diǎn)位移后會(huì)通過Smoothing算法把位移逐步傳遞到內(nèi)部網(wǎng)格節(jié)點(diǎn)。所以UDF里只需要管壁面節(jié)點(diǎn)的位置內(nèi)部網(wǎng)格怎么動(dòng)是求解器的事。Smoothing參數(shù)里我習(xí)慣把Diffusion參數(shù)調(diào)到1.5~2.0。Diffusion參數(shù)越大壁面附近的位移衰減越快網(wǎng)格更傾向于在遠(yuǎn)場(chǎng)吸收變形這對(duì)保護(hù)邊界層質(zhì)量非常關(guān)鍵。如果案例的俯仰角較大再把Remeshing里Minimum Length Scale和Maximum Length Scale設(shè)成當(dāng)?shù)鼐W(wǎng)格尺寸的0.5倍和2倍目標(biāo)偏斜度0.7左右局部網(wǎng)格壞了就讓Fluent自動(dòng)重構(gòu)。2.3 計(jì)算域、網(wǎng)格與邊界層準(zhǔn)備要點(diǎn)網(wǎng)格是整個(gè)動(dòng)網(wǎng)格算例的地基。翼型幾何本身可以用NACA0012標(biāo)準(zhǔn)型值點(diǎn)生成弦長(zhǎng)取1m。計(jì)算域推薦C型或O型拓?fù)溥h(yuǎn)場(chǎng)半徑取20~30倍弦長(zhǎng)太小會(huì)污染氣動(dòng)力系數(shù)太大浪費(fèi)網(wǎng)格量。近壁網(wǎng)格按y≈1來(lái)準(zhǔn)備。粗略估算第一層網(wǎng)格高度可以直接用平板邊界層公式uτ U·√(Cf/2)Cf ≈ 0.0576·Re_x^(-1/5)y? y·μ/(ρ·uτ)以Re3×10?、來(lái)流50m/s估算第一層高度大約在10??m量級(jí)。邊界層內(nèi)網(wǎng)格增長(zhǎng)率建議1.1~1.2尾緣變形區(qū)域x/c從0.7到1.0的弦向網(wǎng)格間距建議不大于0.01c否則變形后網(wǎng)格會(huì)被拉伸得很難看。如果你用的是Fluent Meshing而非ICEM有個(gè)操作細(xì)節(jié)容易踩新建的Group在網(wǎng)格顯示里不出現(xiàn)多半不是網(wǎng)格沒建上而是顯示對(duì)象的勾選沒對(duì)上。在Graphics面板里新建Scene把要顯示的Zone或Group加進(jìn)去并勾選刷新一下就會(huì)顯示。直接在Mesh Display的默認(rèn)設(shè)置里找新Group經(jīng)常找不到。網(wǎng)格生成后先檢查質(zhì)量skewness最好小于0.7minimum orthogonal quality大于0.2。動(dòng)網(wǎng)格案例的網(wǎng)格質(zhì)量余量要比定常算例留得更足因?yàn)樽冃芜^程會(huì)讓網(wǎng)格質(zhì)量持續(xù)下降。3. UDF逐行拆解俯仰加尾緣變形的核心邏輯3.1 為什么剛體旋轉(zhuǎn)要用角度增量而不是絕對(duì)角度這是整個(gè)UDF里最容易寫錯(cuò)的地方也最需要理解清楚。Fluent的DEFINE_GRID_MOTION在每一個(gè)時(shí)間步開始前被調(diào)用它執(zhí)行的是在當(dāng)前網(wǎng)格位置基礎(chǔ)上施加一個(gè)位移增量而不是直接把節(jié)點(diǎn)挪到某個(gè)絕對(duì)位置。如果你寫成“把節(jié)點(diǎn)坐標(biāo)設(shè)置為旋轉(zhuǎn)后的絕對(duì)坐標(biāo)”那么第一個(gè)時(shí)間步網(wǎng)格挪到位第二個(gè)時(shí)間步又基于被挪過的位置再設(shè)置一次位移會(huì)反復(fù)累積幾個(gè)步之后網(wǎng)格就廢了。正確做法是計(jì)算當(dāng)前步的角度增量dθ θ(tΔt) - θ(t)然后把這個(gè)增量對(duì)應(yīng)的節(jié)點(diǎn)位移加到當(dāng)前坐標(biāo)上。寫成代碼就是dx xr·(cos(dθ) - 1) - yr·sin(dθ)dy xr·sin(dθ) yr·(cos(dθ) - 1)其中(xr, yr)是節(jié)點(diǎn)相對(duì)旋轉(zhuǎn)中心的相對(duì)坐標(biāo)。這個(gè)方式無(wú)論dθ多大都能保持剛性旋轉(zhuǎn)的精確性不會(huì)產(chǎn)生小角度近似誤差。還有一個(gè)好處Fluent在一個(gè)物理時(shí)間步內(nèi)可能因?yàn)榫W(wǎng)格重構(gòu)問題多次調(diào)用動(dòng)網(wǎng)格函數(shù)但傳入的time和dtime是同一個(gè)時(shí)間步的這樣dθ每次算出來(lái)都是0不會(huì)重復(fù)疊加位移。用絕對(duì)角度寫法就會(huì)出大問題。3.2 DEFINE_GRID_MOTION完整代碼與注釋下面給出我實(shí)際使用的完整UDF代碼層面做了參數(shù)化處理方便調(diào)到自己的工況。#include udf.h #include dynamesh_tools.h #include math.h #ifndef M_PI #define M_PI 3.14159265358979323846 #endif /* 翼型參數(shù) */ #define CHORD 1.0 /* 弦長(zhǎng) */ #define PITCH_CX 0.25 /* 俯仰旋轉(zhuǎn)中心x坐標(biāo)取1/4弦點(diǎn) */ #define PITCH_CY 0.0 /* 旋轉(zhuǎn)中心y坐標(biāo) */ /* 俯仰運(yùn)動(dòng)參數(shù) */ #define PITCH_AMP 5.0 /* 俯仰振幅單位度 */ #define PITCH_FREQ 2.0 /* 俯仰頻率單位Hz */ /* 尾緣變形參數(shù) */ #define TAIL_X0 0.7 /* 尾緣變形起始位置x/c */ #define TAIL_AMP 0.03 /* 尾緣最大變形量單位m0.03倍弦長(zhǎng) */ #define TAIL_FREQ 4.0 /* 尾緣變形頻率單位Hz */ DEFINE_GRID_MOTION(airfoil_pitch_deform, domain, dt, time, dtime) { Thread *tf DT_THREAD(dt); face_t f; Node *v; int n; real theta, theta_pdt, dtheta; real cos_d, sin_d; real x, y, xr, yr, dxr, dyr; real dx, dy; real xr_norm, shape, dy_deform; /* 標(biāo)記當(dāng)前線程為網(wǎng)格變形線程這一步不能少 */ SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf)); /* 計(jì)算本時(shí)間步的角度增量 */ theta PITCH_AMP * M_PI / 180.0 * sin(2.0 * M_PI * PITCH_FREQ * time); theta_pdt PITCH_AMP * M_PI / 180.0 * sin(2.0 * M_PI * PITCH_FREQ * (time dtime)); dtheta theta_pdt - theta; cos_d cos(dtheta); sin_d sin(dtheta); begin_f_loop(f, tf) { f_node_loop(f, tf, n) { v F_NODE(f, tf, n); if (NODE_POS_NEED_UPDATE(v)) { NODE_POS_UPDATED(v); x NODE_X(v); y NODE_Y(v); /* 1. 剛體俯仰旋轉(zhuǎn)繞旋轉(zhuǎn)中心的位移增量 */ xr x - PITCH_CX; yr y - PITCH_CY; dxr xr * (cos_d - 1.0) - yr * sin_d; dyr xr * sin_d yr * (cos_d - 1.0); dx dxr; dy dyr; /* 2. 尾緣變形在尾緣局部疊加y向變形 */ if (x TAIL_X0) { xr_norm (x - TAIL_X0) / (CHORD - TAIL_X0); /* 二次形狀函數(shù)起始位置變形為0尾緣處變形最大 */ shape xr_norm * xr_norm; dy_deform TAIL_AMP * shape * sin(2.0 * M_PI * TAIL_FREQ * time); dy dy_deform; } NODE_X(v) dx; NODE_Y(v) dy; } } } end_f_loop(f, tf) }代碼本身并不長(zhǎng)但有幾個(gè)關(guān)鍵點(diǎn)需要重點(diǎn)解釋。首先SET_DEFORMING_THREAD_FLAG(THREAD_T0(tf))是必須的它告訴Fluent當(dāng)前線程的網(wǎng)格節(jié)點(diǎn)需要更新位置。如果不設(shè)置UDF雖然會(huì)被調(diào)用但節(jié)點(diǎn)位置可能完全不變化問題是“函數(shù)執(zhí)行了網(wǎng)格紋絲不動(dòng)”很多新手在這個(gè)坑里耗很久。其次NODE_POS_NEED_UPDATE和NODE_POS_UPDATED是配套使用的防重復(fù)更新機(jī)制。在一個(gè)時(shí)間步里同一個(gè)節(jié)點(diǎn)可能被多個(gè)面共享如果不做這個(gè)判斷節(jié)點(diǎn)會(huì)在這個(gè)循環(huán)里被反復(fù)更新位移被疊加多次。這是寫好動(dòng)網(wǎng)格UDF的基本功。第三f_node_loop(f, tf, n)里的n是節(jié)點(diǎn)在當(dāng)前面上的局部編號(hào)每次循環(huán)拿到的v是一個(gè)指向節(jié)點(diǎn)的指針。Fluent允許一個(gè)節(jié)點(diǎn)被多個(gè)面共享但因?yàn)橛辛薔ODE_POS_NEED_UPDATE的判斷共享節(jié)點(diǎn)只被更新一次。編譯時(shí)選擇Compiled UDF不能用Interpreted模式因?yàn)榇a里用了dynamesh_tools.h。編譯成功后在Dynamic Mesh Zones的Deforming區(qū)域的Motion UDF下拉列表里選擇airfoil_pitch_deform。3.3 尾緣變形的形狀函數(shù)與變形范圍設(shè)定代碼里的shape xr_norm * xr_norm也就是從變形起始點(diǎn)x/c0.7到尾緣x/c1.0采用二次函數(shù)過渡。這樣保證了在起始位置變形量及其斜率都為0避免在x0.7處出現(xiàn)幾何突變——如果變形函數(shù)在起始點(diǎn)不光滑那個(gè)位置附近會(huì)產(chǎn)生很大的網(wǎng)格畸變很容易直接負(fù)體積。如果想要更光滑的過渡可以用Hermite型形狀函數(shù)shape xr_norm3 · (6·xr_norm2 - 15·xr_norm 10)這個(gè)函數(shù)在起始點(diǎn)和終點(diǎn)的一階導(dǎo)都是0變形輪廓更接近結(jié)構(gòu)模態(tài)里的懸臂梁一階彎曲振型。我實(shí)際對(duì)比過這個(gè)函數(shù)對(duì)網(wǎng)格質(zhì)量的保護(hù)明顯好于簡(jiǎn)單二次型代價(jià)只是多一行代碼。變形方向這里用了全局y方向是因?yàn)镹ACA0012上下表面本身關(guān)于x軸對(duì)稱尾緣垂直方向變形可以近似用y向表達(dá)。如果要做有彎度的翼型或變形方向沿局部表面法向需要更精細(xì)的處理在face循環(huán)里通過F_AREA(f, tf)取面的面積矢量歸一化得到法向再把變形位移沿法向施加。這會(huì)讓UDF更復(fù)雜但物理意義更準(zhǔn)確。還有一個(gè)經(jīng)驗(yàn)變形量和頻率不要一開始就拉滿。建議先用俯仰UDF單獨(dú)跑通再疊加尾緣變形。疊加時(shí)先給一半振幅確認(rèn)網(wǎng)格沒問題再逐步加大。4. 求解設(shè)置與收斂控制先穩(wěn)后動(dòng)是鐵律4.1 定常初場(chǎng)為什么要先“凍住”翼型算穩(wěn)定動(dòng)網(wǎng)格計(jì)算最忌諱的就是從均勻流場(chǎng)直接啟動(dòng)。如果初始化后就直接開瞬態(tài)、動(dòng)網(wǎng)格第一個(gè)時(shí)間步翼型開始轉(zhuǎn)動(dòng)尾緣開始變形流場(chǎng)會(huì)感受到一個(gè)劇烈的“沖擊”壁面附近必然產(chǎn)生非物理的壓力波輕則前幾個(gè)周期升力系數(shù)亂跳重則直接發(fā)散。我的標(biāo)準(zhǔn)流程是先用定常求解器在平均攻角位置把翼型固定住算一個(gè)穩(wěn)態(tài)流場(chǎng)待殘差降到1×10??以下升阻力系數(shù)不再明顯變化再切換為瞬態(tài)在瞬態(tài)計(jì)算開始的同時(shí)打開動(dòng)網(wǎng)格讓網(wǎng)格運(yùn)動(dòng)在一個(gè)相對(duì)真實(shí)的流場(chǎng)上逐步啟動(dòng)。如果定常計(jì)算本身就很難收斂先解決網(wǎng)格和湍流模型的問題不要指望動(dòng)網(wǎng)格能幫你“兜底”。動(dòng)網(wǎng)格只會(huì)放大初場(chǎng)的不穩(wěn)定不會(huì)修正它。實(shí)際操作中很多老工程師還會(huì)在切瞬態(tài)后先固定翼型再跑幾百步讓定常流場(chǎng)在瞬態(tài)格式下進(jìn)一步穩(wěn)定然后才真正啟動(dòng)動(dòng)網(wǎng)格。這個(gè)方法尤其適合雷諾數(shù)較高、邊界層敏感的算例。4.2 時(shí)間步長(zhǎng)、約化頻率與每步迭代次數(shù)本案例的參數(shù)我建議這樣取項(xiàng)目設(shè)定值說明來(lái)流速度50 m/s低速不可壓縮弦長(zhǎng)1 m參考長(zhǎng)度俯仰頻率2 Hz周期T0.5s約化頻率kπfc/U≈0.126較低約化頻率對(duì)時(shí)間步長(zhǎng)要求不算苛刻時(shí)間步長(zhǎng)0.0005~0.001 sT/500到T/1000每步內(nèi)迭代20~40次視殘差和升力系數(shù)穩(wěn)定情況網(wǎng)格最大位移小于最小網(wǎng)格尺寸的1/3重要的穩(wěn)定性判據(jù)約化頻率k是無(wú)量綱的振蕩頻率公式是kπfc/U。它直接決定了流動(dòng)非定常性的強(qiáng)弱。k0.126屬于低頻大幅振蕩流場(chǎng)能較快響應(yīng)翼型運(yùn)動(dòng)如果k接近0.5甚至更高時(shí)間步長(zhǎng)必須成倍縮小。時(shí)間步長(zhǎng)的選擇除了滿足每個(gè)周期的采樣點(diǎn)數(shù)還要考慮網(wǎng)格位移。一個(gè)時(shí)間步內(nèi)翼型表面節(jié)點(diǎn)移動(dòng)的距離不能超過當(dāng)?shù)刈钚【W(wǎng)格尺寸的三分之一否則Diffusion光順很容易產(chǎn)生負(fù)體積。用這個(gè)判據(jù)反推往往比單純按周期取步長(zhǎng)更有效。并行計(jì)算方面Fluent Launcher啟動(dòng)時(shí)在Parallel選項(xiàng)卡里設(shè)置的Processor進(jìn)程數(shù)一般不超過機(jī)器物理核心數(shù)。動(dòng)網(wǎng)格加UDF的計(jì)算我建議先跑串行或4核以下排查UDF和網(wǎng)格問題確認(rèn)穩(wěn)定后再上大規(guī)模并行否則日志文件里全是網(wǎng)格畸變報(bào)錯(cuò)排查效率極低。4.3 湍流模型與離散格式選型依據(jù)低速翼型繞流壓力基求解器是自然選擇。湍流模型我在這個(gè)案例里推薦兩段式策略先用Spalart-Allmaras模型把計(jì)算框架跑通得到初步結(jié)果后再用SST k-ω模型做正式計(jì)算。理由很直接SA模型只有單一湍流輸運(yùn)方程數(shù)值魯棒性好收斂難度低特別適合在動(dòng)網(wǎng)格調(diào)通階段使用但SA對(duì)強(qiáng)逆壓梯度下的流動(dòng)分離預(yù)測(cè)偏粗糙。翼型俯仰振蕩往往伴隨動(dòng)態(tài)失速尾緣附近的流動(dòng)會(huì)出現(xiàn)周期性分離與再附這時(shí)候SST k-ω對(duì)分離點(diǎn)和再附的捕捉準(zhǔn)確得多代價(jià)是收斂難度上升對(duì)網(wǎng)格質(zhì)量更敏感。離散格式的設(shè)置我有明確偏好壓力用Second Order動(dòng)量用Second Order Upwind湍流量也用Second Order Upwind。瞬態(tài)格式先用First Order Implicit起跑50~100步等流場(chǎng)結(jié)構(gòu)穩(wěn)定后切換到Bounded Second Order Implicit。壓力速度耦合推薦Coupled雖然每個(gè)迭代步成本高一些但整體時(shí)間步內(nèi)收斂更快動(dòng)網(wǎng)格工況下比SIMPLE族算法更穩(wěn)。我不建議在正式計(jì)算時(shí)開Solution Steering的自動(dòng)模式它為了魯棒性會(huì)主動(dòng)降低離散格式迎風(fēng)階數(shù)掩蓋網(wǎng)格運(yùn)動(dòng)帶來(lái)的真實(shí)數(shù)值行為結(jié)果就是“算完看著沒發(fā)散但曲線一塌糊涂”。手動(dòng)控制格式每步監(jiān)控殘差和升力系數(shù)才是最可靠的。5. 實(shí)戰(zhàn)踩坑記錄負(fù)體積、發(fā)散和曲線導(dǎo)出5.1 負(fù)體積網(wǎng)格到底在哪里被“擠爆”了跑動(dòng)網(wǎng)格的人對(duì)這條報(bào)錯(cuò)一定不陌生Negative Cell Volume。第一次遇到時(shí)基本是凌晨?jī)扇c(diǎn)盯著Console里的報(bào)錯(cuò)代碼一臉茫然。根據(jù)我的經(jīng)驗(yàn)負(fù)體積最常出現(xiàn)在兩個(gè)位置一是旋轉(zhuǎn)中心附近的邊界層網(wǎng)格二是尾緣變形起始點(diǎn)x/c0.7附近。前者的機(jī)理是網(wǎng)格在旋轉(zhuǎn)過程中被“壓扁”后者則是變形函數(shù)曲率突變導(dǎo)致的拉伸過度。排查套路要系統(tǒng)看Console里報(bào)出的單元ID在Display面板用Cell ID顯示方式把這幾個(gè)網(wǎng)格高亮出來(lái)先確認(rèn)“爆”在哪回看報(bào)錯(cuò)時(shí)間是物理時(shí)間還是某一步——是所有周期都會(huì)爆還是只在某個(gè)俯仰角附近爆如果是所有周期都爆優(yōu)先減小時(shí)間步長(zhǎng)如果只在大角度時(shí)刻爆說明光順?biāo)惴▉?lái)不及吸收位移需要加大Diffusion參數(shù)或調(diào)整變形形狀函數(shù)確認(rèn)負(fù)體積網(wǎng)格集中在極薄的邊界層區(qū)域時(shí)考慮減少每個(gè)時(shí)間步內(nèi)壁面節(jié)點(diǎn)位移而不是盲目加密網(wǎng)格——網(wǎng)格越密允許的位移反而越小。Debug時(shí)最好先關(guān)閉Remeshing把問題全部歸因到Smoothing身上。如果關(guān)閉重構(gòu)后網(wǎng)格能跑通再開啟Remeshing并仔細(xì)設(shè)置最小和最大尺寸標(biāo)尺。用尾緣變形做調(diào)試時(shí)我會(huì)故意把變形量設(shè)成0先跑一個(gè)32核、幾百步的純俯仰算例確認(rèn)光順參數(shù)沒有大問題再逐步加變形量。這樣能把變量隔離快速定位問題源頭。5.2 升力曲線毛刺與網(wǎng)格更新頻率的關(guān)系動(dòng)網(wǎng)格計(jì)算完成后把升力系數(shù)時(shí)程導(dǎo)出來(lái)經(jīng)常能看到曲線底部有密密麻麻的小鋸齒。試幾個(gè)案例之后你會(huì)明白這通常不是物理現(xiàn)象而是數(shù)值噪聲。毛刺的來(lái)源往往是每個(gè)時(shí)間步開始時(shí)的網(wǎng)格位置更新導(dǎo)致壁面附近的壓力場(chǎng)瞬間被擾動(dòng)。如果每步內(nèi)迭代次數(shù)太少、壓力場(chǎng)尚未充分恢復(fù)就進(jìn)入下一步鋸齒就會(huì)逐級(jí)累積。我的對(duì)策順序是增加每步內(nèi)迭代次數(shù)到30~50次觀察毛刺是否收斂如果仍然有鋸齒把時(shí)間步長(zhǎng)減半代價(jià)是計(jì)算量翻倍調(diào)大Diffusion參數(shù)讓網(wǎng)格位移在空間上更平滑減小局部網(wǎng)格體積突變率檢查Smoothing里的Spring Constant或Diffusion參數(shù)設(shè)置不要同時(shí)開多個(gè)高剛度設(shè)置。我自己的經(jīng)驗(yàn)總結(jié)是升力曲線小幅鋸齒可以接受但如果鋸齒幅值超過升力平均值的1%說明數(shù)值噪聲已經(jīng)大到會(huì)污染高階統(tǒng)計(jì)量必須處理后再繼續(xù)計(jì)算。5.3 Report Definition曲線導(dǎo)出與后處理經(jīng)驗(yàn)很多人在Fluent里定義了Report Definition計(jì)算完了卻不知道怎么把曲線數(shù)據(jù)導(dǎo)出來(lái)。操作其實(shí)很簡(jiǎn)單在樹形菜單的Results → Report Definitions里找到你定義的升力系數(shù)或阻力系數(shù)右鍵選擇Export會(huì)彈出一個(gè)保存CSV或文本文件的對(duì)話框里面就是每個(gè)時(shí)間步的物理時(shí)間與對(duì)應(yīng)量值。如果想要計(jì)算過程中自動(dòng)保存就在創(chuàng)建Report Definition時(shí)勾選Write to File并設(shè)置輸出文件路徑。計(jì)算結(jié)束后直接拿到完整時(shí)程數(shù)據(jù)不用再手動(dòng)點(diǎn)Export。如果要用Origin或matplotlib畫圖我習(xí)慣把第一個(gè)參數(shù)列設(shè)為t/T無(wú)量綱周期數(shù)第二列設(shè)為Cl第三列為Cd第四列為Cm。對(duì)比文獻(xiàn)數(shù)據(jù)時(shí)注意力矩參考點(diǎn)很多文獻(xiàn)的Cm參考點(diǎn)是1/4弦線Fluent里默認(rèn)參考點(diǎn)要自己核對(duì)。需要導(dǎo)出壁面壓力分布時(shí)用File → Export → Solution Data在Location里選翼型壁面Variable里勾選Pressure和Wall Shear可以導(dǎo)出各個(gè)工位的壓力系數(shù)分布。后續(xù)對(duì)所有時(shí)刻的壓力分布做積分還能還原出升力系數(shù)時(shí)程和Report Definition的直出結(jié)果互相校驗(yàn)。最后再多說一句關(guān)于這個(gè)算例我前后調(diào)了很多個(gè)版本最大的體會(huì)是動(dòng)網(wǎng)格發(fā)散十有八九不是求解器不行而是運(yùn)動(dòng)寫得不干凈要么是剛體旋轉(zhuǎn)用了絕對(duì)位移導(dǎo)致累計(jì)誤差要么是變形函數(shù)在節(jié)點(diǎn)處不光滑要么是Deforming區(qū)域和Smoothing參數(shù)配合不合適。尤其是“剛體加柔性變形”這種疊加運(yùn)動(dòng)坐標(biāo)系、增量寫法、起始范圍這三點(diǎn)想清楚整個(gè)案例就成功了一半。如果后面你遇到類似問題可以先用一個(gè)只有俯仰的UDF把網(wǎng)格魯棒性測(cè)利索再漸進(jìn)式地把尾緣變形加進(jìn)來(lái)。動(dòng)網(wǎng)格這個(gè)東西慢就是快急不來(lái)。