中科院計算流體力學(xué)最新講義CFD第講求代數(shù)方程組及網(wǎng)格生成.ppt

計算流體力學(xué)講義 第九講 代數(shù)方程的求解及網(wǎng)格生成 李新亮 ；力學(xué)所主樓219； 82543801,知識點(diǎn)：,1,講義、課件上傳至 (流體中文網(wǎng)） - “流體論壇” -“ CFD基礎(chǔ)理論 ” 講課錄像及講義上傳至網(wǎng)盤 /browse.aspx/.Public,Copyright by Li Xinliang,代數(shù)方程組的求解 網(wǎng)格生成,Copyright by Li Xinliang,2,知識回顧： 有限體積法基本流程,無粘項常用方法 （流過AB邊的通量）： a. 利用周圍點(diǎn)的值，計算出(I+1/2,J) 點(diǎn)處的物理量； 直接利用“差分格式” b. 利用該處的物理量，計算出流過AB邊的流通量 迎風(fēng)型方法需利用“通量分裂技術(shù)”,FVS類：,FDS類：,利用Riemann解,Reimann解： Godunov, Roe, HLL, HLLC,利用坐標(biāo)變換，轉(zhuǎn)化為一維Riemann問題,i,j,i-1,j,i+1,j,i,j+1,i,j+1,知識回顧： 隱格式,代數(shù)方程組的求解,Copyright by Li Xinliang,4,微分方程（組）,代數(shù)方程組,數(shù)值解（離散解）,差分 有限體積,大部分計算量,9.1.1 Gauss 消去法,消元,為了計算穩(wěn)定，通常使用主元消去法 列主元消去法； 全主元消去法 計算量： 乘法： 加法：,9.1 代數(shù)方程組求解的直接法,優(yōu)點(diǎn)： 簡單精確，缺點(diǎn)：計算量大,上三角 矩陣,Copyright by Li Xinliang,5,9.1.2 LU分解法,0,0,Step 1,Step 2,Step k,對角線上不能有0， 計算之前先交換矩陣A的元素，將主值交換到對角線上,Copyright by Li Xinliang,6,回代過程,計算量： 分解O(n3/3), 回代 O(n2),優(yōu)點(diǎn)： 1） 重復(fù)求解 ， 時，僅需一次LU分解，計算量?。?2） LU分解不破壞帶狀稀疏矩陣的性質(zhì)，可大幅減小計算量。,L帶寬的帶狀矩陣：,LU 分解: O(nL) 回代： O(nL),Copyright by Li Xinliang,7,9.1.3 帶狀矩陣求解的追趕法,追趕法：等價于帶狀矩陣的LU分解,例： 三對角矩陣,一般項：,邊界項：,追趕法,令：,Step 1: Step 2: Step 3: Step 4:,計算量： 9n次 （乘法） A為固定值時， 3n次（乘法）,簡單易用，計算量小,Copyright by Li Xinliang,8,9.2 代數(shù)方程組求解的迭代法,9.2.1 Jocabi 及Gauss-Seidel迭代,解出對角元素,Jocabi迭代,Gauss-Seidel迭代,“對角占優(yōu)”,Copyright by Li Xinliang,9,9.2.2 松弛迭代超松弛（SOR）、亞松弛,Step1: 采用Jocabian 或 Gauss-Seidel迭代產(chǎn)生新的值,Step 2: 進(jìn)行松弛,含義： 改變步長,超松弛,精確解,“步子邁大一些”，加快收斂,亞松弛,“步子邁小一些”，穩(wěn)定性好,收斂性： 對角占優(yōu)矩陣，Jocabian及Gauss-Seidel迭代可收斂,Copyright by Li Xinliang,10,舉例： Laplace方程的求解,五點(diǎn)格式,Jacobi迭代,Gauss-Seidel迭代,缺點(diǎn): 每迭代一步，信息只傳遞到周圍網(wǎng)格點(diǎn)，n很大時收斂較慢,n+1,n,n,n,n,n+1,n,n+1,n+1,Copyright by Li Xinliang,11,對稱Gauss-Seidel迭代 （SGS）,n+1,n,n+1,n+1,n+1,n,n,n,n+1,n+1,Step 1,Step 2,特點(diǎn)： 兩次掃描，反復(fù)迭代,Copyright by Li Xinliang,12,9.2.3 交替方向迭代（ADI）方法,例：,Step 1:認(rèn)為 已知 （使用上一步的值）， 求解三對角方程,得到中間步的值,Step 2: 代入中間步的值，求解三對角方程，得到n+1步的值,三對角方程采用追趕法求解，效率較高 在每一條線上采用直接法，信息快速傳遞，有利于收斂,Step 3: 重復(fù)以上兩個步驟，直至收斂,因追趕法實際上是LU分解法，因此又稱LU-ADI方法,n+1,n+1,n+1,n,n+1,n+1,n+1,n,n,n,Copyright by Li Xinliang,13,9.2.4 近似解法 LU-SGS方法,SGS 方法信息傳遞速度仍較慢，需要加速,近似LU分解,Step 1:,Step 2:,優(yōu)點(diǎn)： 不含任何迭代過程，兩步掃描即可完成，效率很高； 缺點(diǎn)： 近似LU分解，結(jié)果不夠準(zhǔn)確,OK， 不迭代,是LU分解的SGS方法，因此成為LU-SGS近似解法,Copyright by Li Xinliang,14,9.2.5 加速收斂的多重網(wǎng)格法,Gauss-Seidel迭代,含義： 線性系統(tǒng)，誤差 滿足同樣的方程,定義誤差：,1） 收斂速度的Fourier分析,增長（收斂）因子,含義：,極端高波數(shù)情況,迭代一次，誤差減小一半,極端低波數(shù)情況,收斂速度趨近于0,Copyright by Li Xinliang,15,策略： 多重網(wǎng)格 粗網(wǎng)格加速低波數(shù)擾動收斂，細(xì)網(wǎng)格加速高波數(shù)收斂,細(xì)網(wǎng)格,粗網(wǎng)格,使用多重網(wǎng)格法求解方程：,迭代方程：,以Jacobian迭代為例,修正方程,Step 1: 在細(xì)網(wǎng)格上迭代一定步數(shù)（無需收斂），得到中間步的值 Step 2: 將修正項 插值到粗網(wǎng)格上，并迭代求解 Step 3: 將求解后的修正項插值到細(xì)網(wǎng)格， 并計算出細(xì)網(wǎng)格上新的值 Step 4: 重復(fù)Step 1-3 直到收斂,修正項,Copyright by Li Xinliang,16,常用方法： V 型 及 W 型迭代,細(xì)網(wǎng)格,粗網(wǎng)格,更粗網(wǎng)格,細(xì)網(wǎng)格,粗網(wǎng)格,更粗網(wǎng)格,V 型迭代 W型迭代,Copyright by Li Xinliang,17,9.2.6 共軛梯度法,1. 求解對稱正定矩陣的共軛梯度法,化代數(shù)方程組問題為極值問題 （設(shè)A為對稱正定矩陣）,的最小值問題,可用最速下降法之類方法求解,例：,的分布,最小值點(diǎn)（0,1）,Copyright by Li Xinliang,18,求極值問題的最速下降法,思路： 沿當(dāng)?shù)靥荻确较蚯斑M(jìn),1）根據(jù)當(dāng)前位置 ， 計算當(dāng)?shù)靥荻确较颍?2） 沿該方向前進(jìn)，使得 達(dá)到極小值,方向,步長,“殘余向量”,特點(diǎn)： 沿當(dāng)?shù)靥荻确较蚯斑M(jìn)，直到不能前進(jìn)為止，然后以按照新的梯度方向前進(jìn)；相鄰路徑方向正交； 缺點(diǎn)：局部最速下降路線并非全局最速下降路線，因而收斂速度并非最優(yōu)。,“最速下降法”示意圖,Copyright by Li Xinliang,19,最速下降法的改進(jìn)：共軛梯度法,最速下降法：每步在一維空間求最優(yōu)解：,改進(jìn)：在二維空間尋求最優(yōu)解 （不再沿當(dāng)?shù)靥荻确较颍?舊路徑方向及當(dāng)前梯度方向所張成的二維空間,尋找該平面內(nèi)的極小值,解出,新線路：,修正,最終得：,特點(diǎn)： 相鄰兩次方向關(guān)于矩陣A正交 （稱為共軛）。,Copyright by Li Xinliang,20,2. 求解一般非奇異方程組的共軛梯度法,A為一般滿秩陣（非對稱正定陣）的情況,正則化方法：,對稱正定陣,Step 1: 設(shè)置初值,Step 2: 迭代推進(jìn),直到殘余向量足夠小為止,Copyright by Li Xinliang,21,9.3 代數(shù)網(wǎng)格生成法,基本思路： 通過代數(shù)方程計算出網(wǎng)格點(diǎn)的位置 優(yōu)點(diǎn)： 靈活、計算量小 缺點(diǎn)： 光滑性差， 過于依賴人工,如圖： 葉柵通道,已知計算域上邊界（紅線）及下邊界（藍(lán)線）的方程為： 和,則網(wǎng)格為：,其中 可控制法向的疏密分布,均勻分布；,在下壁面處密集分布,上下壁面兩側(cè)加密,Copyright by Li Xinliang,22,9.4 橢圓形方程網(wǎng)格生成法,A,B,C,D,E,F,A,B,C,D,E,F,對于如圖單聯(lián)通的計算域，可通過坐標(biāo)變換,變換到圖示矩形計算域,物理空間邊界 計算空間邊界 物理空間內(nèi)點(diǎn) 計算空間內(nèi)點(diǎn),物理空間,計算空間,通常： 給定邊界點(diǎn)的對應(yīng)關(guān)系 （代數(shù)方法） 通過求解方程獲得內(nèi)點(diǎn)的對應(yīng)關(guān)系,方程的邊值問題,橢圓型方程 邊值問題 拋物型方程 雙曲型方程,初邊值問題,橢圓型方程,Copyright by Li Xinliang,23,通常：,或,含義： 給物理空間的每個點(diǎn)找到計算空間的對應(yīng)位置。,注： 由于拓?fù)涞膶?yīng)性，物理空間必須是單聯(lián)通域,如果是多聯(lián)通的，可通過切割，形成單聯(lián)通域,Copyright by Li Xinliang,24,求解,1. 形式變換，改寫成以 為自變量,便于進(jìn)行差分求解,Copyright by Li Xinliang,25,離散化： 中心差分,離散方程,迭代求解： Jacobi, Seidel, SOR,LU-ADI, LU-SGS ,多重網(wǎng)格， ,Copyright by Li Xinliang,26,Step 1: 確定邊界網(wǎng)格,通常采用代數(shù)方法生成 一維網(wǎng)格，容易生成 注意： 1）邊界對應(yīng)關(guān)系容易出錯 2）考慮網(wǎng)格的疏密分布 （翼型尾緣區(qū)，激波區(qū)， 近壁區(qū)),Step 2: 利用上頁的離散方程， 解出全部網(wǎng)格坐標(biāo),不足： 內(nèi)部區(qū)域的網(wǎng)格分布不易控制 無法做到指定區(qū)域網(wǎng)格加密 無法保證網(wǎng)格正交,邊界網(wǎng)格可控，內(nèi)部網(wǎng)格只能“聽天由命”,方案1： 源項P,Q為0 ， 求解Laplace方程,Copyright by Li Xinliang,27,方案2： 設(shè)定源項P,Q 求解Poison方程,源項P,Q對網(wǎng)格的影響 數(shù)值實驗發(fā)現(xiàn)，在某點(diǎn)處加入點(diǎn)源P: P0 使 方向網(wǎng)格線發(fā)散,點(diǎn)源P,P0,P0,在某點(diǎn)處加入點(diǎn)源Q, 可對 方向網(wǎng)格線產(chǎn)生同樣效果,1） 網(wǎng)格線的匯聚,啟發(fā)： 在某條網(wǎng)格線上加入負(fù)的源項，可令網(wǎng)格匯聚,使網(wǎng)格匯聚于,Copyright by Li Xinliang,28,2） 邊界網(wǎng)格的正交（并指定邊界網(wǎng)格間距）,令：,源項在邊界處， 內(nèi)部衰減,利用邊界處網(wǎng)格的正交性及網(wǎng)格間距要求確定系數(shù)P和Q,指定值,基本思路： 以壁面線 處為例,網(wǎng)格線正交,指定法向網(wǎng)格間距,計算第2層網(wǎng)格線 上的坐標(biāo),通過差分計算邊界處的,如需要利用 的信息，可用上一迭代步的值,計算出邊界處的P,Q,根據(jù)指數(shù)衰減原則給出全場的P,Q,具體公式見： 傅德薰 計算空氣動力學(xué)284-286,Copyright by Li Xinliang,29,習(xí)題 9.1 網(wǎng)格生成,通過解橢圓型方程生成NACA0012翼型的網(wǎng)格 要求： 推薦采用圖示的C型網(wǎng)格， 網(wǎng)格點(diǎn)不限； 外邊界的位置不限； 可求解無源