混合形式泊松方程之FEniCS求解

【我已經發了同名西瓜影片，本文僅是影片中Markdown檔案內容，不包括影片中講解】

如果你希望和我一樣Win10上安裝FEniCS，建議採用Win10 + WSL2 + Ubuntu。

FEniCS講座系列

第1講：有限元法求解偏微分方程之FEniCS入門講解【上一講】

第２講要點

將高階偏微分方程變換成1階偏微分方程組

1階偏微分方程組變換成1個變分方程

自然邊界和基本邊界

如何定義混合空間

如何自定義表示式

FEniCS快捷視覺化

混合形式的泊松方程

以泊松方程為例，這是一個二階方程：

現在我們要將其改成一階方程組，引入一個叫

熱通量

（或簡稱

通量

）的輔助向量

，不難得到：

這個

通量

是有物理意義的，那就是所謂的

傅立葉定律

：熱通量

​與溫度

​的負梯度成比例。

這裡我們取​

而已。

而泊松方程實際就是

傅立葉定律＋能量守恆

的結果。如圖：

穩態情況下，根據能量守恆，透過任意閉合曲面​

流出的能量（這裡指熱量）必須等於曲面內熱源發出的能量：

進而

由於閉合區域​

的任意選擇性，所以在整個​

，下式成立

對應的邊界條件則需要改成：

轉換成變分形式

第一個方程是標量方程，選擇標量測試函式

​；第二個方程是向量方程，選擇向量測試函式

​。

方程１× ​ + 方程２· ​，然後對全域

​積分，得到：

參考這個關係，可以具體地定義測試函式和試探函式。

將所有

測試函式元組

​

的集合記作​

；類似地，將所有

試探函式元組

​的集合

記作​：

請注意，這裡和上一講的不同之處在於， ​

上的第一類邊界條件

​進入了變分公式（現在是

自然邊界條件

），而第二類邊界條件條件

​則在

​上進入試探空間​

的定義（現在是

基本邊界條件

）。

將上面的結果整理得到最終的變分方程：

改成標準變分問題表述：尋求​，滿足

求解這個混合形式泊松方程

​

（一個單位矩陣）

​

（Dirichlet邊界）

​

（Neumann邊界）

​

（法嚮導數）

​

（源項）

第一步：為求解域建立網格，並定義函式空間

from dolfin import *# 建立單位矩形網格mesh = UnitSquareMesh（32， 32）# 定義有限元空間和混合空間BDM = FiniteElement（“BDM”， mesh。ufl_cell（）， 1）DG = FiniteElement（“DG”， mesh。ufl_cell（）， 0）W = FunctionSpace（mesh， BDM * DG）

第二步：定義已提供的表示式

# 第一類界值u0 = Constant（0。0）# 定義源函式f = Expression（“10*exp（-（pow（x［0］ - 0。5， 2） + pow（x［1］ - 0。5， 2）） / 0。02）”， degree=2）# 定義滿足G·n = -g的函式G（實際對應公式中的σ）# g = sin（5x）class BoundarySource（UserExpression）：    def __init__（self， mesh， **kwargs）：        self。mesh = mesh        super（）。__init__（**kwargs）    def eval_cell（self， values， x， ufc_cell）：        cell = Cell（self。mesh， ufc_cell。index）        n = cell。normal（ufc_cell。local_facet）        g = sin（5*x［0］）        values［0］ = -g*n［0］        values［1］ = -g*n［1］    def value_shape（self）：        return （2，）G = BoundarySource（mesh， degree=2）

第三步：建立和應用基本邊界條件

# 定義基本邊界條件def boundary（x）：    return x［1］ < DOLFIN_EPS or x［1］ > 1。0 - DOLFIN_EPSbc = DirichletBC（W。sub（0）， G， boundary）

第四步：定義變分問題

# 定義試探函式和測試函式（sigma， u） = TrialFunctions（W）（tau， v） = TestFunctions（W）# 定義變分形式n = FacetNormal（mesh）a = （div（sigma）*v + dot（sigma， tau） - div（tau）*u）*dxL = f*v*dx - u0*dot（n，tau）*ds

第五步：求解並繪圖

# 求解w = Function（W）solve（a == L， w， bc）（sigma， u） = w。split（）# 繪圖plot（sigma）plot（u）

ParaView視覺化

# 將解儲存為VTK格式vtkfile = File（“mixed_poisson/sigma。pvd”）vtkfile << sigmavtkfile = File（“mixed_poisson/u。pvd”）vtkfile << u

FEniCS快捷視覺化

# 繪製向量場圖（如果是向量場的話）# plot（sigma， mode=“glyphs”）plot（sigma）# 繪製位移圖（如果是向量場的話）plot（sigma， mode=“displacement”）# 等高色圖（如果是標量場的話）# plot（u， mode = “contourf”，levels = 40）# plot（u， mode = “color”，shading = “gouraud”）plot（u）# 標量場3d圖（曲面）plot（u， mode=“warp”， linewidths=0）# 等高線（如果是標量場的話）plot（u， mode=“contour”）

END