2
$\begingroup$

I am trying to duplicate https://www.wolfram.com/language/12/nonlinear-finite-elements/navier-stokes-equation.html?product=language Naiver-Stokes type simulation for a charge current density. If I am not wrong, I can use inviscid flow. How can I write wall boundary conditions where the normal velocity is zero?

\[CapitalOmega]=RegionUnion[Rectangle[{0.02,0},{0.1,0.0101}],Rectangle[{0,0.0049},{0.02,0.0101}]];
Show[RegionPlot[\[CapitalOmega], Axes -> False, Frame -> None, 
  ImageSize -> Large], 
 VectorPlot[Evaluate[{1, 0}], {x, y} \[Element] \[CapitalOmega], 
  Sequence[VectorScale -> Small, VectorStyle -> Red, 
   VectorMarkers -> Placed["Arrow", "Start"], 
   VectorPoints -> Table[{0, y}, {y, 0.005, 0.01, 0.001}]]], 
 AspectRatio -> Automatic]
inflowBC = DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, x == 0];
outflowBC = 
  DirichletCondition[{u[x, y] == 0.5, v[x, y] == 0}, x == 0.1];
wallBC = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
  0 < x < 0.1]; (*This is wrong condition. How can I convert to \
inviscid flow?*)
bcs = {inflowBC, outflowBC, wallBC};
{uVel, vVel} = NDSolveValue[{
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] == 0, 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] - 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] == 0, bcs}, {u, 
    v}, {x, y} \[Element] \[CapitalOmega], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.00000025}}];
$\endgroup$
3
  • $\begingroup$ Hi Emrah. You are in 2D. So in that case, v[x,y] is the "normal" velocity along the horizontal wall. Thus half of your wall BC is correct, namely that v[x,y]==0 along the inlet. I guess you don't want u[x,y]=0 there, which is not right for inviscid flow. I suppose you could try setting u[x,y]=1 there, instead? Let me know if that helps. $\endgroup$ Apr 19 '20 at 18:03
  • $\begingroup$ Thanks a lot for the point. Yes, it will be useful to try. Indeed, I also want to try circular shapes. I have to think about more general solution. I know comsol can solve similar systems. The problem does not seem that hard, but may be I am missing some physical points. I will look into potential and stream scalar options, and use Laplacian. $\endgroup$
    – emrah
    Apr 19 '20 at 19:59
  • $\begingroup$ I had a similar problem recently (in fact, am still working on it). I had good success with both potential and streamfunction formulations (I have been learning as I go along). Note that in potential formulation, you can use Neumann boundary conditions for the velocity normal components at boundaries, where you need them. $\endgroup$ Apr 19 '20 at 20:53
3
$\begingroup$

When I plotted the OP's self answered question, I saw vertical flow versus the left to right flow that I was expecting.

ContourPlot[sol[x, y], {x, y} ∈ Ω]
VectorPlot[
 Evaluate@Grad[sol[x, y], {x, y}], {x, y} ∈ Ω]

OP Solution

As @Paul Harrison alluded to in the comments, that simpler formulation may be to use a NeumannValueflux values at the inlet or outlet boundaries. We will define a NeumannValue flux condition at the inlet and an essential DirichletCondition. Also, we will recast the equations in coefficient form so material flows from high potential to low potential. Here is an example workflow.

Ω = 
  RegionUnion[Rectangle[{0.02, 0}, {0.1, 0.0101}], 
   Rectangle[{0, 0.0049}, {0.02, 0.0101}]];
pressurePoint = DirichletCondition[phi[x, y] == 0, x == 0.1];
nv = NeumannValue[1, x == 0];
bcs = {pressurePoint};
sol = NDSolveValue[{Inactive[
       Div][{{-1, 0}, {0, -1}}.Inactive[Grad][phi[x, y], {x, y}], {x, 
       y}] == nv, bcs}, phi, {x, y} ∈ Ω, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.00000025}}];
vel[x_, y_] := Evaluate@Grad[-sol[x, y], {x, y}]
nvel[x_, y_] := {vel[x, y][[2]], -vel[x, y][[1]]}
Show[ContourPlot[sol[x, y], {x, y} ∈ Ω, 
  Sequence[PlotRange -> All, Frame -> None, Axes -> None, 
   Contours -> 40, ColorFunction -> "TemperatureMap", 
   ContourStyle -> Black]], Graphics[RegionBoundary[Ω]],
  AspectRatio -> Automatic, ImageSize -> Medium]
VectorPlot[vel[x, y], {x, y} ∈ Ω, 
 AspectRatio -> Automatic]

Left to right flow

Now, the flow is left to right as I initially expected. We can zoom in on the backward facing step to look at more details of the potential and the flows.

roi = RegionIntersection[ 
   Rectangle[{0.015, 0}, {0.025, 0.01}], Ω];
rmfroi = RegionMember[roi];
bounds = RegionBounds[roi];
vals = sol["ValuesOnGrid"];
mesh = sol["ElementMesh"];
coords = mesh["Coordinates"];
subRange = MinMax@Pick[vals, rmfroi[coords]];
legendBar = 
  BarLegend[{"TemperatureMap", subRange}, 50, LegendLayout -> "Row", 
   LegendLabel -> Style["[phi]", Opacity[0.6`]], 
   LegendMarkerSize -> 500, Charting`TickSide -> Left];
options = {PlotRange -> subRange, 
   ColorFunction -> ColorData[{"TemperatureMap", subRange}], 
   ContourStyle -> Opacity[0.1`], ColorFunctionScaling -> False, 
   Contours -> 30, AspectRatio -> Automatic,
PlotPoints -> 41, FrameLabel -> {"x", "y"}, 
   PlotLabel -> Style["Potential: phi(x,y)", 18], ImageSize -> 500};
Column[{legendBar, 
  ContourPlot[sol[x, y], {x, y} ∈ roi, Evaluate[options]]}, 
 Spacings -> 0]
VectorPlot[vel[x, y], {x, y} ∈ roi, AspectRatio -> Automatic]

Zoomed image

$\endgroup$
0
$\begingroup$

I got the solution by help of a CFD friend. Using stream function makes much more sense. Some background. http://cfd.mace.manchester.ac.uk/twiki/pub/Main/TimCraftNotes_All_Access/fl3-inviscid.pdf

\[CapitalOmega] = 
  RegionUnion[Rectangle[{0.02, 0}, {0.1, 0.0101}], 
   Rectangle[{0, 0.0049}, {0.02, 0.0101}]];
inflowBC = 
  DirichletCondition[phi[x, y] == 5.82688 (y - 0.0049)/0.0101, 
   x == 0];
outflowBC = DirichletCondition[phi[x, y] == 3 y/0.0101, x == 0.1];
wallBC1 = DirichletCondition[phi[x, y] == 3, y > 0.007]; 
wallBC2 = DirichletCondition[phi[x, y] == 0, y < 0.007]; 
bcs = {inflowBC, outflowBC, wallBC1, wallBC2};
sol = NDSolveValue[{D[phi[x, y], x, x] + D[phi[x, y], y, y] == 0, 
   bcs}, phi, {x, y} \[Element] \[CapitalOmega]]
$\endgroup$

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.