1D-связанная переходная диффузия в FiPY с реактивным граничным условием

Я хотел бы решить уравнение нестационарной диффузии для двух соединений A и B, как показано на изображении. Я думаю, что изображение - лучший способ показать мою проблему. Уравнения диффузии и граничные условия.

Как видите, реакция происходит только на поверхности, а поток A равен потоку B. Итак, эти два уравнения связаны только на поверхности. Граничное условие аналогично граничному условию ROBIN, описанному в руководстве Fipy. Однако основным отличием является наличие второй переменной в граничном условии. Кто-нибудь знает, как сформулировать это граничное условие в Fipy? Думаю, мне нужно добавить дополнительный член к граничному условию ROBIN, но я не мог его понять.

Я действительно ценю твою помощь.

Это код, который решает упомянутое уравнение с граничным условием ROBIN @ x = 0.

-D(dC_A/dx) = -kC_A

-D(dC_B/dx) = -kC_B

В этом условии я могу легко использовать граничное условие ROBIN для решения уравнений. Результаты кажутся разумными для этого граничного условия.

"""
Question for StackOverflow
"""
#%%
from fipy import Variable, FaceVariable, CellVariable, Grid1D, TransientTerm, DiffusionTerm, Viewer, ImplicitSourceTerm
from fipy.tools import numerix



#%%
##### Model parameters

L= 8.4853e-4 # m boundary layer thickness
dx= 1e-8 # mesh size 
nx = int(L/dx)+1 # number of meshes
D = 1e-9 # m^2/s diffusion coefficient
k = 1e-4 # m/s reaction coefficient R = k [c_A],
c_inf =  0. # ROBIN general condition, once can think R = k ([c_A]-[c_inf])
c_init = 1. # Initial concentration of compound A, mol/m^3


#%%
###### Meshing and variable definition

mesh = Grid1D(nx=nx, dx=dx)
c_A = CellVariable(name="c_A", hasOld = True,
                    mesh=mesh,
                    value=c_init)
c_B = CellVariable(name="c_B", hasOld = True,
                    mesh=mesh,
                    value=0.)

#%%
##### Right boundary condition 

valueRight = c_init
c_A.constrain(valueRight, mesh.facesRight)
c_B.constrain(0., mesh.facesRight)

#%%
### ROBIN BC requirements, defining cellDistanceVectors
## This code is for fixing celldistance via this link:
## https://stackoverflow.com/questions/60073399/fipy-problem-with-grid2d-celltofacedistancevectors-gives-error-uniformgrid2d
MA = numerix.MA
tmp = MA.repeat(mesh._faceCenters[..., numerix.NewAxis,:], 2, 1)
cellToFaceDistanceVectors = tmp - numerix.take(mesh._cellCenters, mesh.faceCellIDs, axis=1)
tmp = numerix.take(mesh._cellCenters, mesh.faceCellIDs, axis=1)
tmp = tmp[..., 1,:] - tmp[..., 0,:]
cellDistanceVectors = MA.filled(MA.where(MA.getmaskarray(tmp), cellToFaceDistanceVectors[:, 0], tmp))

#%%
##### Defining mask and Robin BC at left boundary 
mask = mesh.facesLeft
Gamma0 = D
Gamma = FaceVariable(mesh=mesh, value=Gamma0)
Gamma.setValue(0., where=mask)
dPf = FaceVariable(mesh=mesh,
                   value=mesh._faceToCellDistanceRatio * cellDistanceVectors)
n = mesh.faceNormals
a = FaceVariable(mesh=mesh, value=k, rank=1)
b = FaceVariable(mesh=mesh, value=D, rank=0)
g = FaceVariable(mesh=mesh, value= k * c_inf, rank=0)
RobinCoeff = (mask * Gamma0 * n / (-dPf.dot(a)+b))

#%%
#### Making a plot
viewer = Viewer(vars=(c_A, c_B),
                     datamin=-0.2, datamax=c_init * 1.4)
viewer.plot()

#%% Time step and simulation time definition
time = Variable()
t_simulation = 4 # seconds
timeStepDuration = .05
steps = int(t_simulation/timeStepDuration)

#%% PDE Equations
eqcA = (TransientTerm(var=c_A) == DiffusionTerm(var=c_A, coeff=Gamma) + 
            (RobinCoeff * g).divergence 
            - ImplicitSourceTerm(var=c_A, coeff=(RobinCoeff * a.dot(-n)).divergence))

eqcB = (TransientTerm(var=c_B) == DiffusionTerm(var=c_B, coeff=Gamma) -
                (RobinCoeff * g).divergence 
            + ImplicitSourceTerm(var=c_B, coeff=(RobinCoeff * a.dot(-n)).divergence))


#%% A loop for solving PDE equations
while time() <= (t_simulation):
    time.setValue(time() + timeStepDuration)
    c_B.updateOld()
    c_A.updateOld()
    res1=res2 = 1e10
    viewer.plot()
    while (res1 > 1e-6) & (res2 > 1e-6):
        res1 = eqcA.sweep(var=c_A, dt=timeStepDuration)
        res2 = eqcB.sweep(var=c_B, dt=timeStepDuration)

pde fipy
person Ali    schedule 06.05.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Это можно решить как полностью неявную систему. Приведенный ниже код упрощает задачу, чтобы получить единый размер домена и коэффициент диффузии. k установлен на 0,2. Он довольно хорошо передает аналитическое решение с некоторыми оговорками (см. Ниже).

from fipy import (
    CellVariable,
    TransientTerm,
    DiffusionTerm,
    ImplicitSourceTerm,
    Grid1D,
    Viewer,
)

L = 1.0
nx = 1000
dx = L / nx
konstant = 0.2
coeff = 1.0

mesh = Grid1D(nx=nx, dx=dx)

var_a = CellVariable(mesh=mesh, value=1.0, hasOld=True)
var_b = CellVariable(mesh=mesh, value=0.0, hasOld=True)

var_a.constrain(1.0, mesh.facesRight)
var_b.constrain(0.0, mesh.facesRight)

coeff_mask = ~mesh.facesLeft * coeff
boundary_coeff = konstant * (mesh.facesLeft * mesh.faceNormals).divergence

eqn_a = TransientTerm(var=var_a) == DiffusionTerm(
    coeff_mask, var=var_a
) - ImplicitSourceTerm(boundary_coeff, var=var_a) + ImplicitSourceTerm(
    boundary_coeff, var=var_b
)

eqn_b = TransientTerm(var=var_b) == DiffusionTerm(
    coeff_mask, var=var_b
) - ImplicitSourceTerm(boundary_coeff, var=var_b) + ImplicitSourceTerm(
    boundary_coeff, var=var_a
)

eqn = eqn_a & eqn_b

for _ in range(5):
    var_a.updateOld()
    var_b.updateOld()
    eqn.sweep(dt=1e10)

Viewer((var_a, var_b)).plot()

print("var_a[0] (expected):", (1 + konstant) / (1 + 2 * konstant))
print("var_b[0] (expected):", konstant / (1 + 2 * konstant))
print("var_a[0] (actual):", var_a[0])
print("var_b[0] (actual):", var_b[0])

input("wait")

Обратите внимание на следующее:

  • Как написано, граничные условия имеют точность только первого порядка, что на самом деле не имеет значения для этой проблемы, но может навредить вам в более высоких измерениях. Могут быть способы исправить это, например, наличие небольшой ячейки рядом с границей или добавление явной поправки второго порядка для граничного условия.

  • Уравнения здесь связаны. В случае разъединения для достижения равновесия, вероятно, потребовалось бы множество итераций.

  • Для достижения равновесия потребовалось несколько итераций, но этого не должно быть. Вероятно, это связано с тем, что решающая программа не сходится должным образом без нескольких попыток. Возможно, у связанных уравнений плохая обусловленность.

person wd15    schedule 20.05.2020