熱伝導計算 2 (非定常計算)

2021年6月6日

はじめに

FEniCS で非定常熱伝導計算を行う。

環境

  • FEniCS 2019.2.0
  • Ubuntu 20.04 LTS

問題

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

コード

heat2.py

from dolfin import *

time = 1.
num_steps = 10
dt = time/num_steps

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.)

u_n = Function(V)
u_n.assign(Constant(300.))

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

vtkfile = File("result.pvd")

T = Function(V, name="T")
T.assign(u_n)
vtkfile << T

t = 0.
for n in range(num_steps):
    t += dt

    solve(a == L, T, bc)
    u_n.assign(T)

    vtkfile << T

実行

$ python3 heat2.py

理論的背景

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

$$\frac{\partial u}{\partial t} - \nabla^2 u = f$$

非定常項を差分法で離散化すると、次のようになる。

$$\frac{u^{n+1} - u^n}{\Delta t} - \nabla^2 u^{n+1} = f^{n+1}$$

これに重み付き残差法を適用して、次式を得る。

$$\int_\Omega uv dx + \int_\Omega \Delta t \nabla u \cdot \nabla v dx = \int_\Omega (u^n + \Delta t f^{n+1}v) dx$$

ここで $u = u^{n+1}$ としている。これを時間ステップを進めながら解く。

コードの解説

from dolfin import *

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

time = 1.
num_steps = 10
dt = time/num_steps

時間ステップと時間刻み幅の設定。

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

1 m x 1 m のメッシュを作って関数空間 V を定義する。

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.)

u_n = Function(V)
u_n.assign(Constant(300.))

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

試行関数 $u$ と試験関数 $v$、右辺の関数 $f$ を定義している。ここでは $u^n$ にあたる u_n を作成し、値を設定している。その後、積分方程式の両辺を定義している。理論式そのままである

vtkfile = File("result.pvd")

T = Function(V, name="T")
T.assign(u_n)
vtkfile << T

未知関数 T を定義し、初期値を設定し、VTK ファイルを書き出している。

t = 0.
for n in range(num_steps):
    t += dt

    solve(a == L, T, bc)
    u_n.assign(T)

    vtkfile << T

時間ループを回しながら方程式を解いて、結果を出力する。u_n を更新する必要があることに注意。