1D-Netzgenerierung für PDE-Lösung
Ich versuche ein System von zwei PDE zu lösen, die von Zeit und Entfernung abhängig sind (H [x, t] und P [x, t]). Ich löse das Problem mit der Linienmethode, möchte aber das Netz selbst generieren und in NDsolve einführen. Das Netz, das ich generieren möchte, ist das folgende
Ich brauche ein Netz wie dieses, weil sich die Werte für eine der Funktionen (P [x, t]) mit der Zeit nur sehr nahe an x = 0 ändern, während sich H [x, t] über den gesamten Bereich 0 <x <ändert xmax. Unten finden Sie ein Beispiel für den Code, den ich verwende
(* Constants *)
f = 38.94; logL = -2;
Ls = 10^logL; a = 0.5;
C1 = 1*^-5; dH = 1*^-6;
Ea = 0.1;
tmax = 40; (* Time in seconds *)
xmax = 10 Sqrt[dH] Sqrt[tmax]; (* Maximum distance to simulate. cm *)
(* PDE system *)
eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] == NeumannValue[Ls Exp[a f Ea ] P[x, t] - Ls Exp[-a f Ea ] H[x, t],
x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == NeumannValue[-Ls Exp[a f Ea ] P[x, t] + Ls Exp[-a f Ea ] H[x, t],
x == 0], P[x, 0] == 1};
(*Solution of the differential equations*)
prec = 7;
msf = 0.001;
sol = NDSolve[{eqsH, eqsP}, {H, P}, {x, 0, xmax}, {t, 0, tmax},
AccuracyGoal -> prec, PrecisionGoal -> prec,
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}] // First //
Quiet;
Kann ich bitte Hilfe beim Erstellen und Einführen des Netzes in NDSolve erhalten? Danke im Voraus !
Antworten
Hier ist ein alternativer Ansatz unter Verwendung eines abgestuften Netzes.
Definieren Sie einige Hilfsfunktionen für ein abgestuftes Netz
Hier sind einige Funktionen, mit denen ich anisotrope 1d- bis 3D-Netze erstellt habe. Nicht alle Funktionen werden verwendet.
(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(* Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=
Quiet@Abs@
FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 100000},
Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]
]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]
Erstellen Sie ein abgestuftes horizontales Netzsegment
Im Folgenden wird ein horizontaler Netzbereich mit 100 Elementen erstellt, wobei die anfängliche Elementbreite 1/10000 der Domänenlänge beträgt.
(*Create a graded horizontal mesh segment*)
(*Initial element width is 1/10000 the domain length*)
seg = leftSegmentGrowth[xmax, 100, xmax/10000];
Print["Horizontal segment"]
rh = pointsToMesh@seg
(*Convert mesh region to element mesh*)
(*Extract Coords from horizontal region*)
crd = MeshCoordinates[rh];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["ListPlot of exponential growth of element size"]
ListPlot[Transpose@mesh["Coordinates"]]
Man kann das exponentielle Wachstum der Elementgröße sehen, wenn die Elementzahl zunimmt.
Konvertieren Sie den Netzbereich in ein Elementnetz und lösen Sie die PDE
Im Allgemeinen konvertiere ich das MeshRegionin ein 'ElementMesh'so, damit ich bei Bedarf Element- und Punktmarkierungen anwenden kann.
(*Solve PDE on graded mesh*)
{Hfun, Pfun} =
NDSolveValue[{eqsH, eqsP}, {H, P}, x ∈ mesh, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}];
(*Animate Hfun solution*)
imgs = Plot[Hfun[x, #], x ∈ mesh,
PlotRange -> {0.9999999, 1.0018}] & /@ Subdivide[0, tmax, 120];
Print["Animation of Hfun solution"]
ListAnimate@imgs
Anhang: Beispiele für anisotrope Vernetzung
Wie ich im Kommentar unten angedeutet habe, zeigt die Aufzählungsliste unten einige Beispiele, bei denen ich anisotropes Quad-Meshing verwendet habe, um scharfe Schnittstellen zu erfassen, die sonst sehr rechenintensiv wären. Der Code ist funktionsfähig, aber nicht optimal, und einige der Funktionen wurden im Laufe der Zeit geändert. Benutzung auf eigene Gefahr
- 2D-stationär
- Mathematica vs. MATLAB: Warum erhalte ich unterschiedliche Ergebnisse für PDE mit nicht konstanter Randbedingung?
- Verbesserung der Konvergenz von Mesh- und NDSolve-Lösungen
- 2D-Transient
- Steuern der dynamischen Zeitschrittgröße in NDSolveValue
- Wie modelliere ich die Diffusion durch eine Membran?
- Massentransport FEM mit Quad Mesh
- NDSolve mit Gleichungssystem mit unbekannten Funktionen, die in verschiedenen Domänen definiert sind
- 3D-Vernetzung
- Erstellen Sie ein abgestuftes Netz
- 3D-stationär
- Wie kann die FEM-Lösung mit NDSolve verbessert werden?
- 3D-FEM-Vektorpotential
Wenn Sie Zugriff auf andere Tools wie COMSOL haben, die über Grenzschichtfunktionen verfügen, können Sie Netze über die Ressourcenfunktion FEMAddOns importieren . Es funktioniert nicht für 3D-Netze, für die zusätzliche Elementtypen wie Prismen und Pyramiden erforderlich sind, die derzeit in Mathematicas FEM nicht unterstützt werden .
Was ist damit?
lst1 = Partition[
Join[Table[0.01*i, {i, 0, 5}], Table[0.1*i, {i, 0, 15}]], 1];
lst2 = Table[{i, i + 1}, {i, 1, Length[lst1] - 1}];
<< NDSolve`FEM`
mesh2 = ToElementMesh["Coordinates" -> lst1,
"MeshElements" -> {LineElement[lst2]}]
(* ElementMesh[{{0., 1.5}}, {LineElement["<" 21 ">"]}] *)
Lassen Sie es uns visualisieren:
mesh2["Wireframe"["MeshElementIDStyle" -> Red]]
Die roten Zahlen zeigen die Netzelemente an. Die Stelle, an der sie sich überlappen, ist tatsächlich die Stelle, an der das Netz zehnmal dichter ist (siehe das vergrößerte Bild unten):
Habe Spaß!