13
|
20466 Views
|
2 Replies
|
14 Total Likes
View groups...
Share
GROUPS:

# FEM Solver for Navier-Stokes equations in 2D

Posted 9 years ago
 I'll show how to use the low level FEM functionality to code up a non-linear Navier-Stokes solver. The documentation explains the details about the low level FEM programming functionality which I use here. Here is the basic idea: After every non-linear iteration we re-create an interpolation function from the now current solution vector and re-insert those into the PDE coefficients and iterate until converged. This will not be insanely efficient but it works on a PDE level. Now, to tackle non-linear problems it's a good idea to get the linear version to work first. In this case this is a Stokes solver. Here is a utility function to convert a PDE into it's discretized version: Needs["NDSolveFEM"] PDEtoMatrix[{pde_, ?___}, u_, r__] := Module[{ndstate, feData, sd, bcData, methodData, pdeData}, {ndstate} = NDSolveProcessEquations[Flatten[{pde, ?}], u, Sequence @@ {r}]; sd = ndstate["SolutionData"][[1]]; feData = ndstate["FiniteElementData"]; pdeData = feData["PDECoefficientData"]; bcData = feData["BoundaryConditionData"]; methodData = feData["FEMMethodData"]; {DiscretizePDE[pdeData, methodData, sd], DiscretizeBoundaryConditions[bcData, methodData, sd], sd, methodData} ]  Next is the problem setup: ? = 10^-3; ? = 1; l = 2.2; h = 0.41; ? = RegionDifference[Rectangle[{0, 0}, {l, h}], ImplicitRegion[(x - 1/5)^2 + (y - 1/5)^2 < (1/20)^2, {x, y}]]; RegionPlot[?, AspectRatio -> Automatic]  ? = { DirichletCondition[p[x, y] == 0., x == l], DirichletCondition[{u[x, y] == 4*0.3*y*(h - y)/h^2, v[x, y] == 0}, x == 0], DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, y == 0 || y == h || (x - 1/5)^2 + (y - 1/5)^2 <= (1/20)^2]}; stokes = { D[u[x, y], x] + D[v[x, y], y], Div[{{-?, 0}, {0, -?}}.Grad[u[x, y], {x, y}], {x, y}] + D[p[x, y], x], Div[{{-?, 0}, {0, -?}}.Grad[v[x, y], {x, y}], {x, y}] + D[p[x, y], y] };  First we generate the system matrices for the Stokes equation: {dPDE, dBC, sd, md} = PDEtoMatrix[{stokes == {0, 0, 0}, ?}, {p, u, v}, {x, y} ? ?, Method -> {"FiniteElement", "InterpolationOrder" -> {p -> 1, u -> 2, v -> 2}, "MeshOptions" -> {"ImproveBoundaryPosition" -> False}}]; linearLoad = dPDE["LoadVector"]; linearStiffness = dPDE["StiffnessMatrix"]; vd = md["VariableData"]; offsets = md["IncidentOffsets"];  You could solve this stationary case, but we move on: The tricky part for non-linear equations is the linearization. For that I am referring you to Chapter 5. uOld = ConstantArray[{0.}, md["DegreesOfFreedom"]]; mesh2 = md["ElementMesh"]; mesh1 = MeshOrderAlteration[mesh2, 1]; ClearAll[rhs] rhs[t_?NumericQ, ut_] := Module[{uOld}, uOld = ut; Do[ ClearAll[u0, v0, p0]; (* create pressure and velocity interpolations *) p0 = ElementMeshInterpolation[{mesh1}, uOld[[offsets[[1]] + 1 ;; offsets[[2]]]]]; u0 = ElementMeshInterpolation[{mesh2}, uOld[[offsets[[2]] + 1 ;; offsets[[3]]]]]; v0 = ElementMeshInterpolation[{mesh2}, uOld[[offsets[[3]] + 1 ;; offsets[[4]]]]]; (* these are the linearized coefficients *) nlPdeCoeff = InitializePDECoefficients[vd, sd, "LoadCoefficients" -> {(* F *) {-(D[u0[x, y], x] + D[v0[x, y], y])}, {-? (u0[x, y]*D[u0[x, y], x] + v0[x, y]*D[u0[x, y], y]) - D[p0[x, y], x]}, {-? (u0[x, y]*D[v0[x, y], x] + v0[x, y]*D[v0[x, y], y]) - D[p0[x, y], y]} }, "LoadDerivativeCoefficients" -> -{(* gamma *) {{0, 0}}, {{? D[u0[x, y], x], ? D[u0[x, y], y]}}, {{? D[v0[x, y], x], ? D[v0[x, y], y]}} }, "ConvectionCoefficients" -> {(*beta*) {{{0, 0}}, {{0, 0}}, {{0, 0}}}, {{{0, 0}}, {{? u0[x, y], ? v0[x, y]}}, {{0, 0}}}, {{{0, 0}}, {{0, 0}}, {{? u0[x, y], ? v0[x, y]}}} }, "ReactionCoefficients" -> {(* a *) {0, 0, 0}, {0, ? D[u0[x, y], x], ? D[u0[x, y], y]}, {0, ? D[v0[x, y], x], ? D[v0[x, y], y]} } ]; nlsys = DiscretizePDE[nlPdeCoeff, md, sd]; nlLoad = nlsys["LoadVector"]; nlStiffness = nlsys["StiffnessMatrix"]; ns = nlStiffness + linearStiffness; nl = nlLoad + linearLoad; DeployBoundaryConditions[{nl, ns}, dBC]; diriPos = dBC["DirichletRows"]; nl[[ diriPos ]] = nl[[ diriPos ]] - uOld[[diriPos]]; dU = LinearSolve[ns, nl]; Print[ i, " Residual: ", Norm[nl, Infinity], " Correction: ", Norm[ dU, Infinity ]]; uOld = uOld + dU; (*If[Norm[ dU, Infinity ]<10^-6,Break[]];*) , {i, 8} ]; uOld ]  You'd then run this: uNew = rhs[0, uOld]; 1 Residual: 0.3 Correction: 0.387424 2 Residual: 0.000752321 Correction: 0.184443 3 Residual: 0.00023243 Correction: 0.0368286 4 Residual: 0.0000100488 Correction: 0.00264305 5 Residual: 3.6416*10^-8 Correction: 0.0000115344 6 Residual: 8.88314*10^-13 Correction: 1.22413*10^-10 7 Residual: 1.50704*10^-17 Correction: 1.08287*10^-15 8 Residual: 1.24246*10^-17 Correction: 6.93036*10^-16  See that the residual and correction converge. And do some post processing: p0 = ElementMeshInterpolation[{mesh1}, uNew[[offsets[[1]] + 1 ;; offsets[[2]]]]]; u0 = ElementMeshInterpolation[{mesh2}, uNew[[offsets[[2]] + 1 ;; offsets[[3]]]]]; v0 = ElementMeshInterpolation[{mesh2}, uNew[[offsets[[3]] + 1 ;; offsets[[4]]]]]; ContourPlot[u0[x, y], {x, y} ? mesh2, AspectRatio -> Automatic, PlotRange -> All, ColorFunction -> ColorData["TemperatureMap"], Contours -> 10, ImageSize -> Large] ContourPlot[v0[x, y], {x, y} ? mesh2, AspectRatio -> Automatic, PlotRange -> All, ColorFunction -> ColorData["TemperatureMap"], Contours -> 10, ImageSize -> Large] ContourPlot[p0[x, y], {x, y} ? mesh1, AspectRatio -> Automatic, PlotRange -> All, ColorFunction -> ColorData["TemperatureMap"], Contours -> 10, ImageSize -> Large] ` Which show the x-, y-velocity components and the pressure. This is a re-post of the original made with authors permission.
2 Replies
Sort By:
Posted 8 years ago
 Could you, please, suggest how to solve non-stationary problem in the similar manner?
Posted 9 years ago
 Please take a look also at the analogical post: Solving 2D Incompressible Flows using Finite Elements