Message Boards Message Boards


Modeling jump conditions in interphase mass transfer

Posted 3 years ago
19 Replies
29 Total Likes

NOTE: Download Full Article as a Notebook from the Attachment Below

enter image description here

Interphase mass transfer operations such as gas absorption or liquid-liquid extraction pose a modeling challenge because the molar species concentration can jump between 2 states at the interface as shown below (from here). enter image description here I wanted to see if I could create a Finite Element Method (FEM) model of jump conditions in the Wolfram Language. I found the results to be reasonable, aesthetically pleasing, and somewhat mesmerizing. The remainder of this post documents my workflow for those that might be interested. I have attached a notebook to reproduce the results.

Preamble on analytical solutions to PDE's

There seems to be quite a few posts where people are trying to find the analytical solution to a system of PDE's. Generally, closed formed analytical solutions only exist in rare-highly symmetric cases. Let us consider the heat equation below.

$$\frac{\partial T}{\partial t}=\alpha \frac{\partial^2 T}{\partial x^2}$$

For the case of a semi-infinite bar subjected to a unit step change in temperature at $x=0$, Mathematica's DSolve[] handles this readily.

u[x, t] /. 
 First@DSolve[{ D[u[x, t], t] == alpha D[u[x, t], {x, 2}], 
    u[x, 0] == 0, u[0, t] == UnitStep[t]}, u[x, t], {x, t}, 
   Assumptions -> alpha > 0 && 0 < x]

So far so good. Now, let us break symmetry by making it a finite bar of length $l$ (See Documentation).

heqn = D[u[x, t], t] == alpha D[u[x, t], {x, 2}];
bc = {u[0, t] == 1, u[l, t] == 0};
ic = u[x, 0] == 0;
sol = DSolve[{heqn, bc, ic}, u[x, t], {x, t}]
(*{{u[x, t] -> 1 - x/l - (2*Inactive[Sum][Sin[(Pi*x*K[1])/l]/(E^((alpha*Pi^2*t*K[1]^2)/l^2)*
         K[1]), {K[1], 1, Infinity}])/Pi}}*)

This little change going from a semi-infinite to finite domain has turned the solution into an unwieldy infinite sum. We should expect that it will only go down hill from here if we add additional complexity to the equation or system of equations. My advice is to abandon the search for an analytical solution quickly because it will likely take great effort and will be unlikely to yield a result. Instead, focus efforts on more productive avenues such as dimensional analysis and numerical solutions.


"All models are wrong, some are useful." -- George E. P. Box "However, many systems are highly complex, so that valid mathematical models of them are themselves complex, precluding any possibility of an analytical solution. In this case, the model must be studied by means of simulation, i.e. , numerically exercising the model for the inputs in question to see how they affect the output measures of performance." -- Dr. Averill Law, Simulation Modeling and Analysis

I find the quotes above help me overcome inertia when starting a modeling and simulation project. Create your wrong model. Calibrate how wrong it is versus a known standard. If it is not too bad, put the model through its paces.

One thing that I appreciate about the Wolfram Language is that I can document a modeling workflow development process from beginning to end in a single notebook. The typical model workflow development process includes:

  • A sketch of the system of interest.
  • Equations.
    • Initial development.
    • Simplification.
    • Non-dimensionalization for better scaling and reducing parameter space.
  • Mathematica implementation.
    • Mesh
      • Boundaries
      • Refinement
    • NDSolve set-up
    • Post-process results
  • Verification/Validation

Mathematica notebooks tend to age well. I routinely resurrect notebooks that are over a decade old and they generally still work.


I did a quick Google search on absorption and came across this figure describing gas absorption in an open access article by Danish et al.

enter image description here

This image looked very similar to an image that I produced in a related post to the Wolfram community regarding porous media energy transport.

enter image description here

The systems look so similar, that we ought to be able to reuse much of the modeling. An area of concern would be for gas absorption where the ratio of the gas diffusion coefficient to liquid diffusion coefficient can exceed 4 orders of magnitude. Such differences often can cause instability in numerical approaches.


System description

For clarity, I always like to begin with a system description. Typically, absorption processes utilize gravity to create a thin liquid film to contact the gas. To reuse the modeling that we did for porous media, we will assume that gravity is in the positive $x$ direction leading us to the image below.

enter image description here

We will assume that the liquid film is a uniform thickness and is laminar (note that for gas liquid contact the liquid velocity is fastest at the interface leading to the parabolic profile shown). We will assume that the gas has a uniform velocity. Further, we will assume that the incoming concentrations of the absorbing species are zero and we will impose a concentration of $C=C_0$ at the lower boundary.

The basic dimensions of the box are shown below. For simplicity, we will make the height and length unit dimensions and set $R$ to be $\frac{1}{2}$.

enter image description here

Balance equations

Dilute species balance

For the purposes of this exercise, we will consider the system to be dilute such that diffusion does not affect the overall flow velocities. Within a single phase, the molar balance of concentration is given by equation (1). We will assume steady-state operation with no reactions so that we can eliminate the red terms.

$${\color{Red}{\frac{{\partial {C}}}{{\partial t}}}} + \mathbf{v} \cdot \nabla C - \nabla \cdot \mathcal{D} \nabla C - {\color{Red}{r^{'''}}} = 0 \qquad (1)$$

Species balance in each phase

For convenience, I will denote the phases by a subscript G and L for gas and liquid with the understanding that these equations could also apply to a liquid-liquid extraction problem. This leads to the following concentration balance equations for the liquid and gas phases.

$$\begin{gathered} \begin{matrix} \mathbf{v}_L \cdot \nabla C_L + \nabla \cdot \left(-\mathcal{D}_L \nabla C_L\right) = 0 & x,y\in \Omega_L & (2*) \\ \mathbf{v}_G \cdot \nabla C_G + \nabla \cdot \left(-\mathcal{D}_G \nabla C_G\right) = 0 & x,y\in \Omega_G & (3*) \\ \end{matrix} \end{gathered}$$

Or in Laplacian form

$$\begin{gathered} \begin{matrix} \mathbf{v}_L \cdot \nabla C_L -\mathcal{D}_L \nabla^2 C_L = 0 & x,y\in \Omega_L & (2*) \\ \mathbf{v}_G \cdot \nabla C_G -\mathcal{D}_G \nabla^2 C_G = 0 & x,y\in \Omega_G & (3*) \\ \end{matrix} \end{gathered}$$

Creating a No-Flux Boundary Condition at the Interface

To prevent the gas species diffusing into the liquid layer and vice versa, I will set the velocities to zero and the diffusion coefficients to a very small value in the other phase. From a visualization standpoint, it will appear that the gas species has diffused into the liquid layer and vice versa, but the flux is effectively zero. To clean up the visualization, we will define plot ranges by gas, interphase, and liquid regions.

Species balance including a thin interphase region

We will define a thin Interphase region between the 2 phases that will allow us to couple the phases in the interphase region via a source term creating the jump discontinuity in concentration as shown in the figure below.

enter image description here

We will modify (2*) and (3*) with the coupling source term as shown below.

$$\begin{gathered} \begin{matrix} \mathbf{v}_L \cdot \nabla C_L -\mathcal{D}_L \nabla^2 C_L - \sigma\left(\Omega \right )k\left(K C_G-C_L \right ) = 0 & x,y\in \Omega & (2) \\ \mathbf{v}_G \cdot \nabla C_G -\mathcal{D}_G \nabla^2 C_G + \sigma\left(\Omega \right )k\left(K C_G-C_L \right ) = 0 & x,y\in \Omega & (3) \\ \end{matrix} \end{gathered}$$

Where $K$ is a vapor-liquid equilibrium constant, $k$ is in interphase mass transfer coefficient (we will make this large because we want a fast approach to equilibrium), and $\sigma$ is a switch that turns on (=1) in the interface region and 0 otherwise.

Dimensional analysis

We will multiply equations (2) and (3) by $\frac{{R^2}}{C_0 \mathcal{D}_G}$ to obtain their non-dimensionalized forms (4) and (5).

$$\begin{gathered} \begin{matrix} C_0\left (\frac{\mathbf{v}_L}{R} \cdot \nabla^* C_{L}^{*} -\frac{\mathcal{D}_L}{R^2} \nabla^{*2} C_{L}^{*} - \sigma\left(\Omega \right )k\left(K C_{G}^{*}-C_{L}^{*} \right ) \right ) = 0\left\| {\frac{{R^2}}{C_0 \mathcal{D}_G}} \right. \\ C_0\left (\frac{\mathbf{v}_G}{R} \cdot \nabla^* C_{G}^{*} -\frac{\mathcal{D}_G}{R^2} \nabla^{*2} C_{G}^{*} + \sigma\left(\Omega \right )k\left(K C_{G}^{*}-C_{L}^{*} \right ) \right ) = 0\left\| {\frac{{R^2}}{C_0 \mathcal{D}_G}} \right. \\ \end{matrix} \end{gathered}$$

$$\begin{gathered} \begin{matrix} \frac{\mathcal{D}_L}{\mathcal{D}_G} \frac{R\mathbf{v}_L}{\mathcal{D}_L} \cdot \nabla^* C_{L}^{*} -\delta \nabla^{*2} C_{L}^{*} - \sigma\left(\Omega \right )\kappa\left(K C_{G}^{*}-C_{L}^{*} \right ) = 0 \\ \frac{R\mathbf{v}_G}{\mathcal{D}_G} \cdot \nabla^* C_{G}^{*} - \nabla^{*2} C_{G}^{*} + \sigma\left(\Omega \right )\kappa\left(K C_{G}^{*}-C_{L}^{*} \right ) = 0\\ \end{matrix} \end{gathered}$$

$$\begin{gathered} \begin{matrix} \delta{Pe}_L\mathbf{v}_{L}^* \cdot \nabla^* C_{L}^{*} -\delta \nabla^{*2} C_{L}^{*} - \sigma\left(\Omega \right )\kappa\left(K C_{G}^{*}-C_{L}^{*} \right ) = 0 & (4) \\ {Pe}_G\mathbf{v}_{G}^* \cdot \nabla^* C_{G}^{*} - \nabla^{*2} C_{G}^{*} + \sigma\left(\Omega \right )\kappa\left(K C_{G}^{*}-C_{L}^{*} \right ) = 0 &(5)\\ \end{matrix} \end{gathered}$$


$$\delta=\frac{\mathcal{D}_L}{\mathcal{D}_G}$$ $$Pe_L=\frac{R\mathbf{v}_L}{\mathcal{D}_L}$$ $$Pe_G=\frac{R\mathbf{v}_G}{\mathcal{D}_G}$$

With a good dimensionless model in place, we can start with our Wolfram Language implementation.

Wolfram Language Implementation

Mesh creation

We start by loading the FEM package.


When I started this effort, I considered co-current flow only. I realized that converting the model to counter-current flow was a simple matter of changing boundary markers. I wrapped the process up in a module that returns an association whose parameters can be used to set up an NDSolve solution. In the counter-current case, I changed the bottom boundary to a wall and the gas inlet concentration to 1.

makeMesh[h_, l_, rat_, gf_, cf_] := 
 Module[{bR, tp, bt, lf, rt, th, interfacel, interfaceg, buf, bnds, 
   rgs, crds, lelms, boundaryMarker, bcEle, bmsh, liquidCenter, 
   liquidReg, interfaceCenter, interfaceReg, gasCenter, gasReg, 
   meshRegs, msh, mDic},
  (* Domain Dimensions *)
  bR = rat h;
  tp = bR;
  bt = bR - h;
  lf = 0;
  rt = l;
  th = h/gf;
  interfacel = 0;
  interfaceg = interfacel - th;
  buf = 2.5 th;

  (* Use associations for clearer assignment later *)
  bnds = <|liquidinlet -> 1, gasinlet -> 2, bottom -> 3|>;
  rgs = <|gas -> 10, liquid -> 20, interface -> 15|>;

  (* Meshing Definitions *)
  (* Coordinates *)
  crds = {{lf, bt}(*1*), {rt, bt}(*2*), {rt, tp}(*3*), {lf, 
     tp}(*4*), {lf, interfacel}(*5*), {rt, interfacel}(*6*), {lf, 
     interfaceg}(*7*), {rt, interfaceg}(*8*)};

  (* Edges *)
  lelms = {{1, 7}, {7, 5}, {5, 4}, {1, 2},
           {2, 8}, {8, 6}, {6, 3}, {3, 4}, 
           {5, 6}, {7, 8}};

  (* Conditional Boundary Markers depending on configuration *)
  boundaryMarker := {bnds[gasinlet], bnds[liquidinlet], 
     bnds[liquidinlet], bnds[bottom], 4, 4, 4, 4, 4, 4} /; cf == "Co";
  boundaryMarker := {4, 4, bnds[liquidinlet], bnds[bottom], 
     bnds[gasinlet], 4, 4, 4, 4, 4} /; cf == "Counter";

  (* Create Boundary Mesh *)
  bcEle = {LineElement[lelms, boundaryMarker]};
  bmsh = ToBoundaryMesh["Coordinates" -> crds, 
    "BoundaryElements" -> bcEle];

  (* 2D Regions *)
  (* Identify Center Points of Different Material Regions *)
  liquidCenter = {(lf + rt)/2, (tp + interfacel)/2};
  liquidReg = {liquidCenter, rgs[liquid], 0.0005};
  interfaceCenter = {(lf + rt)/2, (interfacel + interfaceg)/2};
  interfaceReg = {interfaceCenter, rgs[interface], 0.5*0.000005};
  gasCenter = {(lf + rt)/2, (bt + interfaceg)/2};
  gasReg = {gasCenter, rgs[gas], 0.0005};
  meshRegs = {liquidReg, interfaceReg, gasReg};

  msh = ToElementMesh[bmsh, "RegionMarker" -> meshRegs,
    MeshRefinementFunction -> Function[{vertices, area},
      Block[{x, y},
       {x, y} = Mean[vertices];
        (y > interfaceCenter[[2]] - buf &&
            y < interfaceCenter[[2]] + buf)  ||
         (y < bt + 1.5 buf &&
            x < lf + 1.5 buf)
        , area > 0.0000125, area > 0.01

  mDic = <|
    height -> h,
    length -> l,
    ratio -> rat,
    gapfactor -> gf,
    r -> bR,
    top -> tp,
    bot -> bt,
    left -> lf,
    right -> rt,
    intl -> interfacel,
    intg -> interfaceg,
    intcx -> interfaceCenter[[1]],
    intcy -> interfaceCenter[[2]],
    buffer -> buf,
    mesh -> msh,
    bmesh -> bmsh,
    bounds -> bnds,
    regs -> rgs,
    cfg -> cf

Options[meshfn] = {height -> 1, length -> 1, ratio -> 0.5, 
   gapfactor -> 100, config -> "Co"};
meshfn[OptionsPattern[]] := 
 makeMesh[OptionValue[height], OptionValue[length],
  OptionValue[ratio], OptionValue[gapfactor], OptionValue[config]]

We can create a mesh instance of a co-current flow case by invoking the meshfn[]. I will color the liquid inlet $\color{Green}{Green}$, the gas inlet $\color{Red}{Red}$ and the bottom boundary $\color{Orange}{Orange}$ (the rest of the boundaries are default).

mDicCo = meshfn[config -> "Co"];
 "Wireframe"["MeshElementMarkerStyle" -> Blue, 
  "MeshElementStyle" -> {Green, Red, Orange, Black}, 
  ImageSize -> Large]]

enter image description here

By setting the optional config parameter to "Counter", we can easily generate a counter-current case as shown below (note how the gas inlet shifted to the right side.

mDic = meshfn[config -> "Counter"];
 "Wireframe"["MeshElementMarkerStyle" -> Blue, 
  "MeshElementStyle" -> {Green, Red, Orange, Black}, 
  ImageSize -> Large]]

enter image description here

For the co-current case, the bottom wall and the gas inlet have inconsistent Dirichlet conditions. To reduce the effect, I refined the mesh in the lower left corner as shown below.

enter image description here

I also meshed the interface region finely.

enter image description here

Sometimes it can get confusing to setup alternative boundaries. To visualize the coordinate IDs, you could use something like:

With[{pts = mDic[bmesh]["Coordinates"]}, 
 Graphics[{Opacity[1], Black, 
    Text[Style[ToString[#], Background -> White, 12], #] & /@ 

enter image description here

Solving and Visualization

I have created a module that will solve and visualize depending on the mesh type (co-flow or counter-flow). Hopefully, it is well enough commented that further discussion is not needed.

model[md_, kequil_, d_, pel_, peg_, title_] := 
  Module[{n, pecletgas, por, vl, vg, fac, facg, coefl, coefg, 
    dcliquidinletliquid, dcliquidinletgas, dcinletgas, 
    dcgasinletliquid, dcgasinletgas, dcbottomliquid, dcbottomgas, 
    eqnliquid, eqngas, eqns, ifun, pltl, pltint, pltg, pltarr, sz, 
    grid, lf, rt, tp, bt, interfaceg, interfacel, interfaceCenterY, 
    plrng, arrequil, arrdiff, arrgas, 
   (*localize Mesh Dict Values*)lf = md[left];
   rt = md[right];
   tp = md[top];
   bt = md[bot];
   interfaceg = md[intg];
   interfacel = md[intl];
   interfaceCenterY = md[intcy];
   (*Must swtich gas flow direction for counter-flow*)
   pecletgas = If[md[cfg] == "Co", peg, -peg];
   (*Dimensionless Mass Transfer Coefficient in Interphase Region*)
   n = 10000;
   (*"Porosity" to weight concentration in interphase*)
   por[y_, intg_, intl_] := (y - intg)/(intl - intg);
   (*Region Dependent Properties with Piecewise \
Functions*)(*velocity*)(*Liquid parabolic profile*)
   vl = Evaluate[
     Piecewise[{{{pel d (1 - (y/md[r])^2), 0}, 
        ElementMarker == md[regs][liquid]}, {{pecletgas, 0}, 
        ElementMarker == md[regs][gas]}, {{0, 0}, True}}]];
   (*Gas Uniform Velocity*)
   vg = Evaluate[
     Piecewise[{{{pecletgas, 0}, 
        ElementMarker == md[regs][gas]}, {{pel d (1 - (y/md[r])^2), 
         0}, ElementMarker == md[regs][liquid]}, {{0, 0}, True}}]];
   (*fac switches on mass transfer coefficient in interphase*)
   fac = Evaluate[If[ElementMarker == md[regs][interface], n, 0]];
   (*diffusion coefficients*)(*Liquid*)
   coefl = Evaluate[
     Piecewise[{{d, ElementMarker == md[regs][liquid]}, {1, 
        ElementMarker == md[regs][interface]}, {d/1000000, 
        True} (*Effectively No Flux at Interface*)}]];
   (*Gas*)coefg = 
     Piecewise[{{1, ElementMarker == md[regs][gas]}, {1, 
        ElementMarker == md[regs][interface]}, {d/1000000, 
        True} (*Effectively No Flux at Interface*)}]];
   (*Dirichlet Conditions for Liquid at Inlets*)
   dcliquidinletliquid = 
    DirichletCondition[cl[x, y] == 0, 
     ElementMarker == md[bounds][liquidinlet]];
   dcliquidinletgas = 
    DirichletCondition[cg[x, y] == 0, 
     ElementMarker == md[bounds][liquidinlet]];
   dcgasinletliquid = 
    DirichletCondition[cl[x, y] == 0, 
     ElementMarker == md[bounds][gasinlet]];
   (*Conditional BCs for gas dependent on configuration*)
   dcgasinletgas := 
    DirichletCondition[cg[x, y] == 0, 
      ElementMarker == md[bounds][gasinlet]] /; md[cfg] == "Co";
   dcgasinletgas := 
    DirichletCondition[cg[x, y] == 1, 
      ElementMarker == md[bounds][gasinlet]] /; md[cfg] == "Counter";
   (*Dirichlet Conditions for the Bottom Wall*)
   dcbottomliquid = 
    DirichletCondition[cl[x, y] == 0, 
     ElementMarker == md[bounds][bottom]];
   dcbottomgas = 
    DirichletCondition[cg[x, y] == 1, 
     ElementMarker == md[bounds][bottom]];
   (*Balance Equations for Gas and Liquid Concentrations*)
   eqnliquid = 
    vl.Inactive[Grad][cl[x, y], {x, y}] - 
      coefl Inactive[Laplacian][cl[x, y], {x, y}] - 
      fac (kequil cg[x, y] - cl[x, y]) == 0;
   eqngas = 
    vg.Inactive[Grad][cg[x, y], {x, y}] - 
      coefg Inactive[Laplacian][cg[x, y], {x, y}] + 
      fac (kequil cg[x, y] - cl[x, y]) == 0;
   (*Equations to be solved depending on configuration*)
   eqns := {eqnliquid, eqngas, dcliquidinletliquid, dcliquidinletgas, 
      dcgasinletliquid, dcgasinletgas, dcbottomliquid, dcbottomgas} /;
      md[cfg] == "Co";
   eqns := {eqnliquid, eqngas, dcliquidinletliquid, dcliquidinletgas, 
      dcgasinletliquid, dcgasinletgas} /; md[cfg] == "Counter";
   (*Solve the PDE*)
   ifun = NDSolveValue[eqns, {cl, cg}, {x, y} \[Element] md[mesh]];
   (*Visualizations*)(*Create Arrows to represent magnitude of \
dimensionless groups*)(*Equilibrium Arrow*)
   arrequil = {CapForm["Square"], Red, Arrowheads[0.03], 
     Arrow[Tube[{{1 - 0.0125, 0.025, 1}, {1 - 0.0125, 0.025, kequil}},
        0.005], -0.03]};
   (*Diffusion Arrow*)
   arrdiff = {Darker[Green, 1/2], 
     Arrowheads[0.03, Appearance -> "Flat"], 
     Arrow[Tube[{{-0.025, 0.0, 0.0 .025}, {-0.025, 
         0.5 (1 + Log10[d]/4), 0.025}}, 0.005], -0.03]};
   (*Liquid Peclet Arrow*)
   arrliq = {Blue, Dashed, Arrowheads[1.5 0.03], 
     Arrow[Tube[{{0.0, mDic[top] + 0.025, 0.035}, {pel/50, 
         mDic[top] + 0.025, 0.035}}, 1.5 0.005], -0.03 1.5]};
   (*Conditional Gas Peclet Arrow*)
   arrgas := {Black, Dashed, Arrowheads[1.5 0.03], 
      Arrow[Tube[{{0.0, mDic[bot], 1.035}, {peg/50, mDic[bot], 
          1.035}}, 1.5 0.005], -0.03 1.5]} /; md[cfg] == "Co";
   arrgas := {Black, Dashed, Arrowheads[1.5 0.03], 
      Arrow[Tube[{{mDic[right], mDic[bot], 
          1.035}, {mDic[right] - peg/50, mDic[bot], 1.035}}, 
        1.5 0.005], -0.03 1.5]} /; md[cfg] == "Counter";
   (*Set up plots*)(*Common plot options*)
   plrng = {{lf, rt}, {bt, tp}, {0, 1}};
   SetOptions[Plot3D, PlotRange -> plrng, PlotPoints -> {200, 200}, 
    ColorFunction -> 
     Function[{x, y, z}, Directive[ColorData["DarkBands"][z]]], 
    ColorFunctionScaling -> False, MeshFunctions -> {#3 &}, 
    Mesh -> 18, AxesLabel -> Automatic, ImageSize -> Large];
   (*Liquid Plot*)
   pltl = Plot3D[ifun[[1]][x, y], {x, lf, rt}, {y, interfacel, tp}, 
     MeshStyle -> {Black, Thick}];
   (*Interface region Plot*)
   pltint = 
    Plot3D[ifun[[2]][x, y] (1 - por[y, interfaceg, interfacel]) + 
      por[y, interfaceg, interfacel] ifun[[1]][x, y], {x, lf, rt}, {y,
       interfaceg, interfacel}, 
     MeshStyle -> {DotDashed, Black, Thick}];
   (*Gas Plot*)
   pltg = Plot3D[ifun[[2]][x, y], {x, lf, rt}, {y, bt, interfaceg}, 
     MeshStyle -> {Dashed, Black, Thick}];
   (*Grid Plot*)sz = 300;
   grid = 
    Grid[{{Show[{pltl, pltint, pltg}, 
        ViewProjection -> "Orthographic", ViewPoint -> Front, 
        ImageSize -> sz, Background -> RGBColor[0.84`, 0.92`, 1.`], 
        Boxed -> False], 
       Show[{pltl, pltint, pltg}, ViewProjection -> "Orthographic", 
        ViewPoint -> Left, ImageSize -> sz, 
        Background -> RGBColor[0.84`, 0.92`, 1.`], 
        Boxed -> False]}, {Show[{pltl, pltint, pltg}, 
        ViewProjection -> "Orthographic", ViewPoint -> Top, 
        ImageSize -> sz, Background -> RGBColor[0.84`, 0.92`, 1.`], 
        Boxed -> False], 
       Show[{pltl, pltint, pltg}, ViewProjection -> "Perspective", 
        ViewPoint -> {Above, Left, Back}, ImageSize -> sz, 
        Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False]}}, 
     Dividers -> Center];
   (*Reset Plot Options to Default*)
   SetOptions[Plot3D, PlotStyle -> Automatic];
   pltarr = 
    Grid[{{Text[Style[title, Blue, Italic, 24]]}, {Style[
         "\!\(\*SubscriptBox[\(K\), \(C\)]\)=``, \[Delta]=``, \
\!\(\*SubscriptBox[\(Pe\), \(L\)]\)=``, and \
\!\(\*SubscriptBox[\(Pe\), \(G\)]\)=``", 
         NumberForm[kequil, {3, 2}, NumberPadding -> {" ", "0"}], 
         NumberForm[d, {5, 4}, NumberPadding -> {" ", "0"}], 
         NumberForm[pel, {2, 1}, NumberPadding -> {" ", "0"}], 
         NumberForm[peg, {2, 1}, NumberPadding -> {" ", "0"}]], 
        18]}, {Show[{pltl, pltint, pltg, 
         Graphics3D[{arrequil, arrdiff, arrliq, arrgas}](*,arrequil,
         arrdiff,arrliq,arrgas*)}, ViewProjection -> "Perspective", 
        ViewPoint -> {Above, Left, Back}, ImageSize -> 640, 
        Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False, 
        PlotRange -> {{md[left] - 0.05, md[right]}, {md[bot], 
           md[top] + 0.05}, {0, 1 + 0.1}}]}}];
   (*Return values*){ifun, {pltl, pltint, pltg}, pltarr, grid}];

Options[modelfn] = {md -> mDic, k -> 0.5, dratio -> 1, pel -> 50, 
   peg -> 50, title -> "Test"};
modelfn[OptionsPattern[]] := 
 model[OptionValue[md], OptionValue[k], OptionValue[dratio], 
  OptionValue[pel], OptionValue[peg], OptionValue[title]]

Testing of Meshing and Solving Modules

Now, that we wrapped our meshing and solving work flow into modules, I will demonstrate how to create an instance of a simulation.

(* Create a Co-Flow Mesh *)
mDic = meshfn[config -> "Co"];
(* Simulate and return results *)
res = modelfn[md -> mDic, k -> 0.5, dratio -> 0.1, pel -> 10, 
   peg -> 5, title -> "Co-Flow"];

To visualize a 3D plot with arrows representing the magnitude of dimensionless parameters, we access the third part of the results list.


enter image description here

The solid lines, dashed lines, and dashed-dotted lines represent contours of species concentration in the liquid, gas, and interphase regions, respectively. The $\color{Red}{Red}$ arrow is proportional to (1-K), the $\color{Green}{Green}$ arrow is proportional to the log of the diffusion ratio $\delta$, the $\color{Blue}{Blue}$ arrow is proportional to $Pe_L$, and the $\color{Black}{Black}$ arrow is proportional to $Pe_G$. Multiple views are contained in part 4 of the results list.


enter image description here

Validation (Comparison to another code)

Before continuing, it is always good practice to validate your model versus experiment or at least another code. The other code supports a partition conditions for the concentration jump so that I do not need to create an interface layer. The results are shown below:

enter image description here

The contour plots look very similar to the image in the lower left corner of the grid plot. To be more quantitative, I have highlighted contours at approximately y=-0.15 and y=0.05 in the gas and liquid layers at x=1 corresponding to concentrations of 0.68 and 0.28, respectively. The first part of the results list returns an interpolation function of the liquid and gas species. We can see that we are within a percent of the other code, which is reasonable given that the interface layer is about 1% of the domain. This check gives me good confidence that my model is not too wrong and that I can start to make it useful (i.e., exercising the model by changing parameters).

res[[1]][[2]][1, -0.15] (*0.6769985984321076`*)
res[[1]][[1]][1, 0.05] (* 0.27374616012596314`*)

Generating Animations

I like to animate. For me, animations are the best way to demonstrate how a system evolves as a function of time or parameter changes. We can export an animated gif file to study the effects of dimensionless parameter changes for both flow configurations as shown in the following code. It will take about 30 minutes per animation and about 5 GB of RAM. Undoubtedly, this code could be optimized for speed and memory usage, but you still can create a dozen animations while you sleep.


f = ((#1 - #2)/(#3 - #2)) &; (* Scale for progress bar *)

mDic = meshfn[config -> "Counter"]; (* Create Mesh Instance *)

   modelfn[md -> mDic, k -> kc, dratio -> 1, pel -> 0, peg -> 0, 
     title -> "Counter-Flow"][[3]], {kc, 1, 0.01, -0.01}
    {{"Total progress:", 
      Dynamic[f[kc, 1, 
        0.01, -0.01]]]}, {"\!\(\*SubscriptBox[\(K\), \(C\)]\)=", \
 "AnimationRepetitions" -> \[Infinity]]


I combined the co- (left) and counter-current (right) gif animations for several cases below. Péclet numbers approaching 100 start to look uninteresting visually (all the action is very close to the interface). This should inform the user that perhaps another model is in order with new assumptions to study the small-scale behavior near the interface.

Changing the Equilibrium Constant @ No Flow

enter image description here

As the equilibrium constant, $K$, reduces, the jump condition increases.

Changing the Diffusion Ratio $\delta$ @ No Flow

enter image description here

As the liquid-gas diffusion ratio, $\delta$, decreases, the concentration in the gas layer increases. We also see that the solution does not change much for $\delta<0.01$.

Changing $Pe_L$ @ No Gas Flow

enter image description here

As $Pe_L$ increases, the concentration gradient increases at the interface.

Changing $Pe_G$ @ No Liquid Flow

enter image description here

As $Pe_G$ increases, we see the concentration in the liquid layer decrease for co-flow and increase for counter-current flow. This should make sense since the inlet concentration for co-flow is 0 and 1 for counter-current flow.

Changing the Diffusion Ratio $\delta$ @ Middle Conditions

enter image description here

Again, we do not see much change for $\delta<0.01$. One may have noticed that the concentration in the liquid layer goes up as the diffusion coefficient ratio goes down, which may, at first, seem counterintuitive. The reason for this behavior is that the dimensionless velocity in the liquid layer depends on both $Pe_L$ and $\delta$ so it decreases with decreasing $\delta$.


  • Constructed an FEM model in the Wolfram Language to study concentration jump conditions in interphase mass transfer.
  • Results compare favorably to another code designed to handle jump conditions.
  • Showed several examples of the effect of dimensionless parameter changes on two model flow configurations.
  • Notebook provided.
POSTED BY: Tim Laska
19 Replies
Posted 3 years ago

Nice to see a simulation done by MMA! I always wants to learn fem in MMA. Sadly there wasn't many detailed tutorials around.

This may not be on topic but would you mind telling did you model the same problem in comsol? I'm also a comsol user myself so I recognize the icon immediately.

POSTED BY: Ox Clouding

Yes, I did model in both COMSOL and MMA. I generally don't mention other codes because I don't want to appear as though I am promoting possible competition. I am, however, a firm believer in validation versus other standards whenever possible.

Thank you for the comment!

POSTED BY: Tim Laska
Posted 3 years ago

I see. From your experience, what advantages MMA have over other codes like comsol? I always wonder if MMA solvers are superior in terms of speed or accuracy.

POSTED BY: Ox Clouding

To be honest, I am rather new at the FEM in the Wolfram Language, so I have not optimized my approaches or explored parallelizing. For these small 2D models that I have been creating to demonstrate principles, the solve times are short and the results are similar (let's presume that they are accurate) for both solvers. I really have not had to tune the solvers to match results. To see separation in performance (speed, memory, accuracy, robustness, etc.), I would have to start looking at bigger models.

I believe that one needs many tools in the toolbox and that you need to know the strengths and weaknesses of the tools and play to those strengths and avoid those weaknesses. It is likely that one will need to transition between multiple tools as weaknesses are discovered. Since MMA can do algebra, it is ideal to develop equations and map them between the different solver implementations. My current view (subject to change) is that I would use MMA or COMSOL to determine the dominating physics on a geometrically simple model and apply those simplified physics to a geometrically complex model (i.e., large 3D CAD assembly) in another code.

This is actually a complex question with many dimensions, but I hope this helps.

POSTED BY: Tim Laska
Posted 3 years ago


The strength of MMA is its symbolic abilities. But to do fem modeling in MMA is still a hassle - not because its solvers, but the lack of features compared with other FEM softwares:

  1. Ability to show solving progress, i.e. current time, error plots and plot of latest step solution.

  2. Preset equation models. While heat or diffusion equations are easy to write, models like solid mechanics or N-S equations will take more time and prone to bugs.

  3. Somehow exporting animation is much slower on my computer than comsol.

I wish to see MMA improve in these areas in the future.

POSTED BY: Ox Clouding

I agree with your assessment and the assessment in the Mathematica StackExchange link provided by Kuba. I also agree with Alexi’s comment that it is an And proposition versus an Either Or proposition. For me personally, it is also an And plus a bunch of other tools proposition. I generally have more success integrating a bunch of best in class tools to create a simulation workflow than I do using an integrated tool that does everything. Much of what I do is the simulation of highly non-Newtonian fluid flow. My preference is to use Mathematica and COMSOL to study smaller scale phenomena and other CFD packages to study industrial sized models. If you can reduce the complexity of the physics for the large-scale model with Mathematica or COMSOL, you can dramatically increase your probability of a successful simulation.

Generally, I use COMSOL and Mathematica in complementary fashions. As time marches on, there is a blurring of some of the lines between the tools. For example, Mathematica is expanding their PDE capability, while COMSOL is expanding the portability capability with the COMSOL compiler. The computable document format and Wolfram Cloud did set Mathematica apart with its ability to deploy to a broad audience.

COMSOL is definitely a more mature physics modeling platform. As such, it has some convenience features that I will add in addition to yours.

  • Material library
  • Boundary layer meshing for CFD/non-conformal meshes/multiple reference frames
  • Solving steps
  • An application library of over a thousand worked examples
    • Validation and verification examples

Generally, when I am approached to do new modeling project, I state that I need 3 classes of inputs to begin. Namely, I need a description of the geometry, material properties, and operating conditions. Having a material library allows one to help specify the required inputs or potentially use a proxy material if the client cannot provide the inputs.

For highly shear thinning fluids, it is necessary to create a boundary layer mesh on surfaces otherwise the shear stresses will be overestimated. It is unclear to me that this capability exists in the Wolfram language or if the FEM implementation maintains accuracy with high aspect ratio tetrahedra.

Although COMSOL touts its fully coupled multiphysics capability, it is usually much faster and more robust if you can decouple the physics into a series of solution steps. For example, if you want to study heat transfer in a mixing vessel, you could naïvely conduct a fully coupled transient simulation starting from rest. However, experience shows that you can speed up the process over 100x by first solving the flow field with the steady-state frozen rotor approximation without heat transfer, then fix the flow field and solve the transient heat equation only in the next step and get about 90 to 95% of the correct answer. COMSOL makes it relatively easy to create this sequence of study steps (toggling checkboxes and pointing to a solution to initialize from). It should be noted that most solvers can use this strategy to reduce simulation time. It is not clear to me that this capability exists in the Wolfram language in a straightforward way, but if CFD capability is in the development plan, you need the capability to break the simulation and into a series of steps otherwise simulation times will be slow.

To be an effective modeler, it is likely that you will need to create a workflow consisting of multiple steps with different solver settings. The options space is combinatorially complex. Starting from a good initial guess that you can modify to suit your purpose is much simpler than trying to discover the optimal path from the ground up. This is where an extensive application library of solving workflows can be very beneficial. As the Mathematica modeling capability matures, perhaps they could use their curating ability to create a model lookup function similar to their formula lookup function. For example, ModelLookup[“non-isothermal flow”] would return several types of nonisothermal flow workflows (model data objects similar to formula data objects) to choose from.


What I like to do is integrate multiple tools together into a productive workflow. What is usually more interesting to me is not how the tools compete, but rather how I can get them to cooperate to produce something they cannot do alone.

Flight Simulators.

One of the most successful applications of simulation has been that of a flight simulator, which have been commercially available since 1929 (Link Trainer). The main purpose of a flight simulator is to train the neural nets of the prospective pilot to the real time (i.e. interactive) dynamics of multiple controls in multiple environments. We probably can agree that a flight simulator that allowed the pilot to select the aircraft, one throttle, aileron, rudder, elevator, and flap setting and returned a picture of the aircraft flying into a mountain is essentially useless. Yet that is precisely what is being asked when a client asks for the one-off simulation.

So now, we might ask is it possible to make an interactive tool for the “pilots” (your customers potentially down to the operator level) to explore the effects of multiple controls using Mathematica and another slow tool? Here is one possible idea.

Slider Model with Mathematica


When I think of an interactive slider model, Mathematica’s manipulate is one of the first tools that comes to mind.

"Baking" the Physics

Since CFD simulations are generally quite slow, I am used to postprocessing the simulation in batch mode. So, let us try to make a multidimensional flipbook using MMA in another code. Let us create a series of GIF animations varying the liquid velocity at several gas velocities while fixing the equilibrium constant in the diffusion coefficient (in Mathematica this could take a few hours as I have not optimized the process so you may want to do this overnight). When the simulation is complete, we will import the list of GIF animations.

mDic = meshfn[config -> "Co"];
basename = "varyliquid_";
filenames = 
  Table[basename <> IntegerString[i, 10, 3] <> "." <> "gif", {i, 1, 
     modelfn[md -> mDic, k -> 0.5, dratio -> 0.5, 
       pel -> yy[1, 20, 0, 50, j] // N, 
       peg -> yy[1, 20, 0, 50, i] // N, title -> "Co-Flow"][[3]], {j, 
      1, 20}
     {{"Total progress:", 
        Dynamic[f[j, 1, 20, 1]]]}, {"j=", {Dynamic@j}}}]
   "AnimationRepetitions" -> \[Infinity]], {i, 1, 20}];
files = FileNames["vary*.gif", NotebookDirectory[]];
gifs = Import[#] & /@ files;

It is quite easy to create a flipbook using manipulate so that we can view the effects along the liquid velocity and gas velocity dimensions.

m = Manipulate[
  Show[gifs[[i, j]]], {j, 1, 20, 1}, {i, 1, Length[gifs], 1}]

Co Flow Slider MMA

Just a word of warning. When I used the SaveDefinitions->True option, MMA ballooned up to over 15GB of RAM and took a while. There probably is a more optimal approach.

Likewise, we can conduct a similar process in COMSOL to create a series of GIF animations that can be imported and displayed in manipulate.


What I like about this approach is that almost all solvers can produce a GIF animation, but few or none that I am aware of can create a flipbook along multiple dimensions. Additionally, there could be other clever tricks that one could play with Mathematica using interpolation functions on a sampled mesh, image processing, or machine learning that could be incorporated into manipulate, but I have yet to experiment with those approaches. Now, we have a pathway to distill potentially terabytes of information into a model that an operator can interact with. I find that pretty interesting.

POSTED BY: Tim Laska
Posted 3 years ago

Sorry for not seeing your reply earlier. That is very interesting comparison of the two software.

I mostly do data analysis directly in COMSOL but if there is a easy way to transfer results into MMA that would be nice. (Or even better, transfer equations, geometry and parameters from MMA to COMSOL! )

Currently COMSOL is more integrated with MATLAB. The only code linking MMA and comsol I found is "mathematica-comsol" on Github. Maybe you will be interested to try it.

POSTED BY: Ox Clouding

Thank you for the github link. I will investigate when I have free moment.

In the StackExchange link provided by Kuba above, Alexi does provide a procedure to export COMSOL data so that it can be can be read into Mathematica. I used it successfully to compare Mathematica to COMSOL from my very first answer to an adsorption problem to the Wolfram Community here.

tbl = Import[
conc = tbl[[All, Range[2, 302, 2]]];
ListPlot3D[Transpose@conc, DataRange -> {{0, 2}, {0, 1.5}}, 
 PlotRange -> {{0, 1}, {0, 1.5}}, 
 ColorFunction -> (ColorData["DarkBands"][#3] &), 
 MeshFunctions -> {#2 &}, Mesh -> 20, AxesLabel -> Automatic, 
 MeshStyle -> {Black, Thick}, ImageSize -> Large]

Comsol Solution

Perhaps, one could try to automate the workflow in java on the COMSOL end to facilitate the transfer.

With regard to the mesh, it looks like one might be able to use pyNastran to help bridge the gap between COMSOL and Mathematica. I have used pyNastran in the past to extract nodes and connectivity from an optiStruct model (Nastran format) and it looks as though COMSOL will export a Nastran mesh. GUIs do help facilitate meshing operations.

Parsing a parameter file should be easy, but parsing how all the parameters are used throughout the model would be difficult and prone to errors.

POSTED BY: Tim Laska

Take a look at the extensive documentation “Advanced Solutions to Numerical Differential Equations“ and for NDSolve. There are methods to track solutions (combinations of EvaluationMonitor, StepMonitor, and Dynamic.

There are also good examples of doing solid mechanics.

POSTED BY: W. Craig Carter
Posted 3 years ago

Thanks for the reply. I am aware of MMA capabilities of manually track solutions with code. But I still wish they could make it simpler by adding a simple built-in function like TrackSolving[NDSolve[...]] instead of having us to implement it.

Just like in the past we have to write our own neural network but now we have NetTrain. Why not improve more fundamental functions like NDSolve?

POSTED BY: Ox Clouding

You may find this thread interesting:

POSTED BY: Kuba Podkalicki

@TimLaska, It is perfectly reasonable to mention use of other languages/programs in any of several situations.

(1) Validation of results.

(2) Asking about discrepancies in behavior between the two. (If the Wolfram Language is showing obviously buggy behavior, it is good also to report that directly. Also we do not take well to having the same issue reported repeatedly. But that tends to be an uncommon problem in any forum.)

(3) Discussion about linking to/from other programs.

(4) Speed comparison, when we are faster (just kidding-- if it is an apples-to-apples comparison and WL falls short, that's good for us and others to know about).

(5) Questions to the effect "Does the Wolfram Language have something similar or equivalent to the built-in function XXX from program YYY?"

I'm probably missing a few other valid reasons. In any case, it would be fine for you to add discussion of othe rprograms in this setting.

POSTED BY: Daniel Lichtblau

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

POSTED BY: Moderation Team

Thank you for the nice FEM program. The mass transfer is important in chemistry, but in other fields as well. I have to look more concisely to the analytical formulation of the mas transfer problem as trying to solve some mass transfer problems linked to peritoneal dialysis. Ladislav Lukas

POSTED BY: Ladislav Lukas

Sorry for the late reply. I hope that you may find some additional use for this post.

For the dialysis case, you could have two jump conditions on either side of the membrane as shown below.

Two Jumps Across a Membrane

So, you probably need to model the membrane and add an "interphase" region to both sides for a total of 5 regions. If the model grows too large due to meshing three relatively thin regions, then maybe a structured quad mesh is in order.

POSTED BY: Tim Laska
Anonymous User
Anonymous User
Posted 3 years ago

The strength of MMA is ______.

packages often don't share their math code and the code isn't "good" it does the job.

you can publish you mm code so mm have it, or keep it private for use in a pro product :)

POSTED BY: Anonymous User
Anonymous User
Anonymous User
Posted 3 years ago

(in short i asked if this PDE could be solved as an ODE, but wish to retract question)

POSTED BY: Anonymous User

Hi John,

Here is an attempt to answer your questions.

Generally, one has some variable of interest (e.g., temperature) that DEPENDS on INDEPENDENT variables of space (x, y, z) and time (t). The equations to describe this dependence usually derive from conservation laws (e.g., energy balance, momentum balance, species balance, etc) that are only valid for differential elements of space and time. If your system depends on more than one independent variable or if you are unsure of how many independent variables it depends on, then you are in the realm of PDEs.

As to your second question, if the problem could be represented as an ODE, then the numerical solution would have revealed it in the plots (e.g., families of characteristic contour curves would be identifiable). So, the answer is no, this system can not be represented by an ODE. For an example where I was able to convert a PDE into a simpler ODE based on characteristic lines I observed in the numerical solution of the fully coupled PDE, look here.

Almost all ODEs, PDEs, or systems thereof, cannot be solved analytically. In contrast, with the advent of fast computers and improved algorithms, a much larger class of ODE/PDE can be solved at reasonable speeds. Further, if you start off with the numerical solution, it can give insights on reasonable assumptions for a potential analytical solution.

POSTED BY: Tim Laska
Anonymous User
Anonymous User
Posted 3 years ago

Thank you very much Tim Laska. I shouldn't have asked the question I want to retract it.

My book says a PDE is used when solving a formula containing only partial differentials (not total differentials) which I should have remembered. The discussion of conversion from one to another isn't one I meant to begin.

Perhaps I was wondering when it was that lab measurement taking pushes one into a partial derivative formula rather than total: but I shouldn't start that discussion either :) It would include discussion of derivatives and integrals of several variables which aren't ... any easier discussed than pde per say.

POSTED BY: Anonymous User
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract