Finite Elements For The (Navier) Stokes Equations

Download as pdf or txt
Download as pdf or txt
You are on page 1of 61

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

Noon, 08 December 2011

1 / 61

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

2 / 61

INTRO: Equations of Fluid Motion

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

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

4 / 61

EQUATIONS: The Navier Stokes Equations


The Navier-Stokes equations are the standard for uid motion. Any discussion of uid ow starts with these equations, and either adds complications such as compressibility or temperature, makes simplications such as time independence, or replaces some term in an attempt to better model turbulence or other features. vt v + (v t + )v + p =f (momentum equations)

(v) =0 (continuity equation)

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

EQUATIONS: Unsteady Compressible Navier Stokes

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

EQUATIONS: Steady, Incompressible Navier-Stokes

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

EQUATIONS: The Reynolds Number


If we are solving the full Navier Stokes equations, then the relative magnitude of measures the balance between diusion (smoothing) and nonlinear momentum terms (disruptive). As goes to zero, the character of the physical system, the PDE, and the discretized computer model deteriorate. Smooth solutions become irregular. Laminar ows become urbulent. Computationally, the nonlinear equations become dicult to solve. A time-dependent problem becomes unstable. The Reynolds number is a dimensionless quantity that estimates the dominance of momentum over diusion: Re = ||v ||L ||v ||L =

where L is a characteristic length. Mathematicians tend to solve problems with Re=1.


9 / 61

EQUATIONS: Steady, Incompressible Stokes


If the viscosity is large enough, the nonlinear terms can be neglected, and we have the steady Stokes equations: 2u 2u + )+ x 2 y 2 2v 2v ( 2 + 2 ) + x y u + x ( p =0 x p =0 y v =0 y

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

EQUATIONS: A Problem to Solve


We make take our equations of state to be either the Navier-Stokes or the Stokes equations. We assume we have the value of the kinematic viscosity . We assume we have been given information about a domain , within which the state equations hold. Moreover, we have information about , the boundary of , along which boundary conditions have been specied. Typically, these conditions include walls, inlets and outlets at which the velocity or some component of it is specied. We also need the value of pressure at one point, since it is a potential function, and thus unique only up to an additive constant. Together, this constitutes the mathematical model of the problem.
11 / 61

EQUATIONS: A Test Problem


Here is an example problem for us to solve, a rectangular channel with a square obstacle. Top and bottom are walls, ow enters from the left and exits on the right. Ive already created a grid using mesh2d:

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

EQUATIONS: MESH2D Can Make a Mesh For Us


Darren Engwirdas mesh2d can set up this grid: v = [ 0.0, -1.0; 8.0, -1.0; 8.0, +1.0; 0.0, +1.0; 1.5, -0.5; 1.5, +0.5; 2.5, +0.5; 2.5, -0.5 ]; <-- vertices e = [ 1, 2; 2, 3; 3, 4; 4, 1; 5, 6; 6, 7; 7, 8; 8, 5 ]; <-- connect vertex pairs to form boundaries hdata = []; hdata.hmax = 0.25;

<-- Maximum element size

[ 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

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

14 / 61

FEM: The LBB Condition


You may have heard that, when applying the nite element method to the Navier-Stokes equations for velocity and pressure, you cannot arbitrarily pick the basis functions. The interaction between the momentum and continuity equations can cause a stability problem; an unwary programmer can try to do everything right, and end up computing garbage. The problem that is going on is related to the inf-sup or Ladyzhenskaya-Babuska-Brezzi condition (LBB). Everyone in nite element uid calculations has their favorite way of avoiding the problem. The way we will do it is to use a Taylor-Hood pair of basis functions for the pressure and velocity. In a typical Taylor-Hood scheme, the polynomial degree of the pressure basis functions is one lower than that used for velocities.
15 / 61

FEM: A Linear Grid of Triangles

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

FEM: Pressure Grid

17 / 61

FEM: Pressure Representation

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

FEM: A Linear Basis Function

http://people.sc.fsu.edu/jburkardt/m src/fem2d basis t3 display/fem2d basis t3 display.html

19 / 61

FEM: A Quadratic Grid


Now we will construct a second grid that is a sort of renement of the rst. The set of nodes in this grid will be called vnodes, because they will be associated with velocities. We start by including all the pnodes, but we create a new node at the midpoint of every element edge, and add all these nodes as well. We can look at this procedure as involving two grids, one for pressure and one for velocities. The two grids are nested in an interesting way. The velocities will live on a grid of six-node triangles. These triangles share their vertices with the three-node pressure triangles. But the six-node triangles can be used to dene basis functions i (x, y ) which are 1 at node i, zero at all other nodes, and a quadratic function over each element.
20 / 61

FEM: Velocity Grid

21 / 61

FEM: Velocity Representation

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

FEM: A Quadratic Basis Function

This midside node basis function extends over two elements.


http://people.sc.fsu.edu/jburkardt/m src/fem2d basis t6 display/fem2d basis t6 display.html

23 / 61

FEM: A Quadratic Basis Function

This vertex basis function extends over six elements.


http://people.sc.fsu.edu/jburkardt/m src/fem2d basis t6 display/fem2d basis t6 display.html 24 / 61

FEM: Multiply by Test Functions


We have represented u, v and p in terms of basis functions. To try to determine the coecients in these representations, we multiply the state equations by the appropriate test functions: 2u 2u p ) i =0 + 2) + 2 x y x 2v 2v p (( 2 + 2 ) + ) i =0 x y y u v ( + ) i =0 x y ((

25 / 61

FEM: Integrate Over the Region

We integrate each equation over the region : p 2u 2u + 2) + ) i dx dy =0 2 x y x 2v 2v p (( 2 + 2 ) + ) i dx dy =0 x y y u v ( + ) i dx dy =0 y x ((

26 / 61

FEM: Integrate By Parts / Greens Theorem


We seek to lower the order of dierentiation on u and v : u i u i p + )+ i dx dy = x x y y x v i v i p ( + )+ i dx dy = x x y y y u v ( + ) i dx dy =0 y x ( u i ds n v i ds n

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

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

28 / 61

MAP: Reference Triangle Physical Triangle


The nite element equations are now a linear system of equations. The coecients are dened by integrals involving basis functions and their derivatives. We approximate the integrals using a quadrature rule. The integral approximations can be carried out one element at a time, so we can focus on the problem of estimating an integral over an arbitrary physical triangle. You may have seen approaches in which we approximate the integral by mapping each (x, y ) physical triangle to a single reference triangle (, ), where our quadrature rule is dened. There are complications in this method, especially when derivatives occur in our integrand. I will look at going the other way!

29 / 61

MAP: Quadrature in the Physical Triangle


We adjust the quadrature rule for a physical element T , with vertices (x1 , y1 ), (x2 , y2 ), and (x3 , y3 ). A reference abscissa (, ) is transformed to a physical abscissa (x, y ): x =x1 + x2 + (1 )x3 y =y1 + y2 + (1 )y3 A reference weight becomes a physical weight w : w = ( x1 (y2 y3 ) + x2 (y3 y1 ) + x3 (y1 y2 ) ) which simply multiplies the old weight by the area of T . The integral I(f,T) is approximated by Q(f,T) I (f , T ) Q(f , T ) =
i

wi f (xi , yi )

Derivatives in the integrand dont need any special treatment.


30 / 61

MAP: The Mapping Function


The mapping from reference triangle to physical triangle is: x =x1 + x2 + (1 )x3 y =y1 + y2 + (1 )y3 but this can be written as x y =A + x3 y3

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

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

32 / 61

BASIS: Hard Implementation, Simple Idea


You dont need to memorize a formula for basis functions, but you need to know theres a reasoning behind such formulas. You never want to look at code like the following and say I have no way of understanding this, so Ill believe it and use it til it breaks!
subroutine qbf (x,y,it,in,bb,bx,by,nelemn,nnodes,node,np,xc,yc) integer node(nelemn,nnodes) real xc(np), yc(np) if (in <= 3) then in1 = in; in2 = mod(in,3)+1; in3 = mod(in+1,3)+1 i1 = node(it,in1); i2 = node(it,in2); i3 = node(it,in3) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) t = 1.0+((yc(i2)-yc(i3))*(x-xc(i1))+(xc(i3)-xc(i2))*(y-yc(i1)))/d bb = t*(2.0D+00*t-1.0D+00) bx = (yc(i2)-yc(i3))*(4.0D+00*t-1.0D+00)/d; by = (xc(i3)-xc(i2))*(4.0D+00*t-1.0D+00)/d else inn = in-3; in1 = inn; in2 = mod(inn,3)+1; in3 = mod(inn+1,3)+1 i1 = node(it,in1); i2 = node(it,in2); i3 = node(it,in3); j1 = i2; j2 = i3; j3 = i1 d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) c = (xc(j2)-xc(j1))*(yc(j3)-yc(j1))-(xc(j3)-xc(j1))*(yc(j2)-yc(j1)) t = 1.0D+00+((yc(i2)-yc(i3))*(x-xc(i1))+(xc(i3)-xc(i2))*(y-yc(i1)))/d s = 1.0D+00+((yc(j2)-yc(j3))*(x-xc(j1))+(xc(j3)-xc(j2))*(y-yc(j1)))/c bb = 4.0D+00*s*t bx = 4.0D+00*(t*(yc(j2)-yc(j3))/c+s*(yc(i2)-yc(i3))/d) by = 4.0D+00*(t*(xc(j3)-xc(j2))/c+s*(xc(i3)-xc(i2))/d) end if return end

http://people.sc.fsu.edu/jburkardt/f src/channel/channel.html 33 / 61

BASIS: The Linear Basis Functions


The basis functions for pressure are dened on the three vertex triangle T = {(x1 , y1 ), (x2 , y2 ), (x3 , y3 )}. Basis 1 (x, y ) is 1 at vertex 1, 0 at the other two vertices, and linear over T . Rather than looking up a formula, can we work one out? If 1 (x, y ) is linear, and its zero at nodes 2 and 3, then its zero on the line between them. The slope of the line through (x2 , y2 ) is: y3 y2 s(x3 , y3 ) = x3 x2 and for an arbitrary point (x, y ), the slope is: y y2 s(x, y ) = x x2 We want 1 (x, y ) to be zero if s(x, y ) = s(x3 , y3 ). Lets try y y2 y3 y2 ? 1 (x, y ) = s(x, y ) s(x3 , y3 ) = x x2 x3 x2
34 / 61

BASIS: The Linear Basis Functions


Lets avoid fractions by multiplying through by the denominators: 1 (x, y ) = (y y2 )(x3 x2 ) (y3 y2 )(x x2 ) Notice that 1 (x2, y 2) = 1 (x3, y 3) = 0. What more do we need? Oh yes, we need that 1 (x1, y 1) = 1 Easy! Just normalize this function by its value at (x1 , y1 ): 1 (x, y ) = (y y2 )(x3 x2 ) (y3 y2 )(x x2 ) (y 1 y2 )(x3 x2 ) (y3 y2 )(x1 x2 )
?

Since 1, 2, 3 are arbitrary, we also dened 2 (x, y ) and 3 (x, y )!

35 / 61

BASIS: The Quadratic Basis Functions


Lets symbolize the six node triangle this way: N1 / \ / \ N4 N6 / \ / \ N2----N5----N3 Just as for the linear basis functions, we can nd a linear function which is zero along any line we choose. Therefore, there is a linear function that is zero at N2 and N1 (and hence at N12 as well). Another linear function is zero at N23 and N31, and so on. Will this help us nd a quadratic basis function?
36 / 61

BASIS: The Quadratic Basis Functions


Suppose we want to nd 3 ? There is a linear function g (x, y ) that is zero at N1, N4, and N2. There is a linear function h(x, y ) that is zero at N5 and N6. Therefore, what about 3 (x, y ) = g (x, y ) h(x, y ) Almost, but we need it to be 1 at (x3 , y3 ). Easy again: 3 (x, y ) = g (x, y ) h(x, y ) g (x3 , y3 ) h(x3 , y3 )
?

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

BASIS: A Quintic Basis For Triangles


If we need a 5-th degree polynomial basis, we take a reference triangle and make 6 rows of dots in each direction.

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

BASIS: A Quintic Basis Function


Can we describe the basis function at the blue node?

(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

BASIS: The Quadratic Basis Functions


I dont necessarily want you to ever have to work out the arithmetic involved in computing a quadratic basis function...or its derivatives with respect to x and y . But I do want to convince you that its not magic, it doesnt require a special course in analysis, its really just some carefully thought-out high school geometry! You should be able to see how to construct: cubic elements in a 2D triangle; youll need 10 nodes; quadratic elements in a 3D tetrahedron; again you need 10 nodes; instead of eliminating nodes on 2 lines, you look for 2 planes. There are three cases to consider, a vertex, mid-edge node, or mid-face node. quintic elements in a 4D simplex - see how easy it is?
http://people.sc.fsu.edu/jburkardt/m src/fem basis/fem basis.html 40 / 61

FEM NAVIER STOKES

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

ASSEMBLY: The Pressure Equation


We are looking at node I. Node I has a pressure equation associated with it only if it is actually a vertex of the triangle, what we called a pnode. Lets assume that this is the case. The pressure equation (continuity equation) is (

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

ASSEMBLY: The Horizontal Velocity Equation


Ignoring boundary terms, the horizontal velocity equation is: (

u i u i p + )+ i dx dy = 0 x x y y x

So we always get a contribution to A(I,J): A(I , J) = A(I , J) +

j i j i + ) dx dy x x y y

and if J is a pressure node, we get a contribution to C(I,J): C (I , J) = C (I , J) +

j i dx dy x

Do you see that were just dierentiating the equations with respect to a coecient?
44 / 61

ASSEMBLY: MATLAB Code


% % % Add terms to the horizonal momentum equation. a(iu,ju) = a(iu,ju) + w(quad) * nu ... * ( dbidx(test) * dbjdx(basis) + dbidy(test) * dbjdy(basis) ); if ( 0 < jp ) a(iu,jp) = a(iu,jp) + w(quad) * bi(test) * dqjdx(basis); end % % % Add terms to the vertical momentum equation. a(iv,jv) = a(iv,jv) + w(quad) * nu ... * ( dbidx(test) * dbjdx(basis) + dbidy(test) * dbjdy(basis) ); if ( 0 < jp ) a(iv,jp) = a(iv,jp) + w(quad) * bi(test) * dqjdy(basis); end % % % Add terms to the continuity equation. if ( 0 < ip ) a(ip,ju) = a(ip,ju) + w(quad) * qi(test) * dbjdx(basis); a(ip,jv) = a(ip,jv) + w(quad) * qi(test) * dbjdy(basis); end

http://people.sc.fsu.edu/jburkardt/m src/fem2d stokes sparse/fem2d stokes sparse.m 45 / 61

ASSEMBLY: It All Ends Up as a Linear System


Of course, we dont have separate matrices called A, B and C, so we have to store all these coecients in one big matrix, and we store the coecients of the representations for u, v and p in one big vector. Because we have multiple equations and variables, and a pair of grids, a lot of the programming involves simply guring out where to put things and how to get them back! We still have some boundary conditions to take care of, but thats another side issue. In the end, we wind up with a sparse linear system: Ax =b that we solve for the nite element coecients that give us functional representations of the state variables.
46 / 61

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

47 / 61

IFISS: A MATLAB PDE Solver


Its important to nd and use good software tools that other people have written. It helps you to start solving interesting problems right away, it lets you see how someone has worked out the solution of the underlying software issues, and it gives you a good base from which to add new software features for your own research. IFISS = Incompressible Flow Iterative Solution Solver is a MATLAB package that is a very useful tool for people interested in learning about solving PDEs.

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

IFISS: Customized Problems


IFISS comes with sample problems, which can guide the user in designing a new problem. The domain, and its gridding, are dened by a function such as grids/myow domain.m. The user supplies lists of: vertices that outline the boundary and internal holes; boundary edges for Dirichlet or Neumann conditions; obstacles (v1, v2, ..., vn); Boundary conditions and sources are specied by: myow bc(x,y) returns specied stream function values; myow ow(x,y) returns specied ow values;

50 / 61

IFISS: Customized Problems


The user can dene the PDEs to be solved as well. Actually, this means writing the code to assemble the system matrix. Here is part of the code for the Stokes equations, which should start to look familiar now! The INVJAC and JAC factors arise because these equations are integrated in the reference element, not in their home element.
for j = 1:9 for i = 1:9 ae(:,i,j) = ae(:,i,j) + ae(:,i,j) = ae(:,i,j) + re(:,i,j) = re(:,i,j) + bbxe(:,i,j) = bbxe(:,i,j) bbye(:,i,j) = bbye(:,i,j) end

wght*dpsidx(:,i).*dpsidx(:,j).*invjac(:); wght*dpsidy(:,i).*dpsidy(:,j).*invjac(:); wght*psi(:,i).*psi(:,j).*jac(:); - wght*psi(:,i) .*dpsidx(:,j); - wght*psi(:,i) .*dpsidy(:,j);

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

IFISS: Channel Flow With Obstacle


Here is an IFISS grid for problem NS5. Top and bottom are walls, ow enters from the left and leaves on the right, and theres a square obstacle.

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

IFISS: Channel Flow With Obstacle


Here is how to run IFISS with default data for the obstacle problem: >> setpath >> navier_testproblem <-- sets up MATLAB path <-- request a Navier-Stokes test

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

IFISS: Channel Flow With Obstacle (More Choices)


Now we set the grid size and shape, the velocity and pressure basis functions, and the viscosity: Grid generation for domain with obstacle. grid parameter: 3 for underlying 8x20 grid (default is 4) : return uniform/stretched grid (1/2) (default is uniform) : return Q1-Q1/Q1-P0/Q2-Q1/Q2-P1: 1/2/3/4? (default Q1-P0) : return setting up Q1-P0 matrices... done system matrices saved in obstacle_stokes_nobc.mat ... Incompressible flow problem on obstacle domain ... viscosity parameter (default 1/50) : return
54 / 61

IFISS: Channel Flow With Obstacle (More Choices)


Now we specify some solver options. Picard/Newton/hybrid linearization 1/2/3 (default hybrid) : return number of Picard iterations (default 6) : return number of Newton iterations (default 5) : return nonlinear tolerance (default 1.d-8) : return stokes system ... Stokes stabilization parameter (default is 1/4) : return setting up Q1 convection matrix... done. uniform/exponential streamlines 1/2 (default uniform) : return number of contour lines (default 50) : return
55 / 61

IFISS: Stokes Flow


IFISS displays the Stokes ow used for initialization.

56 / 61

IFISS: Navier Stokes Flow


The nal Navier-Stokes solution shows signicant dierences.

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

FEM NAVIER STOKES

Introduction Equations of Fluid Motion A Finite Element Formulation The Mapping Function Computing Basis Functions Assembling the Matrix IFISS Conclusion

59 / 61

CONCLUSION: The Big Picture


Ive shown you that the nite element method can be applied to the Navier-Stokes equations, and that if you are interested in doing so there are a lot of choices constraints techniques: tools to be familiar with. I hope, at least, Ive given you an idea of the things you can expect if you are interested in pursuing computations involving uid ow. You should be able at least to consider how to deal with some of the problems I have not mentioned, such as boundary conditions, adding temperature eects, or solving the nonlinear system associated with the Navier-Stokes equations.
60 / 61

CONCLUSION: Some References


Howard Elman, Alison Ramage, David Silvester, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford, 2005, ISBN: 978-0198528678, LC: QA911.E39. Howard Elman, Alison Ramage, David Silvester, IFISS, A Matlab Toolbox for Modeling Incompressible Flow, ACM Transactions on Mathematical Software, Volume 33, Number 2, June 2007, Article 14. Anders Logg, Garth Wells, DOLFIN: Automated Finite Element Computing, ACM Transactions on Mathematical Software, Volume 37, Number 2, April 2010, Article 20. Per-Olof Persson, Gilbert Strang, A Simple Mesh Generator in MATLAB, SIAM Review, Volume 46, Number 2, June 2004, pages 329-345.

61 / 61

You might also like