數(shù)值求解波動(dòng)方程 [3]

2022-08-23 15:19 作者:nyasyamorina 0人讀過 | 我要投稿

完美匹配層

????????之前提到的幾種截?cái)喾绞???Dirichlet 邊界,??Neumann 邊界和循環(huán)邊界,? 這幾種截?cái)喾绞蕉加幸粋€(gè)共同特征:? 就是波會(huì)在邊界處反射回來 (循環(huán)邊界可以看作一邊的波被"反射"到另一邊).? 但對(duì)于實(shí)際應(yīng)用來說,? 很多情況下要求邊界不能產(chǎn)生反射,? 以避免反射波對(duì)內(nèi)部區(qū)域造成影響.? 因?yàn)椴淮嬖诜瓷洳?? 對(duì)于內(nèi)部區(qū)域來說,? 波好像消散在無限空間里.

????????有一種實(shí)現(xiàn)無反射"邊界"的技術(shù)叫完美匹配層 (perfectly matched layer),? 簡稱 PML.? 準(zhǔn)確來說這不是一個(gè)邊界條件,? 而是一個(gè)靠近邊界的區(qū)域.? 在對(duì) PML 進(jìn)行展開細(xì)說之前,? 首先要對(duì)波動(dòng)方程進(jìn)行開刀.

? ? ? ? 引入一個(gè)新量 $%5Cvec%20v(x_1%2Cx_2%2C%5Ccdots%2Ct)$ ?使得 $%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3A%3D%5Cvec%5Cnabla%20u$ ,? 代入波動(dòng)方程方程? $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cvec%5Cnabla%5E2u$ 得? $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20t%7D%5Cvec%5Cnabla%5Ccdot%5Cvec%20v$ ,? 即? $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cvec%20v$ .

????????如此一元二階微分方程? $%5Cfrac%7B%5Cpartial%5E2%20u%7D%7B%5Cpartial%20t%5E2%7D%3Dc%5E2%5CDelta%20u$ ?分解為二元一階微分方程組 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B%5Cpartial%5Cvec%20v%7D%7B%5Cpartial%20t%7D%3D%5Cvec%5Cnabla%20u%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cvec%5Cnabla%5Ccdot%5Cvec%20v%5Cend%7Bmatrix%7D%5Cright.$ .

????????題外話:? 可以記 $%5Cvec%20w%3D%5Cbegin%7Bbmatrix%7Du%5C%5C%5Cvec%20v%5Cend%7Bbmatrix%7D$ ,? 則有 $%5Cfrac%7B%5Cpartial%5Cvec%20w%7D%7B%5Cpartial%20t%7D%3D%5Chat%20D%5Cvec%20w$ ,? 其中 $%5Cvec%20D%3D%5Cbegin%7Bbmatrix%7D%26c%5E2%5Cvec%5Cnabla%5Ccdot%5C%5C%5Cvec%5Cnabla%26%5Cend%7Bbmatrix%7D$ ?是反厄米算符.? 一般地,? 如果一個(gè)方程可以寫為 $%5Cfrac%7B%5Cpartial%5Cvec%20w%7D%7B%5Cpartial%20t%7D%3D%5Chat%20D%5Cvec%20w%0A$ ?形式,? 并且 $%5Chat%20D$ ?為反厄米算符,? 那么 $%5Cvec%20w$ ?就會(huì)產(chǎn)生 "類波" 解:? 麥克斯韋方程組,? 薛定諤方程,??Lamé-Navier 方程都符合.? 這也是為什么薛定諤方程又叫薛定諤波動(dòng)方程,? 有些人別整天揪著 "薛定諤波動(dòng)方程" 不放了,? 是薛定諤 "波動(dòng)" 方程,? 不是薛定諤 "波動(dòng)方程".

????????PML 背后的思想就是變換坐標(biāo),? 來看看變換坐標(biāo)如何做到 "吸收" 波.

????????定義坐標(biāo)變換? $x%5Cmapsto%5Ctilde%20x%3Dx%2Bif(x)$ ,? 其中 $%5Cforall%20x%3A%5C%3Af(x)%5Cgeq0$ .??當(dāng)? $f(x)%3D0$ 時(shí),? 坐標(biāo)變換為自身,? 那么? $e%5E%7Bikx%7D%0A$ 的圖像可以展示為

????????改變? $f(x)$ ?后可以得到:

????????可以看到波在? $f(x)%3E0$ ?的區(qū)域內(nèi)呈指數(shù)下降,? 如同波在這個(gè)區(qū)域內(nèi)被 "吸收" 了一樣.? 實(shí)際上,?? $%20e%5E%7Bikx%7D$ ?在施加坐標(biāo)變換 $x%5Cmapsto%5Ctilde%20x%3Dx%2Bif(x)$ ?得到 $e%5E%7Bikx%7D%5Cmapsto%20e%5E%7Bik(x%2Bif(x))%7D%3De%5E%7Bikx%7De%5E%7B-kf(x)%7D$ ,? 其中項(xiàng)? $e%5E%7B-kf(x)%7D$ ?即表現(xiàn)為指數(shù)衰減.

????????當(dāng)使用 PML 包圍內(nèi)部區(qū)域,? 并且在 PML 外進(jìn)行截?cái)?? 既然波在 PML 區(qū)域內(nèi)呈指數(shù)衰減,? 那么只要 PML 足夠厚,? 波在到達(dá)邊界時(shí)可以變得足夠小,? 即使波在邊界處被反射,? 反射波也無法在內(nèi)部區(qū)域產(chǎn)生足夠大的影響.

????????PML 內(nèi)使用的坐標(biāo)變換稱為 PML 變換.? PML 變換與上面展示的坐標(biāo)變換類似:? 由? $%5Ctilde%20x%3Dx%2Bif(x)$ ?得到 $d%5Ctilde%20x%3D%5Cleft(1%2Bi%5Cfrac%7Bdf%7D%7Bdx%7D%5Cright)dx$ .? 在 PML 變換里記 $%5Cfrac%7Bdf%7D%7Bdx%7D%3D%5Csigma(x)$ ,? 但坐標(biāo)變換為 $d%5Ctilde%20x%3D%5Cleft(1%2Bi%5Cfrac%7B%5Csigma(x)%7D%7B%5Comega%7D%5Cright)dx$ ,? 引入頻率 $%5Comega%0A$ ?的原因是上面求得的衰減項(xiàng) $e%5E%7B-kf(x)%7D$ ?與波數(shù) $k$ 有關(guān),? 而頻率和波數(shù)又有? $%5Comega%3Dc%7Ck%7C$ ,? 也就是說衰減速度于頻率成正比,? 這會(huì)造成低頻分量無法有效地被吸收.??在引入頻率后對(duì)應(yīng)的衰減項(xiàng)為? $e%5E%7B-%5Cfrac%7B1%7D%7Bc%7D%5Cint%5Ex%5Csigma(x)dx%7D$ ,? 衰減速度與頻率無關(guān).? 為了確保波在 PML 內(nèi)總是衰減的,? 應(yīng)該有 $%5Cforall%20x%3A%5C%3A%5Csigma(x)%5Cgeq0$ .

????????下面以一維做例子,? 二元一階微分方程組版的一維波動(dòng)方程為 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%5Cend%7Bmatrix%7D%5Cright.$ .? 初值條件? $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Du(x%2C0)%5C%5C%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(x%2C0)%5Cend%7Bmatrix%7D%5Cright.$ ?可以根據(jù)關(guān)系變?yōu)? $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Du(x%2C0)%5C%5Cv(x%2C0)%3D%5Cfrac%7B1%7D%7Bc%5E2%7D%5Cint%5Ex%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D(x'%2C0)dx'%5Cend%7Bmatrix%7D%5Cright.$ .

????????計(jì)算過程以 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D$ ?為例:? 在 PML 內(nèi)有? $v%3De%5E%7Bi(kx-%5Comega%20t)%7D$ ?和 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D-i%5Comega%20v$ ,? 即 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%3D-i%5Comega%20v$ .? 施加 PML 變換? $%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D%5Cmapsto%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%5Ctilde%20x%7D%3D%5Cfrac%7B1%7D%7B1%2Bi%5Cfrac%7B%5Csigma%7D%7B%5Comega%7D%7D%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial%20x%7D$ ,? 得到? $%5Cfrac%7B1%7D%7B1%2Bi%5Cfrac%7B%5Csigma%7D%7B%5Comega%7D%7D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D%3D-i%5Comega%20v$ ,??然后得到? $-i%5Comega%20v%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D-%5Csigma%20v$ ,? 即 $%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20t%7D%3D%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20x%7D-%5Csigma%20v$ .

????????同理可得? $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D%3Dc%5E2%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D-%5Csigma%20u$ .? 在內(nèi)部區(qū)域里,? 應(yīng)該有 $%5Csigma%3D0$ ,? 這時(shí) PML 版本波動(dòng)方程退化為普通的波動(dòng)方程,? 而在 PML 區(qū)域內(nèi)則有 $%5Csigma%3E0$ .

????????這里選取? $%5Cfrac%7Bu_%7Bx%2Ct%2B1%7D-u_%7Bx%2Ct%7D%7D%7B%5CDelta%20t%7D$ ?作為 $%5Cfrac%7B%5Cpartial%20u%7D%7B%5Cpartial%20t%7D$ ?的離散化也不會(huì)造成太大問題 (節(jié)省內(nèi)存),? 那么 PML 版波動(dòng)方程離散化為 $%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Dv_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7B%5CDelta%20t%7D%7B2%5CDelta%20x%7D(u_%7Bx%2B1%2Ct%7D-u_%7Bx-1%2Ct%7D)%2B(1-%5Csigma_x%5CDelta%20t)v_%7Bx%2Ct%7D%5C%5Cu_%7Bx%2Ct%2B1%7D%3D%5Cfrac%7Bc_x%5E2%5CDelta%20t%7D%7B2%5CDelta%20x%7D(v_%7Bx%2B1%2Ct%7D-v_%7Bx-1%2Ct%7D)%2B(1-%5Csigma_x%5CDelta%20t)u_%7Bx%2Ct%7D%5Cend%7Bmatrix%7D%5Cright.$ .? 并且這種實(shí)現(xiàn)不存在數(shù)值穩(wěn)定的參數(shù)選擇,? 最多只能做到較小的時(shí)間數(shù)值不穩(wěn)定.

????????需要注意的是,? 在連續(xù)波動(dòng)方程里,? PML 是無論? $%5Csigma$ ?如何選取都是沒有反射波的,? 但離散化后則不是這樣,? 為了確保離散化后的 PML 仍然沒有反射,? 需要滿足 $%5Csigma$ ?具有連續(xù)的二階或三階導(dǎo)數(shù).

????????這里介紹一個(gè)常用具有連續(xù)二階導(dǎo)數(shù)的? $%5Csigma$ :? 以右邊界為例,? 假設(shè) $x%5Cin(0%2C1)$ ?為 PML,? 那么可以選取? $%5Csigma(x)%3Ds%5Cleft(x-%5Cfrac%7B%5Csin(2%5Cpi%20x)%7D%7B2%5Cpi%7D%5Cright)$ ,? 其中 $s$ ?為 $%5Csigma$ ?的最大值,? 對(duì)于不同位置不同厚度的 PML 可以使用坐標(biāo)變換把這個(gè) σ 映射過去.??剛剛說過,? 盡管 PML 可以快速吸收波,? 但截?cái)嗟倪吔缛詴?huì)產(chǎn)生微小的反射波返回內(nèi)部區(qū)域,? 當(dāng)選取 Dirichlet 或 Neumann 邊值條件時(shí),? 反射波的相對(duì)大小由 $R%3De%5E%7B-%5Cfrac%7BsL%7D%7Bc%7D%7D$ 給出,? 其中 $L$ 為 PML 的厚度.? 那么 PML 系統(tǒng)如下圖所示: