This problem is axially symmetric. I want to solve using the Finite Element method and taking advantage of the symmetry.
(Note: eleField3D
needs to be called with Cartesian coordinates)
Needs["NDSolve`FEM`"]
Clear["Global`*"]
cylinOffset = 0.0125;(*meter*)
cylinRadius = 0.0075;(*meter*)
shellRadius = 0.0233;(*meter*)
scaleLen = 2.10;
cylinVoltage = 100;
shell = ImplicitRegion[0 <= r <= shellRadius && ζ^2 <= (scaleLen shellRadius)^2, {r,ζ}];
halfTorus = ImplicitRegion[(r - cylinOffset)^2 + ζ^2 <= cylinRadius^2 && r <= cylinOffset, {r, ζ}];
reg = ToElementMesh[RegionDifference[shell, halfTorus], AccuracyGoal -> 9, PrecisionGoal -> 9, MaxCellMeasure -> 1*^-8];
RegionPlot[reg, AspectRatio -> scaleLen]
Please check that my 2D (r,z) Cylindrical Laplacian is correct.
solution = NDSolveValue[{(1/r) D[r D[u[r, ζ], r], r] + D[u[r, ζ], {ζ, 2}] == 0,
DirichletCondition[u[r, ζ] == 0, r == shellRadius],
DirichletCondition[u[r, ζ] == 0, ζ^2 == (scaleLen shellRadius)^2],
DirichletCondition[u[r, ζ] == cylinVoltage, r == cylinOffset && ζ^2 <= cylinRadius^2],
DirichletCondition[u[r, ζ] == cylinVoltage, (r - cylinOffset)^2 + ζ^2 == cylinRadius^2 && r <= cylinOffset]},
u, {r, ζ} ∈ reg, Method -> {"FiniteElement"}]
How do I revolve the solution to get the 3D field?
Proposed Answer:
eleFieldCylin[r_, ζ_] = -Grad[solution[r, ζ], {r, ζ}];
eleField3D[x_, y_, z_] = RotationTransform[ArcTan[x, y], {0, 0, 1}][
Insert[eleFieldCylin[Sqrt[x^2 + y^2], z], 0, 2]];
Show[
VectorPlot3D[eleField3D[x, y, z],
{x, -shellRadius, shellRadius}, {y, -shellRadius, shellRadius},
{z, -(scaleLen shellRadius), (scaleLen shellRadius)}],
RegionPlot3D[
DiscretizeRegion[
ImplicitRegion[(cylinOffset - Sqrt[x^2 + y^2])^2 + z^2 <=
cylinRadius^2 && x^2 + y^2 <= cylinOffset^2, {x, y, z}]]]
]
Answer
I think the proposed answer is correct now, here's just some suggestions.
First, you can deduce Laplace equation with the help of Laplacian
to avoid possible mistakes:
Laplacian[u[r, ζ], {r, th, ζ}, "Cylindrical"] == 0
(* Derivative[0, 2][u][r, ζ] + Derivative[1, 0][u][r, ζ]/r +
Derivative[2, 0][u][r, ζ] == 0 *)
Second, eleField3D
can be defined in a clearer and simpler way with TransformField
:
eleFieldCylin[r_, ζ_] = -Grad[solution[r, ζ], {r, th, ζ}, "Cylindrical"];
eleField3D[x_, y_, z_] =
TransformedField["Cylindrical" -> "Cartesian",
eleFieldCylin[r, ζ], {r, th, ζ} -> {x, y, z}]
Third, it may be better to set a region function for VectorPlot3D
, for example:
regionfunc = (Function[{x, y, z}, #1] & )[
0 <= r <= shellRadius && ζ^2 <= (scaleLen*shellRadius)^2 &&
! ((r - cylinOffset)^2 + ζ^2 <= cylinRadius^2 &&
r <= cylinOffset) /. {r -> Sqrt[x^2 + y^2], ζ -> z}];
VectorPlot3D[
(Piecewise[{{#1, regionfunc[x, y, z]}},
{0, 0, 0}] & )@eleField3D[x, y, z]//Evaluate, {x, -shellRadius, shellRadius},
{y, -shellRadius, shellRadius}, {z, -(scaleLen shellRadius), scaleLen shellRadius}]
4th, I think PrecisionGoal
and AccuracyGoal
has no obvious effect here and can be taken away. (Have you read this post?: Is manual adjustment of AccuracyGoal and PrecisionGoal useless?)
Comments
Post a Comment