# How to model diffusion through a membrane?

This is a follow-up on How to handle discontinuity in diffusion coefficient?

Consider diffusion of $$u(t,x)$$ on the domain $$x \in [0,2]$$ with some simple boundary conditions such as $$u(0) = 2, u(2) = 1$$.

Our domain is split into two parts: $$[0,1)$$ on the left and $$(1,2]$$ on the right, with different diffusion coefficients, e.g. $$D^\text{left} = 1, D^\text{right} = 3$$.

The diffusion equation is: $$\partial_t u = \partial_x (D \partial_x u)$$

So far, this is the summary of the linked question.

This time we also have a membrane at $$x=1$$, imposing the following condition on the fluxes at $$x=1$$: $$D^\text{left} \partial_x u^\text{left} = D^\text{right} \partial_x u^\text{right} = d^\text{membrane} (u^\text{right} - u^\text{left})$$

What is the cleanest way to model this with NDSolve? Is there a way to preserve the sharp conditions at $$x=1$$? Perhaps one approximation that could be used is to consider a membrane of finite thickness, having a very high diffusion coefficient of its own. However, this is really a hack. Is it possible to solve the equation on the two half-domains "separately" and couple the boundary conditions at $$x=1$$?

• I don't think this can be solved only with NDSolve. Using my pdetoode to discretize should be relatively straightforward though. Is an answer with pdetoode desired? – xzczd Apr 7 at 12:34
• @Szabolcs You did not specify the initial data. – Alex Trounev Apr 7 at 12:53
• @AlexTrounev Do you mean $u(t=0,x)$? I assumed it does not matter for the solution method what it is. You can make up an initial condition, e.g. $u(0,x) = 2$ if $x < 1$ and $u(0,x) = 1$ if $x>1$. – Szabolcs Apr 7 at 12:59
• @user21 In many cases it is orders of magnitude thinner though. But that is not really the main reason to treat it as infinitely thin. Usually, a physicist would want to understand (and as part of understanding, numerically solve) a simplified model. The goal is not to model a real system as accurately as possible (that's for engineers who work with systems where the physics is well understood already), but to understand the behaviour of a mathematical model. Our starting point is the equations. Giving it a thickness is a deviation from those, and in this sense, a hack. – Szabolcs Apr 8 at 8:22
• @Szabolcs, for anything beyond 1D, yes. Especially for complicated geometries. One might be able to work around this with some low level fem programming and program the interface condition yourself. I for one am not sure how to do it. IMHO Tim's approach is the only one that generalizes. – user21 Apr 9 at 13:17

This answer is a partial response to a comment about extending the approach to more complex geometry. The preliminary results seemed encouraging so I thought I would share my workflow.

I think there are times where one might want to model the membrane region due to difficulties imposing internal boundary conditions. For chemical interphase mass transfer, there can be discontinuities in both coefficients and also the field variable due to phase changes. The characteristic length and timescales of interfacial phenomena are so small, that they are generally assumed to be in equilibrium leading to jumps in concentrations as shown in the following figures.  To use the FEM method it is good to cast your equations into coefficient form as shown FEM Tutorial.

$$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

By doing so, we can use region IDs to toggle equations to be active in some regions and suppressed in others.

For interfacial chemical equilibria, we toggle a source term in the interface region that drives the phase concentrations to their equilibrium values. I posted an article about Modeling jump conditions in interphase mass transfer on the Wolfram Community. In the end, the modeling a thin interface region compared favorably to commercial codes that had support for internal boundary conditions.

What I present here is an approach based on the Acoustic Cloak Monograph to provide an efficient quad mesh for the interface.

# Interface Modeling

Since the interface is a small feature, model sizes can grow to be very large if one tries to isotropically mesh the region. The Acoustic Cloak Monograph, uses high aspect QuadElements to get around this problem. I also make use of the Tensor Product Grid example in the RegionProduct documentation to create 2D regions.

# Helper Functions

I had some difficulty combining some multiple Tri regions with Quad regions so I hacked some functions together. There is probably a better way to do this, but it seems to work.

Needs["NDSolveFEM"];
(* From RegionProduct Documentation *)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
(* Convert RegionProduct to ElementMesh *)
rp2Mesh[rh_, rv_, marker_] := Module[{sqr, crd, inc, msh, mrkrs},
sqr = RegionProduct[rh, rv];
crd = MeshCoordinates[sqr];
inc = Delete /@ MeshCells[sqr, 2];
mrkrs = ConstantArray[marker, First@Dimensions@inc];
msh = ToElementMesh["Coordinates" -> crd,
]
(* Create an annular ElementMesh *)
annularMesh[r_, th_, rh_, rv_, marker_] :=
Module[{r1, r2, th1, th2, anMesh, crd, melms, newcrd},
{r1, r2} = r;
{th1, th2} = th;
anMesh = rp2Mesh[rh, rv, marker];
crd = anMesh["Coordinates"];
melms = anMesh["MeshElements"];
newcrd =
Chop[{#1 Cos[#2], #1 Sin[#2]} & @@@ ({r1 + (r2 - r1) #1,
th1 + (th2 - th1) #2} & @@@ crd), 1*^-7];
ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms]
]
(* Combine and Flatten 2 Tri Meshes *)
combineTriMeshes[mesh1_, mesh2_] :=
Module[{crd1, crd2, newcrd, numinc1, inc, inc1, inc2, mrk, mrk1,
mrk2, elm1, elm2, melms, m},
crd1 = mesh1["Coordinates"];
crd2 = mesh2["Coordinates"];
numinc1 = First@Dimensions@crd1;
newcrd = crd1~Join~ crd2;
inc1 =  ElementIncidents[mesh1["MeshElements"]][];
inc2 =  numinc1 + ElementIncidents[mesh2["MeshElements"]][];
mrk1 = ElementMarkers[mesh1["MeshElements"]][];
mrk2 = ElementMarkers[mesh2["MeshElements"]][];
melms = {TriangleElement[inc1~Join~inc2, mrk1~Join~mrk2]};
m = ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms];
m
]
(* Combine and Flatten 2 Quad Meshes *)
Module[{crd1, crd2, newcrd, numinc1, inc, inc1, inc2, mrk, mrk1,
mrk2, elm1, elm2, melms, m},
crd1 = mesh1["Coordinates"];
crd2 = mesh2["Coordinates"];
numinc1 = First@Dimensions@crd1;
newcrd = crd1~Join~ crd2;
inc1 =  ElementIncidents[mesh1["MeshElements"]][];
inc2 =  numinc1 + ElementIncidents[mesh2["MeshElements"]][];
mrk1 = ElementMarkers[mesh1["MeshElements"]][];
mrk2 = ElementMarkers[mesh2["MeshElements"]][];
m = ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms];
m
]
(* Combine Mixed Quad and Tri Mesh *)
combineMeshes[mesh1_, mesh2_] :=
Module[{crd1, crd2, newcrd, numinc1, inc1, inc2, mrk1, mrk2, elm1,
elm2, melms, m},
crd1 = mesh1["Coordinates"];
crd2 = mesh2["Coordinates"];
numinc1 = First@Dimensions@crd1;
newcrd = crd1~Join~ crd2;
inc1 =  ElementIncidents[mesh1["MeshElements"]][];
inc2 =  ElementIncidents[mesh2["MeshElements"]][];
mrk1 = ElementMarkers[mesh1["MeshElements"]] // Flatten;
mrk2 = ElementMarkers[mesh2["MeshElements"]] // Flatten;
elm1 = mesh1["MeshElements"][[1, 0]];
elm2 = mesh2["MeshElements"][[1, 0]];
melms = Flatten@{elm1[inc1, mrk1], elm2[inc2 + Length[crd1], mrk2]};
m = ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms];
m = MeshOrderAlteration[m, 2];
m
]


# Build a Mixed Element Mesh

Here we will create a half symmetry model of an inner liquid drop, surrounded by a membrane (meshed with flat quads), and gas.

(* Define regions association for clearer assignment *)
regs = <|"inner" -> 10, "membrane" -> 20, "outer" -> 30|>;
(* Geometry Parameters *)
thick = rinner/100.;
rinner = 1.0;
router = rinner + thick;
rmax = 2 rinner;
theta = 180 Degree;
afrac = theta/(360 Degree);
(* Define Mesh Levels *)
nAngular = 120;
(* Elements across the thickness of the membrane *)
(* Angular resolution *)
rv = pointsToMesh[Subdivide[0, 1, nAngular afrac]];
(* Create Membrane Quad Mesh *)
membraneMesh =
annularMesh[{rinner, router}, {0 Degree, 180 Degree}, rh, rv,
regs["membrane"]];
membraneMesh[
"Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> Orange]]
(* Create inner drop mesh based on membraneMesh *)
bmeshinner =
ToBoundaryMesh[
Rectangle[{-rinner, 0}, {rinner, (rinner + router)/2}],
"MaxBoundaryCellMeasure" -> rinner/20];
coordinates =
Join[Select[membraneMesh["Coordinates"], Norm[#] <= rinner &],
Select[bmeshinner["Coordinates"], #[] == 0 &]];
incidents = Partition[FindShortestTour[coordinates][], 2, 1];
innerBoundary =
ToBoundaryMesh["Coordinates" -> coordinates,
"BoundaryElements" -> {LineElement[incidents]}];
innerMesh =
ToElementMesh[innerBoundary, "MeshOrder" -> 1,
"MaxCellMeasure" -> 0.01/4, "SteinerPoints" -> False,
"RegionMarker" -> {{0, rinner/2}, regs["inner"]}];
innerMesh[
"Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> Blue]]

(* Create outer drop mesh based on membraneMesh *)
bmeshouter =
ToBoundaryMesh[Rectangle[{-rmax, 0}, {rmax, rmax}],
"MaxBoundaryCellMeasure" -> rinner/20];
coordinates =
Join[Select[membraneMesh["Coordinates"], Norm[#] >= router &],
Select[bmeshouter["Coordinates"], #[] == 0 &]];
coordinates =
Join[Select[membraneMesh["Coordinates"], Norm[#] >= router &],
Select[
bmeshouter[
"Coordinates"], (! ((-router <= #[] <= router) && #[] ==
0)) &]];
incidents = Partition[FindShortestTour[coordinates][], 2, 1];
outerBoundary =
ToBoundaryMesh["Coordinates" -> coordinates,
"BoundaryElements" -> {LineElement[incidents]}];
outerMesh =
ToElementMesh[outerBoundary, "MeshOrder" -> 1,
"MaxCellMeasure" -> 0.01/4, "SteinerPoints" -> False,
"RegionMarker" -> {{0, (rmax + router)/2}, regs["outer"]}];
outerMesh[
"Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> Green]]

(* Combine Meshes into one *)
mesh = combineTriMeshes[innerMesh, outerMesh];
mesh = combineMeshes[mesh, membraneMesh];
mesh["Wireframe"[
"MeshElementStyle" -> (FaceForm[#] & /@ {Blue, Green, Orange})]]
mesh["Wireframe"[
PlotRange -> {{-rmax/1.75, -router + 0.25}, {0, 0.25}},
"MeshElementStyle" -> (FaceForm[#] & /@ {Blue, Green, Orange})]]  # Set up and Solve Three Region PDE

After creating a 2D mesh, we set up our system of PDEs for gas and liquid concentrations. Note that we introduce small diffusion coefficient, $$dsmall$$, to prevent species from leaking past the membrane.

For simplicity, we will initialize the system at zero concentration and use a Dirichlet condition of 1 for gas concentration on the left wall.

(* Inner Region *)
d1 = 0.1;
(* Outer Region *)
d2 = 3 d1;
(* Membrane Region *)
d3 = 10 d2;
dsmall = d1/10000;
(* Region Dependent Parameters *)
(* Diffusion Coeffiecients *)
di =  With[{d1 = d1, d2 = d2, d3 = d3, dsmall = dsmall},
Piecewise[{{DiagonalMatrix@{d1, d1},
ElementMarker ==
regs["inner"]}, {DiagonalMatrix@{dsmall, dsmall},
ElementMarker == regs["outer"]}, {DiagonalMatrix@{d3, d3},
True}}]];
do =  With[{d1 = d1, d2 = d2, d3 = d3, dsmall = dsmall},
Piecewise[{{DiagonalMatrix@{dsmall, dsmall},
ElementMarker == regs["inner"]}, {DiagonalMatrix@{d2, d2},
ElementMarker == regs["outer"]}, {DiagonalMatrix@{d3, d3},
True}}]];
(* Toggle Source Terms for Interface *)
kappa = 1;
omega = Evaluate[If[ElementMarker == regs["membrane"], kappa, 0]];
kequil = 0.5;
c0 = 1;
tmax = 30;

dcli = DirichletCondition[ui[t, x, y] == 0, x == -rmax];
dcri = DirichletCondition[ui[t, x, y] == 0, x == rmax];
dclo = DirichletCondition[uo[t, x, y] == c0, x == -rmax];
dcro = DirichletCondition[uo[t, x, y] == 0, x == rmax];
ics = {ui[0, x, y] == 0, uo[0, x, y] == 0};
eqni =  D[ui[t, x, y], t] +
Inactive[Div][-di.Inactive[Grad][ui[t, x, y], {x, y}], {x, y}] +
omega (kequil ui[t, x, y] - uo[t, x, y]) == 0;
eqno = D[uo[t, x, y], t] +
Inactive[Div][-do.Inactive[Grad][uo[t, x, y], {x, y}], {x, y}] -
omega (kequil ui[t, x, y] - ui[t, x, y]) == 0;
pdes = {eqni, eqno};
uif = NDSolveValue[
pdes~Join~{dcli, dcri, dclo, dcro}~Join~ics, {ui, uo}, {t, 0,
tmax}, {x, y} \[Element] mesh];
pltfn[u_, t_] :=
Module[{plti, pltinf, plto},
plti = ContourPlot[u[][t, x, y], Element[{x, y}, innerMesh],
AspectRatio -> Automatic, PlotPoints -> All, PlotRange -> {0, c0},
ColorFunction -> "DarkBands"];
pltinf =
ContourPlot[u[][t, x, y], Element[{x, y}, membraneMesh],
AspectRatio -> Automatic, PlotPoints -> All, PlotRange -> {0, c0},
ColorFunction -> "DarkBands"];
plto = ContourPlot[u[][t, x, y], Element[{x, y}, outerMesh],
AspectRatio -> Automatic, PlotPoints -> All, PlotRange -> {0, c0},
ColorFunction -> "DarkBands"];
Show[plto, pltinf, plti]] • This looks like a very promising approach! Well done Tim. – user21 Apr 9 at 9:49
• I did ask not to do this, but you convinced me, so accepting this one ... – Szabolcs Apr 9 at 13:24
• @Szabolcs Thank you for the accept! I was not expecting it. I know that you wanted a membrane free answer, but I thought this was worth sharing versus gathering dust on my hard drive. FEM has its limitations, but it can be very flexible and efficient if you can work within those limitations. – Tim Laska Apr 9 at 19:50
• @user21 Thank you very much! I think it is a good start. Maybe it would be a useful addition to have a region marker preserving ElementMeshJoin to join compatible element types in FEMAddOns? – Tim Laska Apr 9 at 20:06
• This suggestion has been made before - at that time I was hesitant because it very easy to generate faulty meshes with something like this. On the other hand there are certainly cases were such a functionality would make things much easier. Even in the Acoustic Cloak example would benefit from something like that. – user21 Apr 10 at 7:53

Here's a solution using pdetoode to discretize the system in $$x$$ direction. The condition at $$x=1$$ is then straightforwardly introduced in this approach:

{lb = 0, mb = 1, rb = 2, dl = 1, dmem = 2, dr = 3, tmax = 5};

With[{u = u[t, x]}, eq = D[u, t] == # D[D[u, x], x] & /@ {dl, dr};
ic = {u == 2, u == 1} /. t -> 0;
{bcl, bcr} = {{u == 2 /. x -> lb, dl D[u, x] /. x -> mb},
{dr D[u, x] /. x -> mb, u == 1 /. x -> rb }}] ;

points = 25; {gridl, gridr} = Array[# &, points, #] & /@ {{lb, mb}, {mb, rb}};
difforder = 2;

{ptoofuncl, ptoofuncr} = pdetoode[u[t, x], t, #, difforder] & /@ {gridl, gridr};

del = #[[2 ;; -2]] &;
{odel, oder} = MapThread[del@#@#2 &, {{ptoofuncl, ptoofuncr}, eq}];
{odeicl, odeicr} = MapThread[#@#2 &, {{ptoofuncl, ptoofuncr}, ic}];
{odebcl, odebcr} = MapThread[#@#2 &, {{ptoofuncl, ptoofuncr}, {bcl, bcr}}];

linkterm = dmem (ur[t] - ul[t]);
rulel = u -> ul;
ruler = u -> ur;
odebcm = {linkterm == odebcl[] /. rulel, linkterm == odebcr[] /. ruler};

odebc = With[{sf = 1},
Map[sf # + D[#, t] &, Flatten@{odebcl[], odebcr[], odebcm}, {2}]];

sollst = NDSolveValue[{{odel, odeicl} /. rulel, {oder, odeicr} /. ruler,
odebc}, {u /@ gridl // Most, u /@ gridr // Rest, ul, ur}, {t, 0,
tmax}]; // AbsoluteTiming

soll = rebuild[Join[sollst[], {sollst[]}], gridl]
solr = rebuild[Join[{sollst[]}, sollst[]], gridr]
sol = {t, x} \[Function] Piecewise[{{soll[t, x], x < mb}}, solr[t, x]]

Manipulate[Plot[sol[t, x], {x, lb, rb}], {t, 0, tmax}] • To the downvoter, I am interested in what's wrong with my answer, would you please elaborate. I'm not trying to complain here, I just want to improve my answer if possible. – xzczd Apr 11 at 4:02

We can use NDSolve with FEM by changing the variable x->2-x at x>=1 and defining two equations on the same interval (x,0,1), connected for x = 1:

Needs["NDSolveFEM"]; d1 = 1; d2 = 3; dm = 1; reg =
ImplicitRegion[0 <= x <= 1, {x}];
eq = {-d1 Laplacian[u1[t, x], {x}] +
D[u1[t, x], t], -d2 Laplacian[u2[t, x], {x}] + D[u2[t, x], t]};
ic = {u1[0, x] == 2, u2[0, x] == 1};
bc1 = NeumannValue[-dm (u1[t, x] - u2[t, x]), x == 1];
bc2 = NeumannValue[-dm (u2[t, x] - u1[t, x]), x == 1];
bc = DirichletCondition[{u1[t, x] == 2, u2[t, x] == 1}, x == 0];

{U1, U2} =
NDSolveValue[{eq[] == bc1, eq[] == bc2, bc, ic}, {u1, u2},
x \[Element] reg, {t, 0, 2}]


Visualisation

Plot3D[{U1[t, x], U2[t, 2 - x]}, {x, 0, 2}, {t, 0, 2},
AxesLabel -> Automatic]


If the membrane set in an arbitrary point p, then the code should be modified as follows:

Needs["NDSolveFEM"]; d1 = 1; d2 = 3; dm = 1; reg =
ImplicitRegion[0 <= x <= 1, {x}]; p = 1.2; x1 =
x/p; x2 = (2 - x)/(2 - p); k1 = D[x1, x]; k2 = D[x2, x];
eq = {-d1 k1^2 Laplacian[u1[t, x], {x}] +
D[u1[t, x], t], -d2 k2^2 Laplacian[u2[t, x], {x}] +
D[u2[t, x], t]};
ic = {u1[0, x] == 2, u2[0, x] == 1};
bc1 = NeumannValue[-dm k1 (u1[t, x] - u2[t, x]), x == 1];
bc2 = NeumannValue[dm k2 (u2[t, x] - u1[t, x]), x == 1];
bc = DirichletCondition[{u1[t, x] == 2, u2[t, x] == 1}, x == 0];
{U1, U2} =
NDSolveValue[{eq[] == bc1, eq[] == bc2, bc, ic}, {u1, u2},
x \[Element] reg, {t, 0, 2}]


Visualisation

{Table[Plot[{U1[t, x1], U2[t, x2]}, {x, 0, 2}, PlotRange -> All,
PlotLabel -> Row[{"t = ", t}]], {t, 0, 2, .5}],
Plot3D[{U1[t, x1], U2[t, x2]}, {x, 0, 2}, {t, 0, 2},
AxesLabel -> Automatic]} • I came up with something similar before, but this technique doesn't really generalize for other domain shapes. What if the membrane is at x=0.5, not at x=1? I guess one domain can be artificially compressed to make it the same size as the other, and the equations re-scaled accordinagly. But what if we are in 2D and the membrane is circular? – Szabolcs Apr 7 at 17:40
• @Szabolcs There is no general solution for this problem. If you interested in 2D or 3D solution then we could try to make code. But my answer is only for this special 1D case you posted above. – Alex Trounev Apr 7 at 17:55
• For now, I will ask for one generalization only: to put the membrane at an arbitrary position along $[0,2]$ instead of right in the middle. Anything beyond that should go in a separate question. – Szabolcs Apr 8 at 8:02
• @Szabolcs You may want to try DChange. – xzczd Apr 8 at 9:49
• @Szabolcs See update to my answer. – Alex Trounev Apr 8 at 10:58

Version 12.1.1 has a new PDE modeling tutorial on Modeling Mass Transport. The section Interphase Mass Transfer shows an example of how inter phase mass transport can be modeled with NDSolve. The related question of modeling phase change is shown in the section Transient PDEs with Nonlinear Transient Coefficients of the Finite Element Programming tutorial