1、計算流體力學講義計算流體力學講義 第四講第四講 有限差分法（有限差分法（2）李新亮李新亮 ；力學所主樓；力學所主樓219； 82543801 知識點：知識點： 離散誤差的離散誤差的Fourier分析分析; 間斷周圍數值振蕩的原因；間斷周圍數值振蕩的原因； GVC格式；格式； 模型方程向模型方程向N-S方程的推廣；方程的推廣； 1講義、課件上傳至講義、課件上傳至 (流體中文網）流體中文網） - “流體論壇流體論壇” -“ CFD基礎理論基礎理論 ”Copyright by Li Xinliang 3.3 差分格式的進一步分析差分格式的進一步分析1. 耗散與色散誤差耗散與色散誤差2Copyrigh

2、t by Li Xinliang)sin()0 ,(0 xxuxutu精確解1階迎風xu0123456-1-0.500.51Exact1st upwind2nd upwindCopyright by Li Xinliangt=solutions of 1D convection Eq.2階迎風)2/()43(21xuuuujjjjxuuujjj/ )(1數值實驗數值實驗 boundaryperiodicxxxuxutu,2 , 0)sin()0 ,(0時間推進：時間推進： 3步步TVD型型Runge-Kutta, 且且時間步長足夠小（誤差忽略）時間步長足夠小（誤差忽略）空間離散：空間離散： 1

3、階及階及2階迎風格式階迎風格式 （20個網格點）個網格點）實驗觀察到的現象實驗觀察到的現象 兩類誤差：兩類誤差： 振幅誤差振幅誤差 相位誤差相位誤差 （波速誤差）（波速誤差） Copyright by Li Xinliang3對以上對以上“實驗現象實驗現象”進行理論分析進行理論分析半離散分析：半離散分析： 假設時間推進是精確的，僅假設時間推進是精確的，僅分析空間離散帶來的誤差（難度小、常用）分析空間離散帶來的誤差（難度小、常用）全離散分析：全離散分析： 同時分析時、空離散的誤差同時分析時、空離散的誤差 （難度大）（難度大）boundaryperiodicxexfcxfctfikx,2 , 0)

4、0 ,()0(0考查問題：實際上就是普通三角函數，實際上就是普通三角函數，采用復采用復數形式僅僅是為了理論推導數形式僅僅是為了理論推導方便。方便。用實數形式用實數形式 sin(kx), cos(kx)推導形推導形式上略顯繁瑣。式上略顯繁瑣。 精確解：ikxikctctxikeeetxf)(),(差分格式：0jxjfctf（1）xfffjjjx/ )(1)2/()43(21xuuufjjjjx 其他格式 jikxjexu)(假設對于：假設對于：有有jikxjxexku隱含假設隱含假設: 線性差分格式線性差分格式 非線性系統作用于單波，會產生非線性系統作用于單波，會產生多個諧波多個諧波 （2） 差

5、分沒有誤差差分沒有誤差xikkCopyright by Li Xinliang4令：jikxjetutxf)( ),(jikxjxxektutxf/)( ),(（1）式化為：0/)( )( jjikxikxxektucedttud0/)( )( xktucdttud“半離散化半離散化”： 空間導數差分計算，空間導數差分計算，時間方程（常微）精確計算時間方程（常微）精確計算xctkeutu/)0( )( xctkikxikxjjjeuetutxf/)0( )( ),(ixikk如果如果 ， 無誤差無誤差分析分析 （修正波數）與誤差的關系（修正波數）與誤差的關系kirikkk)/(/)0( ),(

6、xkctkxikxctkjijreeutxf 理想情況：理想情況： 的誤差導致解的幅值誤差的誤差導致解的幅值誤差 耗散誤差耗散誤差 的誤差導致解傳播速度的誤差的誤差導致解傳播速度的誤差 色散誤差色散誤差xkkkir , 0rkikjikxjexu)(假設對于：假設對于：有有jikxjxexku 反映了一個波內的點數。反映了一個波內的點數。PPW （波內的點數）（波內的點數）= Copyright by Li Xinliang5jikxjexu)(jikxjxexku耗散、色散誤差分別由修正波數耗散、色散誤差分別由修正波數的實部和虛部決定。的實部和虛部決定。k關鍵參數：關鍵參數： 修正波數修正波

7、數含義：含義： 反應波數（譜）空間內差分的誤差反應波數（譜）空間內差分的誤差jikxkkjefxf)(任意函數：jikxkkjefikxf)(定義：kkfikf求導數，精確解求導數，精確解差分解差分解jikxkkjxjef xk ifF/kkfxkFFourier 分析的任務分析的任務計算出計算出 ，并考差其與，并考差其與 的逼近程度。的逼近程度。kixik 考察格式考察格式分辨率（分辨率（resolution）的重要指標的重要指標 精度：精度： 反映反映 時的情況時的情況 分辨率：網格點數很少（例如波里面只有分辨率：網格點數很少（例如波里面只有6個點）時的性能個點）時的性能對于多尺度問題，分

8、辨率更重要。對于多尺度問題，分辨率更重要。 犧牲精度，提高分辨率犧牲精度，提高分辨率0 xxk 優秀的差分格式，優秀的差分格式，1個個波長里面波長里面6個點個點 即可即可/2精度精度 分辨率分辨率Copyright by Li Xinliang6如何計算修正波數？如何計算修正波數？jikxjexu)(jikxjxexku定義： 方法方法1. 理論計算理論計算 根據差分具體表達式及定義計算根據差分具體表達式及定義計算例例1：xuuujjjx1令jikxjeu 則：jjjjikxiikxxxikikxjjjxexkexeeexxuuu)1 ()(1)(1xk于是：sin)cos1 (1ieki1階

9、迎風階迎風sin,cos1irkk例例2：)2/()43(21xuuuujjjjx2階迎風階迎風)43(2243221iiikxjjjjxeexexuuuuj2/ )43(2iieek2/ )2sinsin4(,2/ )2coscos43(irkkCopyright by Li Xinliang7方法方法2： 數值計算數值計算jikxjexu)(jikxjxexku定義：Step 1）選取計算域）選取計算域0,2 , 計算網格（例如計算網格（例如64，128）Step 2）給定波數給定波數 k, 生成函數值生成函數值Step 3) 調用差分子程序，得到導數值調用差分子程序，得到導數值Step

10、4) 通過通過Fourier反變換，得到譜：反變換，得到譜：假設已有求差分的子程序（黑箱，已知是線性的）假設已有求差分的子程序（黑箱，已知是線性的）x線性黑箱ju強調：研究CFD本身，不能只使用理論手段，還要用數值手段jxjuv).2 , 1(Njeujikxj).2 , 1(NjvjNjikxjkjevNv11根據修正波數的定義，根據修正波數的定義， 有有xvkk Step 5) 改變改變k的值，重復的值，重復2-5， 得到得到 對于對于 的依賴關系。畫圖的依賴關系。畫圖kxk非線性情況會產生高次諧波，造成非線性情況會產生高次諧波，造成 step 4 中隱含的假設無法成立中隱含的假設無法成立

11、 將將Fourier分析手段拓展到非線性系統分析手段拓展到非線性系統 需要研究的課題需要研究的課題隱含條件：只有隱含條件：只有波數為波數為k的那個的那個譜不為譜不為0 （線性（線性系統）系統）Copyright by Li Xinliang8中心差分格式的色散特性中心差分格式的色散特性0： 精確解；精確解； 1 ： 4階普通階普通2： 6階普通；階普通； 3：4階緊致階緊致4： 6階緊致；階緊致； 5：6階超緊致階超緊致迎風差分格式的色散特性迎風差分格式的色散特性0： 精確解，精確解， 1： 2階迎風階迎風2： 5階迎風偏心階迎風偏心 3： 3階迎風緊致階迎風緊致4： 5階迎風緊致階迎風緊致

12、每個波長里面每個波長里面2個網格點，個網格點， 譜方法譜方法的分辨率，的分辨率， 差分法分辨率的極限（只差分法分辨率的極限（只有無窮階精度才能達到）有無窮階精度才能達到）2)int(WavePersPoPPW20階超緊階超緊致格式致格式 接近譜接近譜方法方法Copyright by Li Xinliang9不同差分格式的色散誤差曲線不同差分格式的色散誤差曲線結論：結論： 要求分辨率相同的情況下，要求分辨率相同的情況下， 采用高階格式可放寬空間網格步采用高階格式可放寬空間網格步長，從而減少計算量長，從而減少計算量重要方向：重要方向： 高分辨率差分格式高分辨率差分格式0： 精確解精確解1： 2階迎

13、風階迎風2: 3階迎風階迎風3： 3階迎風緊致階迎風緊致4： 5階迎風緊致階迎風緊致 指定誤差要求的情況下，不同差分格指定誤差要求的情況下，不同差分格式能模擬的最大式能模擬的最大 ( ( 越大，所需網越大，所需網格越少）格越少） /2PPW作業題作業題1： 構造高分辨率差分格式，并進行理論分析及數值實驗構造高分辨率差分格式，并進行理論分析及數值實驗 針對單波方程針對單波方程：0 xutu 對于空間導數，構造出一種不超過對于空間導數，構造出一種不超過6點格式；并進行點格式；并進行Fourier誤差分析，誤差分析，畫出畫出kr,ki的曲線。的曲線。 要求：精度不限；要求：精度不限； 網格基架點數不

14、超過網格基架點數不超過6個；個； 能夠分辨的波數范圍盡量寬；能夠分辨的波數范圍盡量寬； （即（即kr,ki曲線近可能接近準確解）曲線近可能接近準確解） 給出差分的具體表達式，給出差分的具體表達式， 畫出畫出kr,ki的曲線；的曲線； 說明構造格式的階數，并采用本說明構造格式的階數，并采用本PPT第第5頁的方法給出的精度驗證；頁的方法給出的精度驗證； 26154131231jjjjjjjjuauauauauauaxuu16514233241jjjjjjjjuauauauauauaxuu形如：另外，進行如下數值驗證：另外，進行如下數值驗證：)sin()0 ,(2 , 0, 0 xxuxxutu空間

15、采用空間采用20個網格點，采用新構造的差分格式離散；時間推進采用個網格點，采用新構造的差分格式離散；時間推進采用3步步Runge-Kutta方法，時間步長可足夠小（例如方法，時間步長可足夠小（例如0.01）。給出）。給出t=20,50兩個時兩個時刻的數值解，與精確解比較（畫圖），并給出數值解的刻的數值解，與精確解比較（畫圖），并給出數值解的L2模誤差。模誤差。10Copyright by Li Xinliang提示： 1. 如不使用優化技術，則格式構造方法簡單， Taylor展開后解代數方程組即可。 2. 建議嘗試使用優化技術 26154131231jjjjjjjjuauauauauauaxu

16、u例： 假設格式形式如下如果要求其有5階精度，則通過Taylor展開可得到6個方程，6個系數可直接解出。我們要求其有4階精度（當然3階，2階也可），于是Taylor展開只能提供5個方程。6個未知數（a1-a6）， 5個方程； 有1個自由參數。 調整這個自由參數，使得kr,ki曲線最為理想。 如何調整？ 1) 可以人工調整，觀察kr,ki曲線，選取滿意的。 2）可自動調整，設立一個優化目標函數。 例如 調整自由參數，使得該目標函數取最大值。思路：犧牲精度，提高分辨率 05. 0)(:*ik11Copyright by Li Xinliang附錄：附錄： 部分差分格式部分差分格式 j-2 j-1

17、j j+1 0 xuatu表中的迎風差分格式均針對表中的迎風差分格式均針對 a0當a0 時， 右側為“前”） 2） 根據根據GVC的思想構造格式的思想構造格式jjjjjjjjjjjjjjxuuuuwhenxuuuuuuwhenxuuuu11111121)2/()()2/()43(間斷前：快格式；間斷前：快格式；間斷后：慢格式；間斷后：慢格式；格式 GVC23） 改寫成為守恒型改寫成為守恒型非線性情況，通常守非線性情況，通常守恒型效果更好恒型效果更好jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF1111112/12/ )(2/ )3(xFFujjjx/ )(2/12/1NN

18、D格式 GVC2a17Copyright by Li Xinliang 4） a 0 時，時， 同樣思路構造同樣思路構造 （利用對稱性，僅需把下標（利用對稱性，僅需把下標j+k換成換成j-k即可）即可） 00axuatujjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF1111112/12/ )(2/ )3(采用GVC2a (NND 2a)格式的計算結果 消除振蕩xu00.8100.51NND 2aNND 2Exact solution1st upwindDeveloping of initial discontinuity2nd NND(C) Li Xin

19、liangt=0.2使用NND2a （守恒形式）；NND2 （普通形式）及1階迎風格式的計算結果1211121212/12/ )(2/ )3(jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuF將j換成j+118Copyright by Li Xinliang作業題作業題2： 構造更高分辨率的構造更高分辨率的GVC格式格式 對于空間導數，構造出一種不超過對于空間導數，構造出一種不超過6點的點的GVC格式。要求：格式。要求： a. 精度不限；精度不限； b. 網格基架點數不超過網格基架點數不超過6個；個； c. 求解模型方程求解模型方程 1,0axuatu5 . 015 . 00

20、)0 ,(xxxu計算結果間斷盡量保持計算結果間斷盡量保持“銳利銳利”;計算結果振蕩盡量小。計算結果振蕩盡量小。振蕩的定量判據：振蕩的定量判據： 總變差（總變差（Total Variation）：）： 間斷間斷“銳利銳利”的定量判據：的定量判據： 間斷區內的點數？間斷區內的點數？ （自行設計）（自行設計）111NjjjuuTV給出差分格式的表達式、色散給出差分格式的表達式、色散/耗散分析耗散分析 （ki,kr曲線）；曲線）；給出模型方程給出模型方程t=0.2的結果（空間的結果（空間100個網格點，計算域個網格點，計算域0,1，時間推進可采用，時間推進可采用3階階Runge-Kutta方法）；與

21、精確解及方法）；與精確解及NND2a進行比較（畫在同一張圖上）進行比較（畫在同一張圖上）xu00.8100.51NND 2aNND 2Exact solution1st upwindDeveloping of initial discontinuity2nd NND(C) Li Xinliangt=0.2建議： 利用優化方法19Copyright by Li Xinliang 3.5 從模型方程推廣到從模型方程推廣到N-S方程（方程（Euler方程）方程）0 xuatu?0?,0aa格式F+格式F-visf(U)UxtTTpEupuufEuU)(,(,),(2xf(U)xua

22、（教科書第6章）對流項：信息（波）從上游傳至下游對流項：信息（波）從上游傳至下游 上游信息更重要上游信息更重要 迎風差分迎風差分擴散項：擴散項： 信息從中心向周圍擴散信息從中心向周圍擴散 不區分上、下游不區分上、下游 中心差分中心差分迎風差分優點：迎風差分優點： 有效利用信息傳播的方向，增強穩定性有效利用信息傳播的方向，增強穩定性微分與差分方程的影響域N-S方程：單波方程：單波方程單波方程 一個波，容易判斷波傳播方向一個波，容易判斷波傳播方向N-S對流項（對流項（Euler） 方程組：多波問題，方程組：多波問題， 復雜復雜雙曲方程組的原則雙曲方程組的原則 特征分解，找到獨立傳播的波特征分解，找

23、到獨立傳播的波0000 xvtvxtxtxtkkkVVUSSUUAU1常系數矩陣常系數矩陣A的情況的情況 完全解耦，獨立求解完全解耦，獨立求解變系數矩陣變系數矩陣A的情況的情況 局部討論局部討論20Copyright by Li Xinliang1. Jacobian 系數矩陣及其性質系數矩陣及其性質xxUAf(U)Uf(u)Af(U)Uf)(f(U)UUU)f()()(1AUUUf(U)f(U)(xxxUAAUf(U)(重要性質uucucuuu222322311221)3(23010AtttUAAUf(U)(特點：特點： A 可以像常數一樣，和求導運算交換可以像常數一樣，和求導運算交換vis

24、f(U)UxtTTpEupuufEuU)(,(,),(221Copyright by Li Xinliang2. 對流項的分裂對流項的分裂 目的：目的： 確定波傳播方向，便于使用迎風差分確定波傳播方向，便于使用迎風差分方法：方法： 1） 逐點分裂逐點分裂 2） 嚴格特征分裂嚴格特征分裂 1） 逐點分裂逐點分裂 利用性質)(AUUAf(U)xxxSSA1kdiagkdiag2,2kk=+優點：耗散小缺點：導數間斷SSA1SUSUAf1方式 A:2)(2/122kkk) 1(2)(3(232cw特點：特點： 不必進行矩陣運算，計算量小不必進行矩陣運算，計算量小Steger-Warming 分裂分裂

25、wcucuucucuu232221321321)(2)(2) 1()()() 1(2) 1(22)(f)(fffffA: Steger-Warming 分裂分裂22Copyright by Li Xinliang Steger-Warming 具體步驟具體步驟 （以一維為例）（以一維為例）已知已知TEuU),(1） 計算計算2） 計算計算3） 計算計算4） 帶入（帶入（1）式得到）式得到5） 利用不同的迎風格式，分別計算利用不同的迎風格式，分別計算 pu,cucuu321,2)(2/122kkk)3 , 2 , 1(,kkkwcucuucucuu232221321321)(2)(2) 1()(

26、)() 1(2) 1(22)(f) 1(2)(3(232cw)(),(ffff（1）ff ,fffxxff,xxff0a0 xuatu0a（后差，前差）（后差，前差） 6） xxxfff計算計算7） 時間推進時間推進0 xtfU23Copyright by Li Xinliang二維問題的二維問題的steger-Warming 分裂分裂021yxtffUTpEuuvpuu)(,21fTpEvpvuvv)(,22fTEvu,U21fff令：則：WvuvuVvvvuuu2222)(2222421213202413024130430f) 1(2)(3(243cW10) 1(211kcuu12kcuu

27、/1k/2k21kcvv22kcvv22具體使用步驟，具體使用步驟， 以計算以計算 為例為例x1f令令 計算特征值計算特征值 分裂特征值，計算分裂特征值，計算 帶入左式，計算正、負流通矢帶入左式，計算正、負流通矢量量 計算計算0, 1cucuu4321,2)(2/122kkk)4 , 3 , 2 , 1(,kkk)(1ffxxx111fffy2f計算計算 設置設置 ，并注意，并注意1, 0cvcvv4321,對于曲線坐標系對于曲線坐標系021ffUt21111fffyxJJ僅需令yxJJ11, 三維問題同樣處理三維問題同樣處理 二維、三維具體二維、三維具體 公式見傅德薰等公式見傅德薰等計算空氣

28、動力學計算空氣動力學 4.7節節 （158-162）書中公式有一定的排版錯誤，使用前務必書中公式有一定的排版錯誤，使用前務必重新仔細重新仔細推導！推導！24Copyright by Li XinliangB: L-F分裂分裂AUf(U) 2/ )(, 2/ )(*IAAIAA特點： A正特征值 負特征值A2/ )(*UffUA=+缺點：耗散偏大局部局部L-F分裂分裂，每個點上計算，每個點上計算 全局全局L-F分裂分裂，全局（一維）上計算，全局（一維）上計算 cu *足夠大cu *數學性質（光滑性）數學性質（光滑性）最好，但耗散偏大最好，但耗散偏大)(max*cux常數與迎風格式結合，等價于人工

29、粘性與迎風格式結合，等價于人工粘性UfUfUfUfUfUffffxxxxxxxxxxxxxx*0*)(2)(21)()(21)()(21例如，可取例如，可取0 xtf(U)U22*xxtUf(U)U25Copyright by Li Xinliang2,2kkkkkkkkk方式很多=+S-W:L-F:=+2,2*kkkk0, 0kkVan Leer:=+221112) 1() 1(22) 1(cufcufff2121Macfff0f1Ma1MaffcuMa/0f1Ma26Copyright by Li Xinliang 分裂后分裂后 失去了失去了A的性質（可以的性質（可以像常數一樣與求導交換）

30、像常數一樣與求導交換）逐點分裂：逐點分裂： 優點：優點： 無需矩陣運算，計算量小無需矩陣運算，計算量小 缺點：缺點： 分裂后改變了特征方向，分裂后改變了特征方向， 耗散大耗散大xxx(U)f(U)ff(U)SUSUAf1利用了性質)(AUUAf(U)xxx一般情況下：一般情況下：xxxUAUAf變系數，變系數， 不能與導數交換不能與導數交換AAxBxUf AB實質：實質： 沒有做到解耦；沒有做到解耦； 只是把原變量重新組合，組合后波的傳播方向的保證只是把原變量重新組合，組合后波的傳播方向的保證 f+ 向正向正向傳播，向傳播，f-向負向傳播向負向傳播 缺點：缺點： 由于未解耦，各變量的誤差會相互

31、傳遞由于未解耦，各變量的誤差會相互傳遞 27Copyright by Li Xinliang概念澄清：概念澄清： 流通流通矢量分裂本身不帶來耗散矢量分裂本身不帶來耗散， 但其會但其會影響到差分的耗散影響到差分的耗散;舉例：0)(xUftUfff分裂過程),(21),(21UffUfffffxxxUaftUxxx00耗散如果差分格式無耗散（例如都用中心差分），則通量分裂不帶來耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂差分格式耗散分裂后的流場越偏離原先流場，則分裂后的流場越偏離原先流場，則總體耗散越大總體耗散越大fff精確滿足，不引入誤差！如使用低精差分度格式，如使用低精差分

32、度格式， 則對分裂形式敏感則對分裂形式敏感 （推薦使用特征分裂）（推薦使用特征分裂）如使用高精度格式（低耗散），則對分裂形式不敏感如使用高精度格式（低耗散），則對分裂形式不敏感 （可使用逐點分裂）（可使用逐點分裂）28Copyright by Li Xinliang2. 嚴格特征分裂嚴格特征分裂 基架點上凍結系數基架點上凍結系數常系數常系數方程組：方程組：00 xtxtUAUf(U)U0000 xvtvxtxtxtkkkVVUSSUUAU1完全解耦變系數情況變系數情況 局部凍結系數局部凍結系數0jjjxtUAU j-2 j-1 j j+1 在基架點上系數 不變jjxUAjjxUAjA計算：在差

33、分基架點上在差分基架點上Aj 不變，不變， 可按常矩陣處理可按常矩陣處理jj1jSSAj)()(jjjjjjjjxxxxxxVVSVSUSSUSSUA1j1j1jj1j局部凍結系數分別采用后差和前差USVj優點： 嚴格保證（局部）特征方向，數值解質量好；缺點： 大量矩陣運算，計算量大。29Copyright by Li XinliangUSVjUSVj通常寫成守恒型差分通常寫成守恒型差分，計算計算x)(Ufxxjjj/ )()(2/12/1ffUf j-2 j-1 j j+1 在基架點上系數 不變jjxUAjA)(2/12/112/12/1jjjjSVVfUSV2/1j具體步驟：具體步驟： 假

34、設已知假設已知 U， 且針對模型方程（線性單波方程）且針對模型方程（線性單波方程） 已構造出差分格式已構造出差分格式xvvvxvvvjjjjjj/ )(;/ )(2/12/12/12/10 xvatv（1）1） 計算出計算出2/12/112/1,jjjSS教材教材130頁的公式頁的公式(6.1.11-6.1.13), 式中用到各變量在式中用到各變量在j+1/2的值（例如的值（例如 ） 可使用可使用j, j+1 點值的算術平均點值的算術平均 （如（如 ） 或或Roe平均平均 （教材（教材6.4節）；節）； 由由 計算；方法很多，例如前面介紹的計算；方法很多，例如前面介紹的 或或 2/1ju2/ )(12/1jjjuuu30Copyright by Li Xinliang2kkk2*kk 均可本人感覺會限制精度（例如本人感覺會限制精度（例如2階）；但數值實驗沒發現問題。階）；但數值實驗沒發現問題。張樹海張樹海 JCP 2009 對其進行了深入探討；歡迎對該問題進行研究對其進行了深入探討；歡迎對該問題進行研究。2/1j2/1j2kkk2*kk2） 在網格