1
$\begingroup$

I want to solve a simple diffuse equation in a cylinder.

$\nabla\cdot(\hat{c}\nabla\rho_P)=0$

$\rho_P(x,y,\theta)==0,\sqrt{x^2+y^2}=R$

$\rho_P(0,0,\theta)==1,\sqrt{x^2+y^2}=0$

$\rho_P(x,y,2\pi)=\rho_P(x,y,0)$

$\hat{c}=\begin{pmatrix} D_t& 0& 0 \\ 0 & D_t &0\\0 &0 &D_r \end{pmatrix} \quad$

$D_t: transportation\ coefficient$

$D_r: rotation\ coefficient$

Needs["NDSolve`FEM`"]

{R, Difft, Diffr} = {29.5, 1, 0.01};
\[CapitalOmega] = Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi]}}, R];

c = {{Difft, 0, 0}, {0, Difft, 0}, {0, 0, Diffr}};

With[{rhop = rhop[x, y, \[Theta]]}, 
  op1 = Div[c . Grad[rhop, {x, y, \[Theta]}], {x, y, \[Theta]}]];

Subscript[\[CapitalGamma], D1] = 
  DirichletCondition[rhop[x, y, \[Theta]] == 1, 
   x^2 + y^2 == 0 && \[Theta] != 2 \[Pi]];
Subscript[\[CapitalGamma], D2] = 
  DirichletCondition[rhop[x, y, \[Theta]] == 0, 
   x^2 + y^2 > 0 && \[Theta] != 0 && \[Theta] != 2 \[Pi]];
Subscript[\[CapitalGamma], p] = 
  PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] == 2 \[Pi],
    Function[\[Theta], \[Theta] - 2 \[Pi]]];

cylinderMesh = 
  ToElementMesh[\[CapitalOmega], "MaxBoundaryCellMeasure" -> 0.04];

npfun = NDSolveValue[{op1 == 0, Subscript[\[CapitalGamma], p], 
    Subscript[\[CapitalGamma], D1], Subscript[\[CapitalGamma], D2]}, 
   rhop, {x, y, \[Theta]} \[Element] cylinderMesh];

ContourPlot[npfun[x, y, 2 \[Pi]], {x, y} \[Element] Disk[{0, 0}, R], 
 PlotRange -> All]

Obviously, it is a symmetrical function, but when I solve it using NDSolveValue, the result is not symmetric. enter image description here

I think I can solve it by adding

additionalPoints = {{0, 0, 0}, {0, 0, 2 \[Pi]}}
cylinderMesh = 
  ToElementMesh[\[CapitalOmega], "IncludePoints" -> additionalPoints];

But the error: enter image description here 无法产生网络:enable to generate mesh

$\endgroup$
7
  • $\begingroup$ Your second DirichletCondition should be Subscript[\[CapitalGamma], D2] = DirichletCondition[rhop[x, y, \[Theta]] == 0, x^2 + y^2 ==R^2 && \[Theta] != 0 && \[Theta] != 2 \[Pi]] $\endgroup$ Mar 16, 2023 at 9:38
  • $\begingroup$ Region Cylinder suggest a 3D problem, but I'm missing an axial z coordinate!? $\endgroup$ Mar 16, 2023 at 11:05
  • $\begingroup$ but,if i write as Subscript[\[CapitalGamma], D1] = DirichletCondition[rhop[x, y, \[Theta]] == 1, x^2 + y^2 == 0 && \[Theta] != 0 && \[Theta] != 2 \[Pi]];, Subscript[\[CapitalGamma], D2] = DirichletCondition[rhop[x, y, \[Theta]] == 0, x^2 + y^2 == R^2 && \[Theta] != 0 && \[Theta] != 2 \[Pi]]; there is a error: NDSolveValue::bcnop: No position was found on the boundary where x^2+y^2==0&&\[Theta]!=0&\[Theta]!=2 \[Pi] is true, so DirichletCondition[rhop==1,x^2+y^2==0&&\[Theta]!=0&&[Theta]!=2 \[Pi]]. will be ignored $\endgroup$
    – 江蛮子
    Mar 17, 2023 at 0:58
  • $\begingroup$ I think I should add a line(0,0,theta) or points(0,0,2pi)(0,0,0) to the mesh, but i can't find the option. $\endgroup$
    – 江蛮子
    Mar 17, 2023 at 1:04
  • 1
    $\begingroup$ The problem is not well-posed. Subscript[\[CapitalGamma], D2] already determines the solution, you cannot force b.c. at x^2+y^2==0, because it's the removable singularity in cylindrical coordinates. Related: mathematica.stackexchange.com/questions/245195/… $\endgroup$
    – xzczd
    Mar 17, 2023 at 12:44

1 Answer 1

2
$\begingroup$

I modified your code a little bit. You need two periodic boundaryconditions with correct TranslationTransform in all coordinates. Additionally I force NDSolve to use FiniteElement- Method (no need to load Needs["NDSolveFEM"])

I changed the Region \[CapitalOmega] to avoid singularity x^2+y^2->0. Play around with the new small parameter dR.

It is often recommended to rationalize the parameters.

{R, dR, Difft, Diffr} = {29.5, 29.5/100, 1, 0.01}// Rationalize[#, 0]&;
\[CapitalOmega] =RegionDifference[Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi] }},R], Cylinder[{{0, 0, 0}, {0, 0, 2 \[Pi] }}, dR]];

Try

c = {{Difft, 0, 0}, {0, Difft, 0}, {0, 0, Diffr}};

With[{rhop = rhop[x, y, \[Theta]]},op1 = Div[c . Grad[rhop, {x, y, \[Theta]}], {x, y, \[Theta]}]];

\[CapitalGamma]D1 =DirichletCondition[rhop[x, y, \[Theta]] == 1,x^2 + y^2 == dR^2 && 0 < \[Theta] < 2 \[Pi]];
\[CapitalGamma]D2 =DirichletCondition[rhop[x, y, \[Theta]] == 0,x^2 + y^2 == R^2 && 0 < \[Theta] < 2 \[Pi]];
\[CapitalGamma]p = {PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] >= 2 \[Pi] &&dR^2 < x^2 + y^2 < R^2, TranslationTransform[{0, 0, -2 Pi}]],PeriodicBoundaryCondition[rhop[x, y, \[Theta]], \[Theta] <= 0 && dR^2 <x^2 + y^2 < R^2, TranslationTransform[{0, 0, +2 Pi}]]};


npfun = NDSolveValue[{op1 ==0, \[CapitalGamma]p, \[CapitalGamma]D1, \[CapitalGamma]D2}, rhop, {x, y, \[Theta]} \[Element] \[CapitalOmega], 
Method -> {"FiniteElement"}]

reg = ImplicitRegion[dR^2 <= x^2 + y^2 <= R^2 , {x, y}]
Plot3D[ npfun[x, y,  \[Pi]]  , {x, y} \[Element] reg,PlotRange -> All , Mesh -> False]

enter image description here

NDSolve now evaluates a solution which is plausible.

The periodic boundary conditions are fullfiled quite well:

Plot3D[Evaluate[npfun[x, y, 2 \[Pi]] - npfun[x, y, 0]], {x, y}\[Element] reg,PlotRange -> All]

enter image description here

Plot3D[Evaluate[Derivative[0, 0, 1][npfun][x, y, 2 \[Pi]] -Derivative[0, 0, 1][npfun][x, y, 0]], {x, y} \[Element] reg, PlotRange-> All]

enter image description here

I think you might get much better results if you formulate the underlying problem in "Cylindrical"-coordinates!

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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