博士を目指すカービィ

博士学生の独り言

FEniCSでCFD入門①(非圧縮Navier-Stokes方程式)


Opensourceの有限要素法ベースのCFDソルバーであるFEniCSを使ってみようと思い最近触り始めました.

fenicsproject.org

 

WSLとanacondaを用いてFEniCSのインストールを行いました.

開発環境はWSL+anaconda+jupyter notebookです.

↓インストールの詳細に関しては以下の記事を参照しました.

qiita.com

 

今回は有名なカルマン渦列のシミュレーションを行います.

↓カルマン渦列です(日本機械学会より引用https://www.jsme.or.jp/jsme-medwiki/09:1002355).

 

チュートリアルでカルマン渦列のシミュレーションをしようとしたのですがエラーが多く思うように動作しません.

 

↓参考にしたチュートリアルの記事

fenicsproject.org

記事中ではFenicsの内部関数で形状・メッシュの作成を行っていましたが,以下のようなエラーが出て動きませんでした.

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_9077/1743617236.py in <module>
----> 1 import mshr

~/anaconda3/lib/python3.8/site-packages/mshr/__init__.py in <module>
     22 import dolfin
     23 
---> 24 from .cpp import Circle
     25 from .cpp import Ellipse
     26 from .cpp import Rectangle

ImportError: generic_type: type "CSGGeometry" referenced unknown base type "dolfin::Variable"

↓FEniCS projectのオープンスラックで質問したところ以下のような回答が得られました.

どうやら,FEniCS内部で形状とメッシュを作成するmshrは非推奨となっているようです.外部のメッシャーで構築したメッシュをFEniCS内部に読み込むということでしょうか.

 

今回は面倒だったのでカルマン渦用のメッシュをgithubから拾ってきました.このリポジトリ内のdolfin/data/meshes/内のcylinder.xml.gzを使用します.

github.com

 

↓動作を確認したコードを以下に貼っておきます.

from dolfin import *

# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;

T = 5.0            # final time
num_steps = 5000   # number of time steps
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

mesh = Mesh("cylinder.xml.gz")

# Define function spaces
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)

# Define boundaries
inflow   = 'near(x[0], 0.0)'
outflow  = 'near(x[0], 2.4)'
walls    = 'x[1]<0.0 || near(x[1], 0.4)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

# Define inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx  + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx  + inner(sigma(U, p_n), epsilon(v))*dx  + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds  - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
u_file = File('navier_stokes_cylinder/velocity.pvd')
p_file = File('navier_stokes_cylinder/pressure.pvd')

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')
timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')

# Save mesh to file (for use in reaction_system.py)
File('navier_stokes_cylinder/cylinder.xml.gz') << mesh

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Plot solution
    plot(u_, title='Velocity')
    plot(p_, title='Pressure')

    # Save solution to file (XDMF/HDF5)
    u_file << u_
    p_file << p_

    # Save nodal values to file
    timeseries_u.store(u_.vector(), t)
    timeseries_p.store(p_.vector(), t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    print('u max:', u_.vector().get_local().max())

細かいコードの解釈はまだできていないため,今後解釈でき次第追記していきたいと思います!

 

可視化はおなじみのParaviewで行いました.

↓カルマン渦列が見えますね!

見知らぬ人が作成したメッシュなのでかなり荒いメッシュです.

今度は自力でメッシュ作成からFEniCSで流体計算までしてみようと思います!