物質拡散計算

2021年6月6日

はじめに

FEniCS で物質拡散計算を行う。

環境

• FEniCS 2019.2.0
• Ubuntu 20.04 LTS

問題

2 m x 2 m (中心が原点) の領域を考え、初期に中央 0.4 m 角に濃度 1 を置いた場合の濃度分布を求める。

コード

diffusion.py

from dolfin import *

time = 0.5
num_steps = 50
dt = time/num_steps

nx = ny = 50
mesh = RectangleMesh(Point(-1., -1.), Point(1., 1.), nx, ny)
V = FunctionSpace(mesh, "P", 1)

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

u_0 = Expression("abs(x[0]) < 0.2 && abs(x[1]) < 0.2 ? 1. : 0.", degree=1)
u_n = interpolate(u_0, V)

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

vtkfile = File("result.pvd")

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

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

    solve(a == L, u)
    u_n.assign(u)

    vtkfile << u

実行

$ python3 diffusion.py

理論的背景

非定常熱伝導の場合と同じ。温度を濃度と読み換えるだけ。

コードの解説

nx = ny = 50
mesh = RectangleMesh(Point(-1., -1.), Point(1., 1.), nx, ny)
V = FunctionSpace(mesh, "P", 1)

(-1, -1) - (1, 1) の範囲の矩形領域メッシュを作る。

このあと境界条件の設定をしていないが、これは境界を勾配 0 条件にしたいためである。

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

u_0 = Expression("abs(x[0]) < 0.2 && abs(x[1]) < 0.2 ? 1. : 0.", degree=1)
u_n = interpolate(u_0, V)

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

初期条件を式で定義している。この式は C++ の式である。

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

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

    solve(a == L, u)
    u_n.assign(u)

    vtkfile << u

熱伝導計算とほぼ同じ。