一種固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法
【專利摘要】本發(fā)明屬于CFD網(wǎng)格生成技術(shù)領(lǐng)域，具體涉及一種可用于強(qiáng)迫振蕩、六自由度自由飛CFD計(jì)算的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法。該方法包括如下步驟：(1)生成初始網(wǎng)格，并將初始網(wǎng)格分為內(nèi)場(chǎng)與外場(chǎng)兩個(gè)區(qū)間；(2)求解初始網(wǎng)格的源項(xiàng)；(3)旋轉(zhuǎn)內(nèi)場(chǎng)網(wǎng)格，給定旋轉(zhuǎn)軸和轉(zhuǎn)角，剛性旋轉(zhuǎn)內(nèi)部網(wǎng)格，得到新狀態(tài)的內(nèi)部網(wǎng)格；(4)求解橢圓形方程光順外場(chǎng)網(wǎng)格。本發(fā)明所生成的網(wǎng)格保證了外邊界的固定不動(dòng)，同時(shí)保證了壁面處網(wǎng)格的正交性以及法向高度等特性，外場(chǎng)網(wǎng)格通過求解橢圓形方程光順質(zhì)量較好，在大角度轉(zhuǎn)動(dòng)后仍能保證較好的網(wǎng)格質(zhì)量。
【專利說明】
一種固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明屬于CFD網(wǎng)格生成技術(shù)領(lǐng)域，具體涉及一種可用于強(qiáng)迫振蕩、六自由度自 由飛CFD計(jì)算的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法。
【背景技術(shù)】
[0002] 近年來，隨著技術(shù)的進(jìn)步和飛行器性能需求的升級(jí)，基于非定常流動(dòng)的飛行器氣 動(dòng)特性研究越來越成為關(guān)注的熱點(diǎn)，而強(qiáng)迫運(yùn)動(dòng)和自由飛是研究飛行器非定常氣動(dòng)特性的 重要途徑。隨著計(jì)算機(jī)技術(shù)和數(shù)值計(jì)算方法的快速發(fā)展，CFD已經(jīng)成為解決上述非定常問 題的有效途徑，采用CFD計(jì)算非定常動(dòng)態(tài)特性必須首先解決網(wǎng)格運(yùn)動(dòng)的問題。
[0003] 現(xiàn)有結(jié)構(gòu)動(dòng)網(wǎng)格生成方法主要有兩種：剛性旋轉(zhuǎn)法和超限插值法。前者網(wǎng)格與物 面剛性固連，隨物體一起運(yùn)動(dòng)，網(wǎng)格形狀不發(fā)生改變，該方法簡(jiǎn)單，適用于任意外形、任意運(yùn) 動(dòng)，網(wǎng)格質(zhì)量好，完全符合幾何守恒率，但由于網(wǎng)格隨體運(yùn)動(dòng)，使遠(yuǎn)場(chǎng)邊界的位移和速度很 大，導(dǎo)致運(yùn)動(dòng)遠(yuǎn)場(chǎng)邊界處理不易；后者將插值法和剛性旋轉(zhuǎn)法加權(quán)混合的動(dòng)網(wǎng)格生成方法， 接近外邊界時(shí)，插值法網(wǎng)格的權(quán)重，靠近物面時(shí)，剛性旋轉(zhuǎn)法網(wǎng)格權(quán)重，該方法保持外邊界 的靜止，對(duì)復(fù)雜外形具有較好的適應(yīng)能力，但是對(duì)于大幅度的運(yùn)動(dòng)會(huì)出現(xiàn)網(wǎng)格扭曲及負(fù)體 積等現(xiàn)象，如大迎角俯仰運(yùn)動(dòng)。

【發(fā)明內(nèi)容】

[0004] 本發(fā)明的目的在于提供一種固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，對(duì)初始網(wǎng)格固定 外邊界，劃分內(nèi)、外場(chǎng)，同時(shí)通過對(duì)靠近壁面的內(nèi)場(chǎng)剛性旋轉(zhuǎn)、對(duì)外場(chǎng)求解橢圓型方程進(jìn)行 光順，獲得高質(zhì)量的結(jié)構(gòu)動(dòng)網(wǎng)格。
[0005] 為達(dá)到上述目的，本發(fā)明所采取的技術(shù)方案為：
[0006] 一種固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，包括如下步驟：
[0007] (1)生成初始網(wǎng)格，并將初始網(wǎng)格分為內(nèi)場(chǎng)與外場(chǎng)兩個(gè)區(qū)間，Sl為內(nèi)場(chǎng)，s2為外 場(chǎng)；
[0008] (2)求解初始網(wǎng)格的源項(xiàng)P、Q， _] ξχχ+ξ"= (χ ξ2+γξ2)Ρ(ξ, η)
[0010] ηχχ+ nyy= (χ n2+yn2)Q( ξ , η) (1)
[0011] 式中ξ、η為網(wǎng)格點(diǎn)在曲線坐標(biāo)系中的坐標(biāo)，x、y為網(wǎng)格點(diǎn)在物理坐標(biāo)系中的坐 標(biāo)；
[0012] (3)旋轉(zhuǎn)內(nèi)場(chǎng)網(wǎng)格，給定旋轉(zhuǎn)軸和轉(zhuǎn)角，剛性旋轉(zhuǎn)內(nèi)部網(wǎng)格，得到新狀態(tài)的內(nèi)部網(wǎng) 格；
[0013] (4)求解橢圓形方程光順外場(chǎng)網(wǎng)格。
[0014] 所述步驟（1)初始網(wǎng)格由Gridgen/Pointwise、ICEM工程軟件生成；初始網(wǎng)格應(yīng)保 證良好的正交性及光滑過渡特性，同時(shí)保證外場(chǎng)的大小滿足CFD計(jì)算要求，并可適當(dāng)增大 外場(chǎng)。
[0015] 所述步驟（2)初始網(wǎng)格源項(xiàng)求解過程如下：
[0016] 第一步：將式⑴由物理平面轉(zhuǎn)化到計(jì)算平面，即將自變量和因變量對(duì)調(diào)，得到式 (2)：
坐標(biāo)值；
[0019] 初始網(wǎng)格在邊界處滿足正交性條件，見式（3):[0020]
⑶
[0021] 堠
[0022] ⑷：
[0023] 分別用&和r ,點(diǎn)積上式，解出該邊界處源項(xiàng)的如下兩個(gè)表達(dá)式：[0024]
㈤
[0025] ⑹
[0026] 式中，rq，rqq，Γξ， ξ根據(jù)初始網(wǎng)格得坐標(biāo)值差分得到；
[0027] 第二步：得到邊界源項(xiàng)后，內(nèi)部源項(xiàng)采用以邊界源項(xiàng)為基礎(chǔ)的指數(shù)函數(shù)內(nèi)插方式 獲得，見式（7):
[0028]
[0029] (?)
[0030] 式中，ρ(ξ)和rU)分別是由式（5)決定的Ρ(ξ，η_)和Ρ(ξ，n_)，qU)和 sU)分別是由式（6)決定的QU，η_)和QU，η_) ;a，b，c，d都是正常數(shù)，表征由邊 界向流場(chǎng)內(nèi)部源項(xiàng)的衰減速率，較小的常數(shù)表示較慢的衰減。
[0031] 所述步驟（4)，外場(chǎng)網(wǎng)格的內(nèi)外邊界已經(jīng)給定，外場(chǎng)各網(wǎng)格點(diǎn)初值為初始網(wǎng)格，通 過迭代求解方程來獲得高質(zhì)量的外場(chǎng)網(wǎng)格，求解方法可以采用點(diǎn)松弛迭代方法、線松弛迭 代方法。
[0032] 所述點(diǎn)松弛迭代方法的步驟如下：
[0033] 求解方程式（1)，外邊界固定不動(dòng)，內(nèi)邊界通過旋轉(zhuǎn)求出，則η = η_，η_上邊 界條件給定，外場(chǎng)采用〇型網(wǎng)格，在ξ =c〇nst邊界上采用周期性邊界條件；設(shè)1_，】_分 別為ξ，η方向上的網(wǎng)格節(jié)點(diǎn)總數(shù)，兩個(gè)方向上的步長(zhǎng)均取為1，方程中涉及到的偏導(dǎo)數(shù)均 采用中心差分格式，即：
[0034]
[0035]
[0036]
[0037]
[0038]
[0039] 逐點(diǎn)超松弛方法步驟如下：
[0040] 第一步：先采用Gauss - Seidel改進(jìn)迭代，計(jì)算的循環(huán)安排是j = 2, 3,… ,Jmxl為外循環(huán)，i = 1，2, 3，··· I _為內(nèi)循環(huán)，則方程⑴在（i，j)點(diǎn)盡量采用最新值的 Gauss-Seidel改進(jìn)迭代的第n場(chǎng)的差分方程為：
[004 ^ 衝
[0042] 式中，gn，g22, g12也盡量采用最新值算得，則可求得2^ ,將之作為中間結(jié)果；: [00431
[0044] 第二步：引入超松弛迭代因子，將$與前一步解?丨，進(jìn)行加權(quán)平均，得到逐點(diǎn)超 松弛迭代公式：
[0046] w的取值范圍為[0, 2]。[0047] 本發(fā)明所取得的有益效果為：
[0045] CU)
[0048] 本發(fā)明將初始網(wǎng)格分為靠近物面的內(nèi)場(chǎng)以及帶外邊界的外場(chǎng)兩部分，由于壁面網(wǎng) 格對(duì)法向第一層高度、正交性等要求較高，因此對(duì)內(nèi)場(chǎng)采取剛性旋轉(zhuǎn)的方法生成動(dòng)網(wǎng)格，保 證了初始網(wǎng)格在壁面處的網(wǎng)格質(zhì)量，同時(shí)外邊界固定不動(dòng)，外場(chǎng)通過求解橢圓型方程進(jìn)行 光順的方式得到，其中方程中的源項(xiàng)采用初始網(wǎng)格的源項(xiàng)，在方程迭代求解中減少了一層 內(nèi)迭代，縮短了網(wǎng)格生成時(shí)間。所生成的網(wǎng)格保證了外邊界的固定不動(dòng)，同時(shí)保證了壁面處 網(wǎng)格的正交性以及法向高度等特性，外場(chǎng)網(wǎng)格通過求解橢圓形方程光順質(zhì)量較好，在大角 度轉(zhuǎn)動(dòng)后仍能保證較好的網(wǎng)格質(zhì)量。
【附圖說明】
[0049] 圖1為網(wǎng)格內(nèi)外場(chǎng)劃分示意圖；
[0050] 圖2為物理平面上的網(wǎng)格示意圖；
[0051] 圖3為計(jì)算平面上的網(wǎng)格示意圖；
[0052] 圖4為采用Gridgen軟件生成的初始網(wǎng)格示意圖；
[0053] 圖5為旋轉(zhuǎn)30度示意圖；
[0054] 圖6為旋轉(zhuǎn)180度示意圖；
[0055] 圖7為三維OREX飛船外形效果示意圖。
【具體實(shí)施方式】
[0056] 下面結(jié)合附圖和具體實(shí)施例對(duì)本發(fā)明進(jìn)行詳細(xì)說明。
[0057] 本發(fā)明所述固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法包括如下步驟：
[0058] 1、生成初始網(wǎng)格，并將初始網(wǎng)格分為內(nèi)場(chǎng)與外場(chǎng)兩個(gè)區(qū)間，見圖1，其中si為內(nèi) 場(chǎng)，s2為外場(chǎng)，ξ、η分別為二維曲線坐標(biāo)系的兩個(gè)坐標(biāo)符。由于網(wǎng)格生成軟件的普及， 初始網(wǎng)格一般由Gridgen/Pointwise、ICEM等工程軟件生成，以NACA0012翼型為例，采用 Gridgen軟件生成的初始網(wǎng)格見圖4。
[0059] 初始網(wǎng)格應(yīng)保證良好的正交性及光滑過渡特性，同時(shí)保證外場(chǎng)的大小滿足CFD計(jì) 算要求，并可適當(dāng)增大外場(chǎng)。
[0060] 2、求解初始網(wǎng)格的源項(xiàng)，以二維網(wǎng)格為例，源項(xiàng)即式（1)中的P、Q，為后續(xù)外場(chǎng)求 解橢圓型方耜摁他源頂佰"
[0061]
[0062]
[0063] 式中ξ、η為網(wǎng)格點(diǎn)在曲線坐標(biāo)系中的坐標(biāo)，x、y為網(wǎng)格點(diǎn)在物理坐標(biāo)系中的坐 標(biāo)，在求解橢圓型方程（式（1))光順網(wǎng)格的過程中，源項(xiàng)起著控制網(wǎng)格的正交性以及過渡 特性，且需再加入一層內(nèi)迭代來求解，由于初始網(wǎng)格具有較好質(zhì)量，為了縮短計(jì)算時(shí)間，后 續(xù)源項(xiàng)可直接采用初始網(wǎng)格的源項(xiàng)而不需內(nèi)迭代。初始網(wǎng)格源項(xiàng)求解過程如下：
[0064] 第一步，將式（1)由物理平面轉(zhuǎn)化到計(jì)算平面，即將自變量和因變量對(duì)調(diào)，得到式 (2)，物理平面與計(jì)算平面上得網(wǎng)格轉(zhuǎn)換見圖2、圖3。
[0065] ⑵
[0066] 式中：+方，其中?,)為單位正交向量，X，y為網(wǎng)格點(diǎn)在物理平面內(nèi)坐標(biāo)值
[0067]
[0068]
[0069]
[0070]
[0071] (3)
[0072]
[0073] (4)
[0074] 分別用&和r ,點(diǎn)積上式，可解出該邊界處源項(xiàng)的如下兩個(gè)表達(dá)式：[0075]
(5)
[0076] (6)
[0077] 式中，rq，rqq，Γξ， ξ可根據(jù)初始網(wǎng)格得坐標(biāo)值差分得到。
[0078] 第二步，得到邊界源項(xiàng)后，內(nèi)部源項(xiàng)采用以邊界源項(xiàng)為基礎(chǔ)的指數(shù)函數(shù)內(nèi)插方式
獲得，見式（7):
[0079]
[0080] (7)
[0081] 式中，ρ(ξ)和rU)分別是由式（5)決定的Ρ(ξ，η_)和Ρ(ξ，n_)，qU)和 sU)分別是由式（6)決定的QU，η_)和QU，η_)。a，b，c，d都是正常數(shù)，表征由邊 界向流場(chǎng)內(nèi)部源項(xiàng)的衰減速率，較小的常數(shù)表示較慢的衰減。
[0082] 類似的，可以求出另一邊界（ξ = const)上的源項(xiàng)分布。
[0083] 3、旋轉(zhuǎn)內(nèi)場(chǎng)網(wǎng)格。給定旋轉(zhuǎn)軸和轉(zhuǎn)角，剛性旋轉(zhuǎn)內(nèi)部網(wǎng)格，得到新狀態(tài)的內(nèi)部網(wǎng) 格，由于外場(chǎng)外邊界固定不動(dòng)，此時(shí)已經(jīng)固定了外場(chǎng)網(wǎng)格的內(nèi)外邊界（η = n_，n_)。
[0084] 4、求解橢圓形方程光順外場(chǎng)網(wǎng)格。外場(chǎng)網(wǎng)格的內(nèi)外邊界已經(jīng)給定，外場(chǎng)各網(wǎng)格點(diǎn) 初值為初始網(wǎng)格，通過迭代求解方程來獲得高質(zhì)量的外場(chǎng)網(wǎng)格，求解方法可以采用點(diǎn)松弛 迭代方法、線松弛迭代方法等，為了使迭代穩(wěn)定或加速收斂可以采用松弛法進(jìn)行，方程中各 項(xiàng)采用中心差分形式，以二維NACA0012翼型網(wǎng)格為例，采用逐點(diǎn)超松弛迭代方法的步驟推 導(dǎo)如下：
[0085] 求解方程式（1)，由于外邊界固定不動(dòng)，內(nèi)邊界通過旋轉(zhuǎn)求出，則η = η_，η_ 上邊界條件給定，由于外場(chǎng)采用0型網(wǎng)格，在ξ =Const邊界上采用周期性邊界條件，而不 必特別指定邊界條件。設(shè)I_，J_分別為ξ，η方向上的網(wǎng)格節(jié)點(diǎn)總數(shù)，兩個(gè)方向上的步長(zhǎng) 均取為1，方程中涉及到的偏導(dǎo)數(shù)均采用中心差分格式，即：
[0086]
[0087]
[0088]
[0089]
[0090] (8)
[0091] 逐點(diǎn)超松弛（SOR)方法步驟如下：
[0092] 第一步，先采用Gauss - Seidel改進(jìn)迭代，盡可能利用最新的值進(jìn)行計(jì)算，計(jì)算的 循環(huán)安排是j = 2, 3,…，J_ i為外循環(huán)，i = 1，2, 3,…I _為內(nèi)循環(huán)，則方程（1)在（i，j) 點(diǎn)盡量采用最新值的Gauss-Seidel改進(jìn)迭代的第n場(chǎng)的差分方程為：
[0093] M
[0094] 式中，gn，g22, g12也盡量采用最新值算得，則可求得2$，將之作為中間結(jié)果g : [00951
cm)
[0096] 第二步，引入超松弛迭代因子，將〗=與前一步解?="進(jìn)行加權(quán)平均，得到逐點(diǎn)超 松弛迭代公式：
[0097] (11)
[0098] w的取值范圍為[0, 2]。當(dāng)w = 1時(shí)就是Gauss-Seidel改進(jìn)迭代。
[0099] 三維動(dòng)網(wǎng)格生成方法可依此類推。
[0100] NACA0012翼型旋轉(zhuǎn)30度后生成的網(wǎng)格見圖5,旋轉(zhuǎn)180度后生成的網(wǎng)格見圖6。
[0101] 三維OREX飛船外形旋轉(zhuǎn)80度后生成的網(wǎng)格見圖7。
【主權(quán)項(xiàng)】
1. 一種固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，其特征在于：包括如下步驟： (1) 生成初始網(wǎng)格，并將初始網(wǎng)格分為內(nèi)場(chǎng)與外場(chǎng)兩個(gè)區(qū)間，si為內(nèi)場(chǎng)，s2為外場(chǎng)； (2) 求解初始網(wǎng)格的源項(xiàng)P、Q， ξ ΧΧ+ ξ yy= (x ξ^+Υξ^)Ρ( ξ , η) Πχχ+η"= (χ n"+Yn")Q( ξ , η) (1) 式中ξ、η為網(wǎng)格點(diǎn)在曲線坐標(biāo)系中的坐標(biāo)，x、y為網(wǎng)格點(diǎn)在物理坐標(biāo)系中的坐標(biāo)； (3) 旋轉(zhuǎn)內(nèi)場(chǎng)網(wǎng)格，給定旋轉(zhuǎn)軸和轉(zhuǎn)角，剛性旋轉(zhuǎn)內(nèi)部網(wǎng)格，得到新狀態(tài)的內(nèi)部網(wǎng)格； (4) 求解楠圓形方程光順外場(chǎng)網(wǎng)格。2. 根據(jù)權(quán)利要求1所述的固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，其特征在于：所述步驟 (1) 初始網(wǎng)格由Gridgen/化intwise、ICEM工程軟件生成；初始網(wǎng)格應(yīng)保證良好的正交性及 光滑過渡特性，同時(shí)保證外場(chǎng)的大小滿足C抑計(jì)算要求，并可適當(dāng)增大外場(chǎng)。3. 根據(jù)權(quán)利要求1所述的固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，其特征在于：所述步驟 (2) 初始網(wǎng)格源項(xiàng)求解過程如下： 第一步：將式（1)由物理平面轉(zhuǎn)化到計(jì)算平面，即將自變量和因變量對(duì)調(diào)，得到式 似：獄 式中；產(chǎn)=沿· +癡，其中?*，/為單位正交向量，X，y為網(wǎng)格點(diǎn)在物理平面內(nèi)坐標(biāo)值； 筑 &2=弓.巧/ 二為 + 分 &2=弓 初始網(wǎng)格在邊界處滿足正交性條件，見式（3);式中，r。，r。。，r &，r W根據(jù)初始網(wǎng)格得坐標(biāo)值差分得到； 第二步：得到邊界源項(xiàng)后，內(nèi)部源項(xiàng)采用W邊界源項(xiàng)為基礎(chǔ)的指數(shù)函數(shù)內(nèi)插方式獲得， 見式（7);式中，ρ(ξ)和Ηξ)分別是由式妨決定的Ρ(ξ，llmJ和Ρ(ζ，nmJ，qU)和3(ξ) 分別是由式（6)決定的9(ξ，llmJ和QU，nmJ ;a，b，c，d都是正常數(shù)，表征由邊界向流 場(chǎng)內(nèi)部源項(xiàng)的衰減速率，較小的常數(shù)表示較慢的衰減。4. 根據(jù)權(quán)利要求1所述的固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，其特征在于：所述步驟 (4)，外場(chǎng)網(wǎng)格的內(nèi)外邊界已經(jīng)給定，外場(chǎng)各網(wǎng)格點(diǎn)初值為初始網(wǎng)格，通過迭代求解方程來 獲得高質(zhì)量的外場(chǎng)網(wǎng)格，求解方法可W采用點(diǎn)松弛迭代方法、線松弛迭代方法。5. 根據(jù)權(quán)利要求4所述的固定外邊界的結(jié)構(gòu)動(dòng)網(wǎng)格生成方法，其特征在于：所述點(diǎn)松 弛迭代方法的步驟如下： 求解方程式（1)，外邊界固定不動(dòng)，內(nèi)邊界通過旋轉(zhuǎn)求出，則η = nmw nmm上邊界條 件給定，外場(chǎng)采用0型網(wǎng)格，在ξ = const邊界上采用周期性邊界條件；設(shè)分別為 ξ，η方向上的網(wǎng)格節(jié)點(diǎn)總數(shù)，兩個(gè)方向上的步長(zhǎng)均取為1，方程中涉及到的偏導(dǎo)數(shù)均采用 中必差分格式，即；逐點(diǎn)超松弛方法步驟如下： 第一步；先采用Gauss - Seidel改進(jìn)迭代，計(jì)算的循環(huán)安排是j = 2, 3,…，Jm。、1為外 循環(huán)，1 = 1，2,3，."1。。,為內(nèi)循環(huán)，則方程（1)在（1，如點(diǎn)盡量采用最新值的6曰1133-561(161 改進(jìn)迭代的第η場(chǎng)的差分方程為：式中，gii，拓2, g。也盡量采用最新值算得，則可求得說：/，將之作為中間結(jié)果;t.i:第二步；引入超松弛迭代因子，將；與前一步解為:7|，進(jìn)行加權(quán)平均，得到逐點(diǎn)超松弛 迭代公式：(11) W的取值范圍為[0, 2]。
【文檔編號(hào)】G06F17/50GK105843977SQ201510023282
【公開日】2016年8月10日
【申請(qǐng)日】2015年1月16日
【發(fā)明人】苗萌, 盧鳳翎, 周禹, 閔昌萬
【申請(qǐng)人】北京臨近空間飛行器系統(tǒng)工程研究所, 中國(guó)運(yùn)載火箭技術(shù)研究院