I have the following system of diffusion equations with these boundary conditions.

Diffusion Equations

Boundary Conditions

The second condition is a Dirichlet condition at x approaching infinity. How do I implement that? I've tried replacing it with 'b' and then applying the limit but that hasn't worked.

eqn1 := D[cA[x, t], t] - dA D[cA[x, t], {x, 2}];
eqn2 := D[cB[x, t], t] - dB D[cB[x, t], {x, 2}];

{ncA, ncB} = NDSolveValue[{
    eqn1 == 
     NeumannValue[kc/dA cA[x, t] - ka/dA cB[x, t], (t > 0) && (x == 0)],
    DirichletCondition[cA[x, t] == cAbulk, t == 0],
    DirichletCondition[cA[x, t] == cAbulk, (t > 0) && (x == \[Infinity])],
    eqn2 == 
     NeumannValue[-(kc/dB) cA[x, t] + ka/dB cB[x, t], (t > 0) && (x == 0)], 
    DirichletCondition[cB[x, t] == 0, t == 0],
    DirichletCondition[cB[x, t] == 0, (t > 0) && (x == \[Infinity])]
    }, {cA, cB}, {x, 0, 1}, {t, 0, 1}];

I'm not sure if I should ask this separately or here but on a related note, I've never seen the Neumann Condition written as an inhomogeneity to the original differential equation like it does in the help file (or in other posts on stack exchange). I've consulted a few books on the subject. Is it mathematically valid or some quirk of Mathematica?

Thank you in advance.

Edit: The code does run when I evaluate it at the limit of NDSolveValue which is x=1. Nevertheless my question is, can I apply it at infinity directly?

  • 2
    $\begingroup$ “my question is, can I apply it at infinity directly?” No, there's no built-in support for b.c. at infinity in NDSolve. As shown in the answered question under the tag boundary-condition-at-infinity, the usual treatments are: 1. use a large enough number to approximate infinity. 2. Use proper coordinate transformation to map the infinite domain to a finite one. $\endgroup$ – xzczd Feb 17 at 4:22
  • $\begingroup$ Thank you. Is there a resource on the way Mathematica uses eqn==NeumannValue? Is it standard practice for PDEs in mathematics? $\endgroup$ – Walser Feb 17 at 4:37
  • 2
    $\begingroup$ The document of NeumannValue is clear enough in my view, just check it by pressing F1. Still, if you have difficulty in understanding its usage, consider avoiding DirichletCondition and NeumannValue and turn to the old good TensorProductGrid method for spatial discretization. (Given that you're dealing with a regular domain, it's a reasonable choice. ) $\endgroup$ – xzczd Feb 17 at 5:00

"Infinite domains" with anisotropic meshing

I will demonstrate an approach that extends the domain by a factor of a thousand with a small increase in computational cost through anisotropic meshing. I will verify the approach using heat transfer in an infinite rod where the analytical solution is known.

Helper functions

I will use some of the following helper functions to build the mesh segments.

(*Import required FEM package*)
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] := 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=(*Quiet@*)
    firstElmHeight[x0, xf, n, r] - fElm, {r, 0.0001, 10000/fElm}, 
    Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := 
 Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
reflectRight[pts_] := 
 With[{rt = ReflectionTransform[{1}, {Last@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

Create anisotropic mesh

Here, I create a graded mesh on the finite domain where the mesh size geometrically grows from the $x=0$ wall. Then, I create an infinite domain where the 1st element size matches the finite domain's last element size and grows geometrically out to infinity ($x=1000$).

(*Create mesh segment with finer mesh at x=0*)
seg = leftSegmentGrowth[1, 100, 1/1000];
Print["Region of interest domain"]
r = pointsToMesh@seg
(*Infinite segment where 1st element matches last element of finite \
seginf = meshGrowthByElm0[1000, 50, Last@seg - seg[[-2]]];
Print["Extended domain"]
Print["Full domain"]
regfull = pointsToMesh@extendMesh[seg, seginf]
(*Extract Coords from region*)
crd = MeshCoordinates[r];
crdinf = MeshCoordinates[regfull];
(*Create element mesh*)
mesh = ToElementMesh[crd];
meshinf = ToElementMesh[crdinf];
Print["Coordinate positions on a log scale"]
ListLogPlot[Transpose@meshinf["Coordinates"], PlotRange -> All]

Computational mesh

Set up and solve the heat transfer problem for "infinite rod."

We can use the HeatTransferPDEComponent to set up our PDE system for heat transfer an infinite rod. Initially, the system starts at a temperature of $\theta=0$. Instantaneously, we increase the wall's temperature to ($\theta(0,0)=1$). To avoid inconsistency in the boundary and initial conditions, we will ramp up the temperature rapidly using a $Tanh$ function. We will set up systems for both the finite and the "infinite" domains for comparison purposes.

(*Use HeatTransferPDEComponent to build heat transfer operator*)
vars = {Θ[t, x], t, {x}};
pars = <|"MassDensity" -> 1, "SpecificHeatCapacity" -> 1, 
   "ThermalConductivity" -> 1|>;
bcs = {DirichletCondition[Θ[t, x] == Tanh[5000 t], 
    x == 0], DirichletCondition[Θ[t, x] == 0, x == 1]};
bcsinf = {DirichletCondition[Θ[t, x] == Tanh[5000 t], 
    x == 0], 
   DirichletCondition[Θ[t, x] == 0, x == 1001]};
ics = Θ[0, x] == 0;
eqn = HeatTransferPDEComponent[vars, pars] == 0;
Tfun = NDSolveValue[{eqn, bcs, ics}, Θ, {t, 0, 20}, 
   x ∈ mesh, StartingStepSize -> 0.000001];
Tfuninf = 
  NDSolveValue[{eqn, bcsinf, ics}, Θ, {t, 0, 20}, 
   x ∈ meshinf, StartingStepSize -> 0.000001];

Comparison of FEM solutions with the analytical solution.

The following equation gives the analytical solution for the infinite rod heat transfer problem:

$$\Theta \left( {t,x} \right) = Erfc\left( {\frac{x}{{\sqrt {4\alpha t} }}} \right);\alpha = \frac{k}{{\rho {{\hat C}_p}}}$$

We can compare the solutions with the following code:

(*Animate solution over a short timescale*)
frames = Table[
   Plot[{Tfun[t, x], Tfuninf[t, x], If[t == 0, 0, Erfc[x/Sqrt[4 t]]]},
     x ∈ r, 
    PlotStyle -> {Thick, 
      Directive[Dashed, Thickness[0.01], Black, Opacity[1]], 
      Directive[Thickness[0.02], Green, Opacity[0.5]]}, 
    PlotLegends -> {"\!\(\*SubscriptBox[\(x\), \(nax\)]\)=1", 
      "\!\(\*SubscriptBox[\(x\), \(max\)]\)=1000", "Analytical"}, 
    PlotRange -> {0, 1.01}, PlotLabel -> t], {t, 0, 0.2, 0.002}];
(*Animate solution over a long timescale*)
frames = Table[
   Plot[{Tfun[t, x], Tfuninf[t, x], If[t == 0, 0, Erfc[x/Sqrt[4 t]]]},
     x ∈ r, 
    PlotStyle -> {Thick, 
      Directive[Dashed, Thickness[0.01], Black, Opacity[1]], 
      Directive[Thickness[0.02], Green, Opacity[0.5]]}, 
    PlotLegends -> {"\!\(\*SubscriptBox[\(x\), \(nax\)]\)=1", 
      "\!\(\*SubscriptBox[\(x\), \(max\)]\)=1000", "Analytical"}, 
    PlotRange -> {0, 1.01}, PlotLabel -> t], {t, 0, 20, 0.2}];

Short timescale

Longtime scale

The finite domain solution starts to break down when the temperature front of the analytical solution reaches the end of the domain ($t=0.08 s$). The infinite domain finite element solution tracks the analytical solution quite well over a long time scale. We achieve this by the addition of only 50 additional elements.

The approach is quite robust. Here, I extended the domain to be 10,000 times larger than the region of interest and simulated out to 10 million seconds without any issue.

10 million seconds simulation

Application to the OP problem

Now that we have verified that the anisotropic meshing approach is quite robust over a huge domain and large timescales, we will apply it to the OP problem.

Following @user21's links and advice in his answer, we will use the MassTransportPDEComponent to recast and solve the OP problem. We will use the same parameters as @user21's as they were not provided in the OP. The workflow follows:

(*Define some parameters*)
dA = 1; dB = 1; ka = 1; kc = 1; cAbulk = 1; xf = 1; large = 1000; \
nelm = 100; ninf = 50;
vars = {{cA[t, x], cB[t, x]}, t, {x}};
pars = <|"DiffusionCoefficient" -> {{dA, 0}, {0, dB}}, 
   "BoundaryCondition1" -> <|
     "AmbientConcentration" -> {kc/ka cB[t, x], ka/kc cA[t, x]}, 
     "MassTransferCoefficient" -> {ka, kc}|>, 
   "BoundaryCondition2" -> <|"MassConcentration" -> {cAbulk, 0}|>|>;
ops = MassTransportPDEComponent[vars, pars];
Γflux = 
  MassTransferValue[x == 0, vars, pars, "BoundaryCondition1"];
Γcond = 
  MassConcentrationCondition[x == xf, vars, pars, 
Γcondinf = 
  MassConcentrationCondition[x == xf + large, vars, pars, 
ics = {cA[0, x] == cAbulk, cB[0, x] == 0};
{cAfun, cBfun} = 
  NDSolveValue[{ops == Γflux, Γcond, 
    ics}, {cA, cB}, {t, 0, 10}, x ∈ mesh];
{cAfuninf, cBfuninf} = 
  NDSolveValue[{ops == Γflux, Γcondinf, 
    ics}, {cA, cB}, {t, 0, 10}, x ∈ meshinf];
frames = Table[
   Plot[{cAfun[t, x], cBfun[t, x], cAfuninf[t, x], cBfuninf[t, x]}, 
    x ∈ mesh, PlotRange -> All, 
    PlotLegends -> {"CA", "CB", "CAinf", "CBinf"}, 
    PlotLabel -> t], {t, 0, 2, 0.01}];

Concentration plots

With the chosen parameters, the finite domain begins to deviate from the "infinite domain" at about $t = 0.2 seconds $.

  • 1
    $\begingroup$ Great answer, as always! If you want to you could add a request for anisotropic meshes to the list under Future enhancements for the finite element method $\endgroup$ – user21 Feb 18 at 5:21
  • $\begingroup$ @user21 thank you very much! I will see what I can pull together. It should be a handy addition to the toolbox. $\endgroup$ – Tim Laska Feb 18 at 15:05
  • $\begingroup$ @Walser, Tim made a proposal to add anisotropic meshing capabilities to the FEM mesh generation here. If you think this is worthwhile, consider giving this proposal an upvote $\endgroup$ – user21 Feb 22 at 10:31

Here is something to get you started.

dA = 1; dB = 1; ka = 1; kc = 1; cAbulk = 1;
eqn1 = D[cA[x, t], t] - 
   Inactive[Div][dA Inactive[Grad][cA[x, t], {x}], {x}];
eqn2 = D[cB[x, t], t] - 
   Inactive[Div][dB Inactive[Grad][cB[x, t], {x}], {x}];

large = 10;
{ncA, ncB} = 
  NDSolveValue[{eqn1 == 
     NeumannValue[kc cA[x, t] - ka cB[x, t], x == 0], 
    cA[x, 0] == cAbulk,
    DirichletCondition[cA[x, t] == cAbulk, x == large], 
    eqn2 == NeumannValue[-kc cA[x, t] + ka cB[x, t], x == 0],
    cB[x, 0] == 0,
    DirichletCondition[cB[x, t] == 0, x == large]
    }, {cA, cB}, {x, 0, large}, {t, 0, 1}];
Plot3D[ncA[x, t], {t, 0, 1}, {x, 0, large}, PlotRange -> All]

enter image description here

There is a Mass Transport tutorial and a related Mass Transport and Boundary Conditions Guide page that seem relevant.

Concerning the NeumannValue, yes that is common in the derivation of the equations. The NeumannValue is essentially a boundary integral that comes up after the application of an integration of parts. This is described briefly in the NeumannValue ref page and more elaborately in the Solving PDE with FEM tutorial. The Finite Element Usage Tips tutorial also has some examples how this relates to other boundary conditions.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.