譜元方法課程

譜元方法課程第1頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@21、正交函數(shù)系與譜近似1.1正交函數(shù)系與正交多項(xiàng)式1.2函數(shù)的Fourier展開1.3Chebyshev譜逼近(離散函數(shù)的Fourier展開)1.4插值函數(shù)的導(dǎo)數(shù)1.5二維函數(shù)的Chebyshev譜近似第2頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@31.1正交函數(shù)系與正交多項(xiàng)式設(shè)函數(shù)系{j(x)}，如滿足條件，則稱函數(shù)系{j(x)}，在區(qū)間[a,b]上關(guān)于權(quán)函數(shù)w(x)正交，當(dāng)j(x)是代數(shù)多項(xiàng)式時，稱為正交多項(xiàng)式。如Chebyshev多項(xiàng)式，

第3頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@4Chebyshev多項(xiàng)式的正交性第4頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@51.2函數(shù)的Fourier展開函數(shù)f(x)在區(qū)間[a,b]上滿足Dirichlet條件，且加權(quán)平方可積，即對于權(quán)函數(shù)w(x)，有，則f(x)可在[a,b]上以w(x)為權(quán)函數(shù)，按正交多項(xiàng)式φn(x)展成廣義Fourier級數(shù)，設(shè){j(x)}為在節(jié)點(diǎn){xk,k=0,1,…,N}上的正交函數(shù)系，權(quán)為{wk>0,k=0,1,…,N}，設(shè)f(x)為在節(jié)點(diǎn)系{xk,k=0,1,…,N}上取值的已知函數(shù)，則函數(shù)f(x)的Fourier展開定義為，第5頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@61.3Chebyshev譜逼近(離散函數(shù)的Fourier展開)取節(jié)點(diǎn)系{xk,k=0,1,…,N}為N階Chebyshev多項(xiàng)式的極值點(diǎn)(稱Gauss-Lobatto點(diǎn))，即，權(quán)系數(shù)則上面的展開式變?yōu)椋耸揭嗫煽醋鱢(x)在配置點(diǎn)上的插值,g(x)為插值函數(shù)。第6頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@71.4插值函數(shù)的導(dǎo)數(shù)對微分得：插值函數(shù)在配置點(diǎn)的導(dǎo)數(shù)可表示為，第7頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@81.5二維函數(shù)的Chebyshev譜近似對二維函數(shù)u(ξ,η)，在標(biāo)準(zhǔn)正方形單元內(nèi)，定義ξ,η方向節(jié)點(diǎn)系{ξj，j=0,1,…,Nx}和{ηk，k=0,1,…,Ny}，第8頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@92譜方法求解微分方程2.1解微分方程的加權(quán)殘量法MWR

（MethodofWeightResiduals)2.2譜方法求解微分方程

2.3小結(jié)第9頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@102.1解微分方程的加權(quán)殘量法MWR（MethodofWeightResiduals)基本思想MWR的步驟MWR的基本方法第10頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@11加權(quán)殘量法基本思想(1)設(shè)微分方程邊值問題為，加權(quán)余量法，是求微分方程形式如下的近似解。a=[a0,a1,…,aN]T為待定向量，0,1,…,N為為基函數(shù)，且是一組線性無關(guān)的函數(shù)第11頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@12加權(quán)殘量法基本思想(2)“殘量”定義如下，顯然只有當(dāng)試函數(shù)為邊值問題的精確解時，余量(2.31)才為零。加權(quán)殘量法適當(dāng)選取待定向量a=[a0,a1,…,aN]T，使得“殘量”極小。通常某一權(quán)函數(shù)系{wi(x),i=0,1,…,N}，使得權(quán)函數(shù)與“殘量”正交，來確定待定向量。

第12頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@13MWR的步驟①選取一組滿足要求的基函數(shù)②構(gòu)造試函數(shù)并滿足所有邊界條件；③選取一組權(quán)函數(shù)④運(yùn)用MWR準(zhǔn)則，得到的代數(shù)方程組；⑤解上述代數(shù)方程組，確定待定參數(shù)。第13頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@14MWR的基本方法按權(quán)函數(shù)的不同有以下幾種基本方法Galerkin法(GalerkinMethod)配置法(CollocationMethod)Tan方法(TanMethod)

類似于Galerkin法，Tan方法中，試函數(shù)不需要滿足邊界條件，增加邊界約束系數(shù)來滿足邊界條件

第14頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@152.2譜方法求解微分方程

以非線性Burger方程為例來說明譜方法的離散過程,并區(qū)別Galerkin法,Tau方法和Collocation方法。Burger方程如下:

第15頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@16ChebyshevGalerkin方法

(GalerkinMethod方程的近似解表達(dá)為，

由MWR準(zhǔn)則得,

第16頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@17ChebyshevTau方法(TanMethod方程的近似解表達(dá)為，

由MWR準(zhǔn)則得,

邊界約束,第17頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@18ChebyshevCollocation方法(CollocationMethod)選取節(jié)點(diǎn)系為N階Chebyshev多項(xiàng)式的極值點(diǎn),

方程的近似解用下式表示，由Collocation方程得,第18頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@192.3小結(jié)一種求解偏微分方程的高階精度數(shù)值方法。屬于求解偏微分方程加權(quán)殘量法的特例。譜方法使用高階正交多項(xiàng)式作為展開函數(shù)。譜方法最受人青睞的優(yōu)越性在于它具有“無窮階”的收斂速度，其確切含義為，若原微分方程的解無限可微，則由適當(dāng)?shù)淖V方法所得到的近似解對原問題的收斂速度比1/N的任何冪都更快，這里N是所取基函數(shù)的個數(shù)。Kreiss和Oliger,OrszagGottlieb和Orszag第19頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@203譜元方法求解微分方程回顧譜方法———使用加權(quán)殘量法譜元方法基本思想常微分方程邊值問題的Galerkin變分原理偏微分方程的Galerkin變分原理Galerkin逼近解常微分方程元素矩陣的形成偏微分方程元素矩陣的形成極坐標(biāo)系下環(huán)形區(qū)域的譜元方法總剛度矩陣的形成第20頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@213.1回顧譜方法———使用加權(quán)殘量法

選取一組滿足要求的基函數(shù)

構(gòu)造試函數(shù)，滿足所有的邊界條件

選取一組權(quán)函數(shù)

利用MWR準(zhǔn)則，得到ai的代數(shù)方程

求解上述代數(shù)方程組，確定待定參數(shù)第21頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@223.2譜元方法基本思想使用有限元方法的思想，將求解域分成若干子域

采用譜方法中譜近似技術(shù)代替有限元中的插值函數(shù)

采用有限元中等參單元,亦可逼近復(fù)雜求解邊界

基函數(shù)的選用在一個單元上選取在整個區(qū)域上選取計算精度不是太高比較高單元的選擇三角形，四邊形等四邊形提高精度方法網(wǎng)格加密不改變插值次數(shù)增加基函數(shù)多項(xiàng)式的次數(shù)第22頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@233.3常微分方程邊值問題的Galerkin變分原理考察二點(diǎn)邊值問題上述問題的Galerkin變分問題是：求，使得，其中第23頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@243.4偏微分方程的Galerkin變分原理考察二維橢圓邊值問題上述問題的Galerkin變分問題是：求，使得，第24頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@253.5Galerkin逼近解設(shè)Vh是V的有限維子空間，當(dāng)h->0時，Vh的維數(shù)無限增加，直到充滿V為止，那么Galerkin逼近解，使得，設(shè)是Vh的基函數(shù)系，那么由于，故有，其中得由于的任意性，得，即得線性代數(shù)方程組，可以證明此代數(shù)方程組的解唯一存在，即可解出Galerkin逼近解。第25頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@263.6常微分方程元素矩陣的形成

考察以下微分方程用譜方法求解此微分方程,先將求解區(qū)域分成若干單元,設(shè)單元總數(shù)為N,考察第i個單元ei,將ei轉(zhuǎn)化為標(biāo)準(zhǔn)求解單元,設(shè)x為總體坐標(biāo),為標(biāo)準(zhǔn)求解單元內(nèi)的局部坐標(biāo),并設(shè)

第26頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@27插值點(diǎn)和插值函數(shù)在第i個元素內(nèi)選取Nx+1個點(diǎn)作為u(x)和f(x)的插值點(diǎn),本文選取階Chebyshev多項(xiàng)式的極值點(diǎn)作為插值點(diǎn),即

,u(x)和f(x)可以表達(dá)為，第27頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@28元素矩陣形成對一維邊值問題,在ei單元內(nèi)線性方程式變?yōu)?/p>

其中，上面線性方程組變?yōu)椋?8頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@29的推導(dǎo)過程第29頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@303.7偏微分方程元素矩陣的形成具體考察Helmholtz方程，其中，作以下坐標(biāo)變換在標(biāo)準(zhǔn)正方形單元內(nèi)，使用偽譜Chebyshev譜逼近u(ξ,η)，則元素內(nèi)線性方程組式變?yōu)?，?0頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@31總體坐標(biāo)系與局部坐標(biāo)系第31頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@32元素矩陣第32頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@333.8極坐標(biāo)系下環(huán)形區(qū)域的譜元方法在極坐標(biāo)系下，Helmholtz方程表達(dá)為，將極坐標(biāo)系下環(huán)形段單元，變換為直角坐標(biāo)系下的正方形單元，坐標(biāo)變換關(guān)系如下，元素矩陣的形式為，第33頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@34極坐標(biāo)坐標(biāo)系與局部坐標(biāo)系第34頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@35元素矩陣式中各變量的表達(dá)式如下，第35頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@363.9總剛度矩陣的形成節(jié)點(diǎn)影響元素集，節(jié)點(diǎn)所在的元素的并集影響節(jié)點(diǎn)，對某點(diǎn)，其影響元素內(nèi)的所有節(jié)點(diǎn)相互影響元素集，兩個節(jié)點(diǎn)和影響元素之交集對節(jié)點(diǎn)和，總體矩陣元素為，由于在一個元素e內(nèi)計算，故稱它為單元矩陣，上式表明矩陣有貢獻(xiàn)的只是()影響元素集所有單元矩陣相應(yīng)的元素。第36頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@374復(fù)雜區(qū)域的譜元方法采用譜逼近的等參數(shù)單元總體坐標(biāo)系下導(dǎo)數(shù)的計算坐標(biāo)變換矩陣及變換行列式局部坐標(biāo)系下偏導(dǎo)數(shù)的計算總體坐標(biāo)系下偏導(dǎo)數(shù)的計算復(fù)雜求解區(qū)域的譜元方法第37頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@38總體坐標(biāo)與局部坐標(biāo)系第38頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@394.1采用譜逼近的等參數(shù)單元總體坐標(biāo)與局部坐標(biāo)(如圖所示)等參單元，插值函數(shù)表達(dá)式和坐標(biāo)變換表達(dá)式一致函數(shù)的插值表達(dá)為，坐標(biāo)變換的表達(dá)式為，導(dǎo)數(shù)表達(dá)式第39頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@404.2總體坐標(biāo)系下導(dǎo)數(shù)的計算坐標(biāo)變換矩陣(Jacobi矩陣)及變換行列式(Jacobi行列式)微元面積、微元弧長

(為常數(shù))時，(為常數(shù))時第40頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@41局部坐標(biāo)系下偏導(dǎo)數(shù)的計算由插值式對ξ求偏導(dǎo)得，上式在(ξm,ηn)處的值為，比照函數(shù)得，則局部坐標(biāo)系下的偏導(dǎo)數(shù)表達(dá)為，第41頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@42總體坐標(biāo)系下偏導(dǎo)數(shù)的計算由插值函數(shù)得，由Jacobi矩陣得總體坐標(biāo)系下，偏導(dǎo)數(shù)計算式為，第42頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@434.3復(fù)雜求解區(qū)域的譜元方法考察Helmholtz方程則元素矩陣為，第43頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@44元素矩陣各變量表達(dá)式第44頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@455非定常不可壓縮流動與自然對流換熱問題求解不可壓縮流動控制方程封閉空腔內(nèi)自然對流換熱基本控制方程解非定常不可壓縮Navier-Stokes方程的時間分裂法時間分裂法求解自然對流問題時間分裂法求解Navier-Stokes方程的誤差分析第45頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@465.1非定常不可壓縮流動控制方程非定常不可壓縮流動無因次Navier-Stokes方程為，以頂蓋移動的方腔驅(qū)動流為例，定義無因次參量如下，第46頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@475.2封閉空腔內(nèi)自然對流換熱基本控制方程封閉空腔內(nèi)自然對流換熱基本控制方程第47頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@48無因次量定義對正方形封閉空腔，無因次參數(shù)的定義如下，對空氣來說，第48頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@495.3解非定常不可壓縮Navier-Stokes方程的時間分裂法Navier-Stokes可寫成如下形式，對時間積分得，設(shè)為三步中間值，則分三步求解的過程為：非線性步、壓力步、粘性步第49頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@50非線性步非線性步使用顯式4階Runge-Kutta方法，替代傳統(tǒng)時間分裂法的Adams-Bashforth方法，使本時間步僅與上一時間步有關(guān)，此步中代入初始條件，不使用邊界條件。第50頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@51壓力步假設(shè)滿足不可壓縮條件，且在邊界上法向上的投影為零，即，得,此步為解帶有第二類邊界條件的壓力Poisson方程，解出壓力Poisson方程后，再由上式求解。第51頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@52粘性步利用Crank-Nicolson格式離散粘性步，得Helmholtz方程，施以固體壁面邊界條件即可求解第52頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@535.4時間分裂法求解自然對流問題求解自然對流問題的過程同Navier-Stokes方程，也是將一個時間步分為三步求解，所不同的是在一個時間步內(nèi)除了解Navier-Stokes外，還需求解能量方程，能量方程中沒有壓力項(xiàng)，因此在第二步壓力Poisson方程求解中與能量方程無關(guān)。第53頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@54非線性步非線性步仍使用顯式4階Runge-Kutta方法，第54頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@55壓力步為解帶有第二類邊界條件的壓力Poisson方程，解出壓力Poisson方程后，再由(5.34)式求解。第55頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@56粘性步利用Crank-Nicolson格式離散粘性步，得Helmholtz方程，施以固體壁面邊界條件即可求解。第56頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@575.5時間分裂法求解Navier-Stokes方程的誤差分析利用時間分裂法求解Navier-Stokes方程基于以下兩個假設(shè)，

滿足不可壓縮約束，即滿足連續(xù)性方程，在邊界的法線上的投影為零，高精度壓力邊界條件描述，壓力Poisson方程變?yōu)椋?7頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@586面向?qū)ο笞V元方法程序設(shè)計面向?qū)ο蟪绦蛟O(shè)計思想面向?qū)ο笞V元方法程序設(shè)計節(jié)點(diǎn)類標(biāo)準(zhǔn)單元類單元類整體區(qū)域類譜元方法程序設(shè)計中的幾個問題插值函數(shù)積分的存儲節(jié)點(diǎn)的編碼及總體矩陣的形成Navier-Stokes方程和能量方程的求解過程第58頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@596.1面向?qū)ο蟪绦蛟O(shè)計思想按人們通常的思維方式建立模型，設(shè)計盡可能自然地表現(xiàn)求解方法的軟件。為此，必須建立直接表現(xiàn)問題中事物與事物相互聯(lián)系的概念，建立適合人們思維方式的描述方法；在面向?qū)ο蟮姆椒ㄖ兄饕幸韵赂拍睿簩ο蟆⑾鬟f、類和繼承、方法等。對象和消息傳遞是表示事物與事物間的關(guān)系；類和繼承是適應(yīng)人們思維方式對事物的描述方法；方法是作用在對象上的各種操作。第59頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@606.2面向?qū)ο笞V元方法程序設(shè)計面向?qū)ο蟮某绦蛟O(shè)計包括面向?qū)ο蠓治?、面向?qū)ο笤O(shè)計和面向?qū)ο髮?shí)現(xiàn)對譜元方法來說，從譜元方法分析問題的思路出發(fā)，找出描述譜元方法的對象類（數(shù)據(jù)和方法的結(jié)合體），建立對象類的關(guān)系，并集成這些對象。設(shè)計節(jié)點(diǎn)類、標(biāo)準(zhǔn)單元、單元類、總體區(qū)域類。這些類中用屬性來描述類的特性，用方法來實(shí)現(xiàn)某些操作。有些類中的屬性可能又需要用對象類來描述，方法的作用對象也可能是某一類的實(shí)體。第60頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@61節(jié)點(diǎn)類主要的數(shù)據(jù)及方法包含節(jié)點(diǎn)的屬性，如節(jié)點(diǎn)的編號、節(jié)點(diǎn)的位置坐標(biāo)、節(jié)點(diǎn)的邊界屬性、微分方程右端項(xiàng)及邊界條件值、所要求解的變量及其導(dǎo)數(shù)值。節(jié)點(diǎn)類定義，classCGridPoint{protected://節(jié)點(diǎn)的屬性標(biāo)準(zhǔn)節(jié)點(diǎn)，編碼，坐標(biāo)，邊界屬性，方程右端項(xiàng)；壓力，速度及其導(dǎo)數(shù)；public://節(jié)點(diǎn)類方法屬性設(shè)置，邊界設(shè)置，壓力和速度設(shè)置；屬性取值，邊界屬性取值；}第61頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@62標(biāo)準(zhǔn)單元類包括了在各求解單元中的常量。類中的數(shù)據(jù)及方法主要有：單元內(nèi)網(wǎng)格的劃分(含兩方向的節(jié)點(diǎn)數(shù)及總節(jié)點(diǎn)數(shù))、標(biāo)準(zhǔn)插值節(jié)點(diǎn)、插值函數(shù)的導(dǎo)數(shù)矩陣及有關(guān)標(biāo)準(zhǔn)插值函數(shù)的積分。標(biāo)準(zhǔn)單元類定義，classCStdElement{protected://標(biāo)準(zhǔn)單元類的屬性網(wǎng)格屬性，標(biāo)準(zhǔn)插值節(jié)點(diǎn)，插值函數(shù)的導(dǎo)數(shù)，插值函數(shù)積分；public://標(biāo)準(zhǔn)單元類方法標(biāo)準(zhǔn)單元類的構(gòu)造，Chebyshev多項(xiàng)式及其導(dǎo)數(shù)；插值函數(shù)的積分，屬性設(shè)置及取值；}第62頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@63單元類主要數(shù)據(jù)及方法包含元素序號、單元內(nèi)節(jié)點(diǎn)劃分(和標(biāo)準(zhǔn)單元內(nèi)劃分相對應(yīng))、標(biāo)準(zhǔn)單元類、節(jié)點(diǎn)類和總體坐標(biāo)與局部坐標(biāo)的變換關(guān)系。單元類定義，classCElement{protected://單元類的屬性單元及網(wǎng)格屬性，標(biāo)準(zhǔn)單元類，坐標(biāo)變換，單元矩陣；public://單元類方法單元構(gòu)造，標(biāo)準(zhǔn)單元類的構(gòu)造，節(jié)點(diǎn)類構(gòu)造；元素矩陣的形成，屬性設(shè)置及取值；}第63頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@64總體區(qū)域類主要數(shù)據(jù)及方法包含區(qū)域內(nèi)單元數(shù)、總節(jié)點(diǎn)數(shù)、節(jié)點(diǎn)類、標(biāo)準(zhǔn)單元類、總體矩陣等?？傮w區(qū)域類定義，classCDomain{protected://整體區(qū)域類的屬性區(qū)域的屬性，單元及網(wǎng)格屬性，單元類，標(biāo)準(zhǔn)單元類；總體矩陣，對角元位置，邊界節(jié)點(diǎn)序號；public://整體區(qū)域類方法單元構(gòu)造，標(biāo)準(zhǔn)單元類的構(gòu)造，節(jié)點(diǎn)類構(gòu)造；元素矩陣的形成，屬性設(shè)置及取值；}第64頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@656.2譜元方法程序設(shè)計中的幾個問題插值函數(shù)積分的存儲節(jié)點(diǎn)的編碼及總體矩陣的形成對稱帶狀矩陣的一維存儲元素的編碼影響元素和影響節(jié)點(diǎn)對角元的位置及一維數(shù)組的長度總剛度矩陣第一類邊界條件的處理第65頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@66插值函數(shù)積分的存儲在形成元素矩陣時，要計算如下積分，

可以看出，

在(i,j,k)變化時，不必占用(Nx+1)3個存儲單元，這樣既節(jié)約了存儲空間又節(jié)省了計算時間。Bijk所占的空間可由下式計算，

Bijk的檢索按如下方式進(jìn)行，先將i,j,k按由大到小排列，設(shè)i>j>k，則Bijk在數(shù)組中的位置為，第66頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@676.3Navier-Stokes方程和能量方程的求解過程構(gòu)造標(biāo)準(zhǔn)單元stdElement，求出對各單元不變的量。構(gòu)造求解域pGeometry，對求解域進(jìn)行網(wǎng)格劃分，同時構(gòu)造節(jié)點(diǎn)類pDot[NG]，構(gòu)造各求解單元pElement[Ne]，對節(jié)點(diǎn)類中U、V、Φ賦初值，賦零值，Φ在熱邊界上賦值1，在冷邊界上賦0，其余兩邊為第二類邊界條件，在形成單元矩陣時自動滿足。利用Runge-Kutta法求解、、，將求解結(jié)果存于、、解壓力Poisson方程，第二類邊界條件，解出壓力后更新、、求解關(guān)于、、的Helmholtz方程，邊界條件第一或第二類。解出Helmholtz方程后更新、、，即完成了一個時間步的求解。第67頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@687計算實(shí)例常微分方程求解Helmholtz和Poisson的求解簡單邊界Helmholtz方程求解復(fù)雜邊界Poisson的求解移動頂蓋方腔驅(qū)動非定常流動封閉方腔中的自然對流第68頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@697.1常微分方程求解【例1】用譜方法求解下面方程該問題的精確解的表達(dá)式為：

將求解域[0,-1]等分為10個求解單元,單元內(nèi)使用5階Chebyshev多項(xiàng)式進(jìn)行譜逼近。數(shù)值計算結(jié)果與精確解的比較如下,比較表明譜元方法具有相當(dāng)高的精度。第69頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@70例1數(shù)值計算結(jié)果與精確解的比較如下序號x數(shù)值解解析解絕對誤差0123456789100.01-3.83295107661418e-0180.1559246312085910.3082929811551880.4535843017646150.5883485562994740.7092408909868030.8130550556913160.8967554397068940.957507402575330.99270559687679100.1559246312085960.3082929811551960.4535843017646250.5883485562994880.7092408909868240.8130550556913390.8967554397069170.9575074025753540.9927055968768061-3.83295107661418e-018-4.88498130835069e-015-8.04911692853239e-015-9.99200722162641e-015-1.45439216225896e-014-2.04281036531029e-014-2.33146835171283e-014-2.38697950294409e-014-2.34257058195908e-014-1.62092561595273e-014-2.22044604925031e-016第70頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@71【例2】【例2】用譜元方法求解下面方程該問題的精確解的表達(dá)式為：

將求解域[0,-1]等分為4個求解單元,單元內(nèi)使用5階Chebyshev多項(xiàng)式進(jìn)行譜逼近。數(shù)值計算結(jié)果與精確解的比較如下,比較表明譜元方法具有相當(dāng)高的精度第71頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@72例2數(shù)值計算結(jié)果與精確解的比較如下0123456789100.01-8.3294087758479e-0190.0147662965508840.0286795455946950.04087816195710.0504833997185180.0565905580149630.0582599254786710.0545073764403450.0442945201363430.0265183090547675.337544225295e-018

00.0147662965612980.0286795455705450.0408781618890680.050483399757920.0565905580149630.0582599255415590.0545073763178480.0442945199587630.0265183091540420

-8.3294087758479e-019-1.041408452529e-0112.4150008381962e-0116.8031476285757e-011-3.9402599238958e-011-3.816391647149e-016-6.2888715068876e-0111.2249663089436e-0101.7757908338245e-010-9.9275830611756e-0115.337544225295e-018

第72頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@737.2簡單邊界Helmholtz方程求解具體考察Helmholtz方程，精確解為，采用矩形單元采用曲邊四邊形單元第73頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@74兩種單元求解結(jié)果比較(簡單邊界)表7-1(采用矩形單元的求解結(jié)果)表7-2(采用曲邊四邊形單元的求解結(jié)果)第74頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@75復(fù)雜邊界Poisson的求解考察方程邊界條件由精確解給出采用環(huán)形段單元采用曲邊四邊形單元第75頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@76兩種單元求解結(jié)果比較(復(fù)雜邊界)表7-3(采用環(huán)形段單元的求解結(jié)果)表7-4(采用曲邊四邊形單元的求解結(jié)果)第76頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@777.3移動頂蓋方腔驅(qū)動非定常流動是計算流體力學(xué)典型數(shù)值實(shí)驗(yàn)?zāi)Ｐ?，在很多情況下都是用它來檢驗(yàn)數(shù)值方法的有效性；除移動頂蓋外，初始速度場為零；移動頂蓋上速度u=1,v=0，其余邊界速度為零；取9個求解單元，使用8階Chebyshev多項(xiàng)式插值；計算了Re=100的情況，取兩種時間步長T1=0.004，T2=0.005，計算2000和1600步，各需12和9.6小時第77頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@78移動頂蓋方腔驅(qū)動流物理模型第78頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@79豎直中心線上的速度U,水平中心線上的速度V圖7-3豎直中心線上的速度U圖7-4水平中心線上的速度V第79頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@80方腔內(nèi)中心點(diǎn)的速度U和速度V隨時間的變化圖7-5方腔內(nèi)中心點(diǎn)的速度U圖7-6方腔內(nèi)中心點(diǎn)的速度V第80頁，共94頁，2023年，2月20日，星期四2023年5月7日Tel：82663537e-mail：glqin@81封閉方腔中的自然對流計算中網(wǎng)格單元數(shù)為3x3，網(wǎng)格內(nèi)部的插值函數(shù)采用8階Che