熱伝導計算

2021年6月6日

はじめに

FEniCS で熱伝導計算を行う。

環境

  • FEniCS 2019.2.0
  • Ubuntu 20.04 LTS

問題

1 m x 1 m の領域を考え、左端を 300 K、右端を 500 K とした場合の温度分布を求める。

コード

heat.py

from dolfin import *
import matplotlib.pyplot as plt

nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "P", 1)

def bcf_left(x, on_boundary):
    return on_boundary and near(x[0], 0.)

def bcf_right(x, on_boundary):
    return on_boundary and near(x[0], 1.)

bc_left = DirichletBC(V, Constant(300.), bcf_left)
bc_right = DirichletBC(V, Constant(500.), bcf_right)
bc = [bc_left, bc_right]

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)

a = dot(grad(u), grad(v))*dx
L = f*v*dx

T = Function(V, name="T")
solve(a == L, T, bc)

plot(T)
plot(mesh)
plt.show()

vtkfile = File("result.pvd")
vtkfile << T

実行

$ python3 heat.py

plot

ParaView

理論的背景

熱伝導方程式は、次のような形で表される。

$$-\nabla^2 u = f$$

ここで、計算領域を $\Omega$ とし、$x$ は $\Omega$ 上の点として、$u(x)$を未知関数とする。$f(x)$ は規定された関数である。

有限要素法では、重み付き残差式という積分方程式を作る (もともと変分法から出てきたので変分式と呼ばれたりする)。方程式に試験関数 (test function) $v$ を重みとしてかけて積分する。

$$\int_\Omega (\nabla^2 u)v dx = \int_\Omega fv dx$$

これを変形して、部分積分を用いたり、ガウスの発散定理を用いたり、$v$ が境界で 0 であることを用いたりして、結局次式が得られる。

$$\int_\Omega \nabla u \cdot \nabla v dx = \int_\Omega fv dx$$

これを、もとの方程式に対して弱形式 (weak form) という。というのも、$u$ に求められる連続性が 1 段落ちているからである。ここで、試験関数 $v$ の空間を試験空間 (test space) $\hat{V}$、$u$ の空間を試行空間 (trial space) $V$、$u$ を試行関数 (trial function) と呼ぶ。この問題は、ある境界条件において、任意の $v$ に対して上式が成り立つような $u$ を見つける問題になる。有限要素法では、これを離散化して解くことになる。

$V$ と $\hat{V}$ は一般的には分けて書くが、実際には同じものを用いる場合がある。

コードの解説

from dolfin import *
import matplotlib.pyplot as plt

ライブラリのインポート。

nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
V = FunctionSpace(mesh, "P", 1)

1 m x 1 m のメッシュを作って関数空間 V を定義する。ここで、FunctionSpace() の "P" は、普通の内挿関数 (補間関数) を意味する。1 は 1 次式の意味である。

def bcf_left(x, on_boundary):
    return on_boundary and near(x[0], 0.)

def bcf_right(x, on_boundary):
    return on_boundary and near(x[0], 1.)

bc_left = DirichletBC(V, Constant(300.), bcf_left)
bc_right = DirichletBC(V, Constant(500.), bcf_right)
bc = [bc_left, bc_right]

境界条件を作成している。上で定義している関数は、境界であれば True を返す。on_boundary は領域の境界であれば True になる。x には座標が入ってくるので、左端は on_boundary and near(x[0], 0.) となる。near() は等しいという意味である (浮動小数点では厳密に等しくはなりにくいため、近ければ等しいとする)。境界条件は別々に定義して、リストとして結合できる。

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)

a = dot(grad(u), grad(v))*dx
L = f*v*dx

試行関数 $u$ と試験関数 $v$、右辺の関数 $f$ (ここでは 0) を定義している。関数空間は試行空間、試験空間ともに V としている。その後、積分方程式の両辺を定義している。積分記号はないものの、弱形式そのものである。

T = Function(V, name="T")
solve(a == L, T, bc)

未知関数 T を定義し、それと境界条件を指定して、方程式を解いている。

plot(T)
plot(mesh)
plt.show()

matplotlib による計算結果とメッシュの描画。

vtkfile = File("result.pvd")
vtkfile << T

VTK ファイルの書き出し。ParaView で確認できる。