我试图在由gmsh在Fipy中创建的网格中求解Meinhart模型。然而,我不知道如何添加零通量边界条件。下面,你可以找到我的代码。我想知道线路是否:
#u.constrain(0, where=mesh.exteriorFaces)
#u.faceGrad.constrain(mesh.faceCenters, where=mesh.exteriorFaces)
应激活,以确保零通量边界条件
我在:https://github.com/usnistgov/fipy/issues/674中读到“FiPy的边界默认情况下没有通量”。但是,我不知道它是否也适用于由gmsh选项创建的网格
如果不需要激活或向我的代码中添加额外的行来设置零通量边界条件,当我有更复杂的网格(例如多边形、圆或其他不规则形状)时,它也可以工作
谢谢
"""
Emacs Editor
This is a temporary script file.
"""
# -*- coding: utf-8 -*-
"""
@author: Irbin B.
"""
'''Solving Meinhart model 2D'''
# 1. Libraries
import fipy as fi # Finite volume method's package
# 2. Building the domain (Gmsh square)
## 2.1. Domain lenght
nx = 100
ny = nx
## 2.2. Gmsh config
mesh = fi.Gmsh2D('''
Side = 100;
CellSize = 1;
Point(1) = {0, 0, 0, CellSize};
Point(2) = {0, Side, 0, CellSize};
Point(3) = {Side, Side, 0, CellSize};
Point(4) = {Side, 0, 0, CellSize};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(6) = {1, 2, 3, 4};
Plane Surface(7) = {6};
''' % locals())
## 2.3. Adding initial conditions values (Random)
noise_u = fi.GaussianNoiseVariable(mesh=mesh,
mean=0.5,
variance=0.05).value
noise_v = fi.GaussianNoiseVariable(mesh=mesh,
mean=0.5,
variance=0.05).value
# 3. Zero-Flux boundary conditions
BCs = (fi.FixedFlux(faces=mesh.facesRight, value=0.),
fi.FixedFlux(faces=mesh.facesLeft, value=0.),
fi.FixedFlux(faces=mesh.facesTop, value=0.),
fi.FixedFlux(faces=mesh.facesBottom, value=0.))
# 4. Defining the variables
u = fi.CellVariable(name = "u",
mesh=mesh,
value=0.,
hasOld=True)
u[:] = noise_u
#u.constrain(0, where=mesh.exteriorFaces)
#u.faceGrad.constrain(mesh.faceCenters, where=mesh.exteriorFaces)
v = fi.CellVariable(name = "v",
mesh = mesh,
value = 0.,
hasOld = True)
v[:] = noise_v
#v.constrain(0, where=mesh.exteriorFaces)
#v.faceGrad.constrain(mesh.faceCenters, where=mesh.exteriorFaces)
# 5. Defining the parameters
Da = 1.
Db = 100
alpha = -0.005
beta = 10
# 6. Creating the system of PDEs
equ = fi.TransientTerm(var=u) == fi.DiffusionTerm(coeff=Da, var=u) + u - u**3 - v + alpha
eqv = fi.TransientTerm(var=v) == fi.DiffusionTerm(coeff=Db, var=v) + (u - v) * beta
eqn = (equ & eqv)
# 7. Solving the PDEs and showing the results in figures
timeStepDuration = .1
steps = 100
for step in range(steps):
u.updateOld()
v.updateOld()
eqn.sweep(dt = timeStepDuration)
print(step+1)
if __name__== '__main__':
uviewer = fi.Viewer(vars=u, datamin=0., datamax=.7)
vviewer = fi.Viewer(vars=v, datamin=0., datamax=.7)
默认情况下,FiPy的边界条件对于所有网格都不是通量。这是以单元为中心的有限体积法的一个基本特征
相关问题 更多 >
编程相关推荐