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

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

1 голос
/ 06 мая 2020

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

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

Я очень признателен за вашу помощь.

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

-D (dC_A / dx) = -kC_A

-D (dC_B / dx) = -kC_B

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

1 Ответ

0 голосов
/ 20 мая 2020

Это можно решить как полностью неявную систему. Приведенный ниже код упрощает задачу, чтобы получить единый размер домена и коэффициент диффузии. 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")

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

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

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

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

...