Coordenadas cilíndricas en FEM

Si encuentras alguna parte que no entiendes puedes comentarlo y haremos todo lo necesario de ayudarte rápidamente.

Solución:

Esto parece ser un problema de flujo impulsado por la tapa. Estoy de acuerdo con la perspectiva de @ user21 de que debes resolver esto en coordenadas cartesianas. Debería simplificar la especificación de la condición de contorno. Dado que el sistema está cerrado, deberá definir la presión en un nodo. Usé OpenCascade para construir el medio cilindro. Aquí está el flujo de trabajo.

(* Load Required Packages *)
Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
(* Use OpenCascade To Make Half Sym Geometry *)
pp = Polygon[0, 0, -1, 0, 0, 1, 1, 0, 1, 1, 0, -1];
shape = OpenCascadeShape[pp];
axis = 0, 0, 0, 0, 0, 1;
sweep = OpenCascadeShapeRotationalSweep[shape, axis, -Pi];
(* Create Mesh *)
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> "Length" -> .075, 
   "IncludePoints" -> 0, 0.5, -1];
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
(* Create PDE System *)
ClearAll[μ]
op = Inactive[
              Div][(-μ, 0, 0, 0, -μ, 0, 0, 
                    0, -μ.Inactive[Grad][
         u[x, y, z], x, y, z]), x, y, 
              z] + 
     D[p[x, y, z], x], 
        Inactive[
              Div][(-μ, 0, 0, 0, -μ, 0, 0, 
                    0, -μ.Inactive[Grad][
         v[x, y, z], x, y, z]), x, y, 
              z] + 
     D[p[x, y, z], y],
        Inactive[
              Div][(-μ, 0, 0, 0, -μ, 0, 0, 
                    0, -μ.Inactive[Grad][
         w[x, y, z], x, y, z]), x, y, 
              z] + 
     D[p[x, y, z], z], 
    D[u[x, y, z], x] + 
     D[v[x, y, z], y] + 
     D[w[x, y, z], z] /. μ -> 1;
pde = op == 0, 0, 0, 0;
bcs = DirichletCondition[
        u[x, y, z] == 1, v[x, y, z] == 0., w[x, y, z] == 0., 
    z == 1.],
      DirichletCondition[
        u[x, y, z] == 0, v[x, y, z] == 0., w[x, y, z] == 0., 
        z == -1. ;
(* Solve PDE *)
xVel, yVel, zVel, pressure = 
    NDSolveValue[pde, bcs, u, v, w, p, x, y, z ∈ mesh, 
      Method -> "FiniteElement", 
          "InterpolationOrder" -> u -> 2, v -> 2, w -> 2, p -> 1];
(* Visualize Solution *)
surf = "YStackedPlanes", 0, "ZStackedPlanes", -1, 1;
Show[SliceContourPlot3D[
    [email protected]xVel[x, y, z], yVel[x, y, z], zVel[x, y, z], 
    surf, x, y, z ∈ mesh, PlotPoints -> 50, 
    BoxRatios -> Automatic, ColorFunction -> "TemperatureMap"], 
  ImageSize -> Medium, ViewPoint -> Front]
DensityPlot3D[
  Norm[xVel[x, y, z], yVel[x, y, z], zVel[x, y, z]], x, y, 
      z ∈ mesh, BoxRatios -> Automatic, 
  ColorFunction -> "TemperatureMap", ViewAngle -> 0.3669386546105606`,
  ViewPoint -> 3.7435513617679828`, 1.2106476957796874`, 
   0.9258298223054351`, 
 ViewVertical -> 0.27079048490259205`, 0.14735018657087556`, 
   0.9512940848148628`]
SliceVectorPlot3D[xVel[x, y, z], yVel[x, y, z], 
    zVel[x, y, z], surf, x, y, z ∈ mesh, 
 VectorPoints -> 20,
   VectorColorFunction -> "BrightBands", BoxRatios -> Automatic, 
 ViewPoint -> Front]

Cualitativamente, concuerda con el modelo COMSOL que preparé.

Aquí hay una versión en coordenadas cartesianas para comenzar:

reg = Cylinder[0, 0, 0, 0, 0, 1, 1];

a = IdentityMatrix[3];
stokesFlowOperator = Inactive[Div][
     a.Inactive[Grad][u[x, y, z], x, y, z], x, y, z] - 
    D[p[x, y, z], x], 
   Inactive[Div][a.Inactive[Grad][v[x, y, z], x, y, z], x, y, z] -
     D[p[x, y, z], y], 
   Inactive[Div][a.Inactive[Grad][w[x, y, z], 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];
[CapitalGamma]D = 
   DirichletCondition[u[x, y, z] == 1., v[x, y, z] == 0., 
     w[x, y, z] == 0., x == 1], 
   DirichletCondition[u[x, y, z] == 0., v[x, y, z] == 0., 
     w[x, y, z] == 0., x < 1], 
   DirichletCondition[p[x, y, z] == 0, x == -1 && y == 0 && z == 1];

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[reg];

xVel, yVel, zVel, pressure = 
  NDSolveValue[stokesFlowOperator == 0, 0, 0, 
      0, [CapitalGamma]D, u, v, w, p, x, y, z [Element] mesh, 
   Method -> "FiniteElement", 
     "InterpolationOrder" -> u -> 2, v -> 2, w -> 2, p -> 1];

Debería pensar más en las condiciones de contorno, especialmente la condición de presión.

rmf = RegionMember[MeshRegion[mesh]];
Quiet[VectorPlot3D[xVel[x, y, z], yVel[x, y, z], zVel[x, y, z], 
  Evaluate[Sequence @@ Join[x, y, z, mesh["Bounds"]*1.01, 2]],
   VectorStyle -> "Arrow3D", VectorColorFunction -> "TemperatureMap", 
  VectorScale -> Tiny, Scaled[0.4], None, VectorPoints -> 9, 9, 9,
   Axes -> None, Boxed -> False, 
  RegionFunction -> (rmf[#1, #2, #3] &)], 
 InterpolatingFunction::femdmval]

Reseñas y puntuaciones

Eres capaz de respaldar nuestro cometido poniendo un comentario y puntuándolo te damos la bienvenida.

---

---

Utiliza Nuestro Buscador