Finite Elements For The (Navier) Stokes Equations
Finite Elements For The (Navier) Stokes Equations
Finite Elements For The (Navier) Stokes Equations
Numerical Analysis Seminar .......... John Burkardt Department of Scientic Computing Florida State University http://people.sc.fsu.edu/jburkardt/presentations/pitt 2011.pdf
1 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
2 / 61
Your introduction to the nite element method probably started with the Poisson equation on a rectangle; its important to see how those ideas can be extended to harder problems in a more general geometry. Well look at the Navier-Stokes equations for uid ow; even in their simple form, there are some unexpected features. Well consider the question of what terms to keep, how to dene a mesh, how to set up the basis, and to assemble the matrix. Well nish with a look at IFISS, a MATLAB program that makes it easy to do some sophisticated computations for PDEs in 2D.
3 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
4 / 61
v is the velocity vector; p is the pressure; is the uid density; is the dynamic viscosity; f represents body forces such as gravity.
5 / 61
Here are the time-dependent compressible Navier Stokes equations, in 2D Cartesian coordinates: 2u 2u u p u u ( 2 + 2 ) + u + v + =0 t x y x y x v 2v v p 2v v ( 2 + 2 ) + u + v + = g t x y x y y u v + + =0 t x y
6 / 61
EQUATIONS: Simplications
We might be interested in steady state ow, in which case we can drop the time derivatives. If we assume that density is constant, we let the velocity absorb a multiple of , call it mass velocity, and still represent it as (u, v ). The pressure p can absorb the gravity force, and can be rescaled by the constant density, and we still call it p. The dynamic viscosity can be rescaled by density to yield the kinematic viscosity .
7 / 61
Now we have the steady incompressible Navier Stokes equations: 2u 2u u p u +v + =0 + 2) + u 2 x y x y x 2v 2v v v p ( 2 + 2 ) + u +v + =0 x y x y y u v + =0 x y ( The viscosity multiplies the nice Poisson operator, which represents the tendency of momentum to spread out or diuse.
8 / 61
Because these equations are linear in u, v and p, they are much easier to work with. They are also useful even if we want to solve the Navier-Stokes equations, because they can give us a reasonable starting solution.
10 / 61
In a nite element approach, we divide the region into small elements that are only inuenced by their immediate neighbors. Typically, the solution is sought at the vertices of the elements. For our uid problem, our approach will be a little more complicated!
12 / 61
[ p, t ] = mesh2d ( v, e, hdata );
http://www.mathworks.com/matlabcentral/leexchange/25555-mesh2d-automatic-mesh-generation http://people.sc.fsu.edu/jburkardt/classes/fem 2011/fem meshing.pdf 13 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
14 / 61
Using a program like MESH2D, we can take a description of a region , perhaps outlined by a set of vertices, and produce a set of nodes which can be triangulated so that triplets of nodes dene elements. As you have probably seen before, such a triangulation allows us, in a natural way, to dene linear basis functions i (x, y ), which are 1 at node i, 0 at all other nodes, and linear over each element. For reasons I will explain in a minute, lets call the nodes weve just created pnodes. This name is meant to suggest that these nodes are associated with the pressure variable.
16 / 61
17 / 61
We need to represent our variables as linear combinations of basis functions. The easy case is the pressure p. We can take this to be a linear combination of piecewise linear basis functions i (x, y ),
pnodes
p=
i=1
ci i (x, y )
where the i-th basis function is associated with the i-th pnode.
18 / 61
19 / 61
21 / 61
Our velocities will similarly be represented using the quadratic functions. Since velocity is a vector, we can think of it as having components (u, v ). Our representation can then be written:
vnodes
u=
i=1 vnodes
ai i (x, y )
v=
i=1
bi i (x, y )
22 / 61
23 / 61
25 / 61
26 / 61
The right hand sides are only interesting (nonzero) for nodes on the boundary where a normal inow or outow condition is allowed. Right now, we dont need things to get any more interesting, so well assume the right hand sides are all zero!
27 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
28 / 61
29 / 61
wi f (xi , yi )
and so the inverse map is easy to construct: (y2 y3 ) (x x3 ) (x2 x3 ) (y y3 ) (x1 x3 ) (y2 y3 ) (y1 y3 ) (x2 x3 ) ((y1 y3 ) (x x3 ) + (x1 x3 ) (y y3 ) = (x1 x3 ) (y2 y3 ) (y1 y3 ) (x2 x3 ) = The denominator is a multiple of the area of the triangle. You can construct a two-way map between any pair of triangles.
31 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
32 / 61
http://people.sc.fsu.edu/jburkardt/f src/channel/channel.html 33 / 61
35 / 61
The product of two linear functions is, of course, quadratic. Pick any node on the six node triangle, and you can cover the other ve nodes with two straight lines. End of story!
37 / 61
For each node, we must dene a degree 5 polynomial in x and y that is 1 at that node, and zero at the other 20. A quintic function could be dened as the product of ve linear functions. A linear function covers up points on a line. Find ve straight lines that cover all the other nodes!
38 / 61
(x, y ) = (x) (y ) (y 0.2) (x + y 1) (x + y 0.8) This will be zero at all the red nodes. Of course, we have to normalize the function to get a value of 1 at the blue node. To nd the linear factors for any node, simply walk to each boundary, and note the parallel lines you cross.
39 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
41 / 61
ASSEMBLY
When its time to assemble the matrix, we have to keep in mind that we have three variables to worry about and two related grids. To assemble the equation associated with a variable at a given node, we have to consider all the elements that include that node, all the nodes in those elements, and all the variables associated with those nodes. You can see there can be a lot of bookkeeping! But at some point, were looking at the equation for node I, and considering contributions from variables dened at node J. These contributions get added to the (I , J) matrix element, and if we want to, we can call this element A(I,J) for horizontal velocity, B(I,J) for vertical velocity, and C(I,J) for pressure.
42 / 61
u v + ) i dx dy = 0 x y
So, oddly enough, pressure itself doesnt show up in the pressure equation! However, node J will contribute to coecients A(I,J) and B(I,J) for the horizontal and vertical velocities: A(I , J) = A(I , J) +
B(I , J) = B(I , J) +
j i dx dy x j i dx dy y
43 / 61
u i u i p + )+ i dx dy = 0 x x y y x
j i j i + ) dx dy x x y y
j i dx dy x
Do you see that were just dierentiating the equations with respect to a coecient?
44 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
47 / 61
http://www.cs.umd.edu/elman/iss.html
48 / 61
IFISS: Features
IFISS includes built-in solvers for 2D versions of: the the the the Poisson equation convection-diusion equation Stokes equations Navier-Stokes equations
The user can specify the geometry and the boundary conditions, and time dependence. The package uses MATLABs sparse storage structure; it can use MATLABs sparse direct solver, but also can invoke iterative solvers, including GMRES and multigrid methods. IFISS oers a variety of mixed nite element bases for ow: Stable rectangular: Q2 Q1 or Q2 P1 ; Stabilized rectangular: Q1 P0 or Q1 Q1 ;
49 / 61
50 / 61
for i=1:3 bxe(:,i,j) = bxe(:,i,j) - wght*chi(:,i) .* dpsidx(:,j); bye(:,i,j) = bye(:,i,j) - wght*chi(:,i) .* dpsidy(:,j); end end
51 / 61
Yes, IFISS uses quadrilateral elements, not triangles! The pressures are piecewise constant (asterisks at centers) and the velocities are piecewise linear (vertices of quadrilaterals).
52 / 61
specification of reference Navier-Stokes problem. choose specific example (default is cavity) 1 Channel domain 2 Flow over a backward facing step 3 Lid driven cavity 4 Flow over a plate 5 Flow over an obstacle : 5
53 / 61
56 / 61
57 / 61
IFISS: Overview
If youre interested in uid ow, IFISS is a great place to start. You can learn a lot just by looking at how it is put together. You can easily set up new 2D problems (new geometry, boundary conditions, source terms) You can also use it as a starting point for new algorithms you are interested in.
58 / 61
Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion
59 / 61
61 / 61