Skip to main content

differential equations - Mismatch between Mathematica and COMSOL in 3D FEM problem


I would like to solve an advection-diffusion problem on a torus domain. There are three Dirichlet conditions: One at the inlet (concentration $c=0$), one at the outlet ($c=0.5$) and one at the wall ($c=1$).



The geometry is set up in the following code. Note that RegionCapillary is the full solution domain, RegionCapillary2D is a 2D projection for convenient plotting and the list points goes through the centre of the tube. In the plot below RegionCapillary is orange and RegionCapillary2D is blue.


radiusMajor = 4;
radiusMinor = 1;
RegionCapillary = 
  ImplicitRegion[( radiusMajor - Sqrt[x^2 + y^2])^2 + z^2 <= 
     radiusMinor^2 && x >= 0, {x, y, z}];
RegionCapillary2D = 
  ImplicitRegion[( radiusMajor - Sqrt[x^2 + y^2])^2 <= radiusMinor^2 &&
     x >= 0, {x, y}];
RP = RegionPlot3D[

   DiscretizeRegion /@ {RegionCapillary, RegionCapillary2D} // 
    Evaluate, Axes -> True, AxesLabel -> {x, y, z}, 
   PlotStyle -> {Opacity[0.4], Opacity[1]}];
LP = ListPointPlot3D[
   points = 
    Reverse@Table[{radiusMajor Cos[\[Phi]], radiusMajor Sin[\[Phi]] , 
       0}, {\[Phi], -Pi/2, Pi/2,  Pi/1000.0}], 
   PlotStyle -> PointSize[0.03]];
Show[RP, LP]


enter image description here


I solve the problem very similarly to user21's reply to my question (here). There are two parts to the problem: Stokes flow and then an advection-diffusion in which the previous (Stokes) solution provides the fluid velocity.


The Stokes problem is:


<< NDSolve`FEM`
StokesOp := {
   -Inactive[Laplacian][u[x, y, z], {x, y, z}] + D[p[x, y, z], x],
   -Inactive[Laplacian][v[x, y, z], {x, y, z}] + D[p[x, y, z], y],
   -Inactive[Laplacian][w[x, y, z], {x, y, z}] + D[p[x, y, z], z],
   Div[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}]};


inletCondition[x_, y_, z_] := x == 0 && y >= 0
outletCondition[x_, y_, z_] := x == 0 && y < 0
capillaryCondition[x_, y_, z_] := x > 0

bcs := {
   DirichletCondition[p[x, y, z] == 1, inletCondition[x, y, z]],
   DirichletCondition[p[x, y, z] == 0, outletCondition[x, y, z]],
   DirichletCondition[{u[x, y, z] == 0., v[x, y, z] == 0., 
     w[x, y, z] == 0.}, capillaryCondition[x, y, z]]};


VelMMA[x_, y_, z_] := Sqrt[
 xVel[x, y, z]^2 + yVel[x, y, z]^2 + zVel[x, y, z]^2]

Clear[xVel, yVel, pressure];
AbsoluteTiming[{xVel, yVel, zVel, pressure} = 
   NDSolveValue[{StokesOp == {0, 0, 0, 0}, bcs}, {u, v, w, 
     p}, {x, y, z} \[Element] RegionCapillary, 
    Method -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 5 10^-3}}, 

    "ExtrapolationHandler" -> {0. &, "WarningMessage" -> False}];]

The advection-diffusion problem is:


Pe = 0.7*10^4;

eqn = Pe {xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]}.Grad[
      c[x, y, z], {x, y, z}] - Laplacian[c[x, y, z], {x, y, z}];

bcs = {
   DirichletCondition[c[x, y, z] == 0, inletCondition[x, y, z]],

   DirichletCondition[c[x, y, z] == 0.5, outletCondition[x, y, z]],
   DirichletCondition[c[x, y, z] == 1, capillaryCondition[x, y, z]]
   };

AbsoluteTiming[
 conMMA = NDSolveValue[{eqn == 0, bcs}, 
    c, {x, y, z} \[Element] RegionCapillary, 
    Method -> {"FiniteElement", "InterpolationOrder" -> {c -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 10^-3}}, 
    "ExtrapolationHandler" -> {0. &, "WarningMessage" -> False}];]


The velocity magnitude and concentration are plotted like this:


DensityPlot[
 Sqrt[xVel[x, y, 0]^2 + yVel[x, y, 0]^2 + zVel[x, y, 0]^2], {x, 
   y} \[Element] DiscretizeRegion@RegionCapillary2D, PlotRange -> All,
  PlotLabel -> "velocity magnitude", PlotPoints -> 100, 
 ColorFunction -> "TemperatureMap", AspectRatio -> 2, 
 PlotLegends -> Automatic]
DensityPlot[
 conMMA[x, y, 0], {x, y} \[Element] 

  DiscretizeRegion@RegionCapillary2D, PlotRange -> All, 
 PlotLabel -> "concentration", PlotPoints -> 100, 
 ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, 
 AspectRatio -> 2]

enter image description here


I have compared these results with an identical problem setup in COMSOL, finding that COMSOL's solution looks smoother and more plausible: enter image description here


If you trace the velocity and concentration along the centreline which goes right through the tube, you can see that while both solutions satisfy $c=0.5$ at the outlet (i.e. the right hand side of the ListPlots), Mathematica has a very weird behaviour for the diffusion problem. The Stokes problem is largely the same between the two software packages: enter image description here I am puzzled by this discrepancy and would very much appreciate some help. My suspicion is that I am specifying the boundary conditions in Mathematica incorrectly.


EDIT:


With a finer mesh on both Mathematica and COMSOL, the solutions do converge somewhat, but the behaviour near the outlet of the Mathematica solution is still somewhat puzzling: enter image description here



EDIT 2: Do things get better for lower Péclet number?


Yes. I made a comparison for Pe=10^3(original version is Pe = 0.7*10^4), finding a better convergence and less weirdness near the boundary:


enter image description here enter image description here


The weird effects at the boundary of the DensityPlot are just interpolation artifacts from importing Comsol results from text file.


EDIT 3: A more complex example


OK, if I'm doing this, might as well go all in. Here is a more complex example of the tube, based on a looping segment with quite irregular cross-sections. Here it is plotted with centreline as well as planes which I use in inequalities to set boundary conditions in Mathematica:


enter image description here


A major difference to the above examples is that at the outlet there is a Neumann condition, $\boldsymbol{n}\cdot \nabla c = 0$, which in Mathematica is NeumannValue[0, pred]. Here is a plot of concentration and velocity accross the centreline, just as in the original examples:


enter image description here


You can see that once again the Stokes problem is a superb match, but the advection-diffusion problem has some serious discrepancy between solutions. Still, it is interesting that the curves are so similar and only appear to be offset. The negative concentration produced near the inlet by Mathematica is however unphysical.



Here is a pretty visualisation of with SliceContourPlot, which shoes how similar the cross-section profiles look:


enter image description here


EDIT 4: A conclusion on 3D FEM with imported geometries


I ended up working a lot on similarly complex imported 3D geometries as above, switching completely to Comsol was essential. The boundary conditions were getting so tricky it really helped having a visual interface. The meshing control on boundary layers was also helpful, and their constructive geometry (boolean operations on discrete regions) works ok. I ended up being very careful about ensuring flux conservation in Comsol (The best way to evaluate advective and diffusive flux is via the tds.ncflux_c and tds.ndflux_c). My simulation results to similar problems to edit 3, made in Comsol and checked very carefully, can be seen in papers here and here. In fact, in the first linked paper we even looked fluxes at cross-sections of sub-meshes ($N_j$) of quite complex vasculature, so a visual interface was enormously helpful for that:


enter image description here´


As for Mathematica, I gave up on trying to solve complex advection-diffusion problems on discrete 3D domains, simply because the meshing got too complicated. If you import a discrete BoundaryMesh in MMA, the boundary cannot easily refined: As it does not come from predefined MMA regions, options like MaxCellMeasure and AccuracyGoal cannot be used on the boundary (related question here). Maybe the workflow with an external mesher like gmsh is the way to go here, but I haven't tried.




Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Adding a thick curve to a regionplot

Suppose we have the following simple RegionPlot: f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}] Now I'm trying to change the curve defined by $y=g[x]$ into a thick black curve, while leaving all other boundaries in the plot unchanged. I've tried adding the region $y=g[x]$ and playing with the plotstyle, which didn't work, and I've tried BoundaryStyle, which changed all the boundaries in the plot. Now I'm kinda out of ideas... Any help would be appreciated! Answer With f[x_] := 1 - x^2 g[x_] := 1 - 0.5 x^2 You can use Epilog to add the thick line: RegionPlot[{y < f[x], f[x] < y < g[x], y > g[x]}, {x, 0, 2}, {y, 0, 2}, PlotPoints -> 50, Epilog -> (Plot[g[x], {x, 0, 2}, PlotStyle -> {Black, Thick}][[1]]), PlotStyle -> {Directive[Yellow, Opacity[0.4]], Directive[Pink, Opacity[0.4]],
Post a Comment
Read more