Introduction to CFd analysis.FEM4CFD-Week1.pptx

Introduction to CFD analysis . Grid genaration, application of CFD in various disciplines.

1. A short course on FEM for CFD Week 1 : Basics of FEM and FreeFEM , and Poisson equation Dr Chennakesava Kadapa Email: [email protected] 1

2. Contents Introduction to the course Poisson’s equation Navier-Stokes equation Introduction to Finite Element Method Introduction to the FreeFEM library Examples 2

3. Delivery mode Lecture Tutorial Exercise 3

4. 4 References More references at https://github.com/chennachaos/awesome-FEM4CFD

5. Incompressible Navier-Stokes equation   Acceleration term Viscous term Pressure gradient Local Acceleration Convective Acceleration Body force 5

6. Course structure Poisson’s equation Steady Advection-Diffusion equation Unsteady Advection-Diffusion equation Stokes equation Navier-Stokes equation   6 Week 1 Week 2 Week 3 Week 4,5,6 Week 4,5,6

7. Poisson’s equation 7

8. Poisson’s equation Is an elliptic PDE Many applications Diffusion, heat transfer, electrostatics, magnetostatics Poisson’s equation is where is the (scalar) field, is the source term and is the Laplacian operator. Laplacian operator ( )   8

9. Poisson’s equation (Cont’d) We can also write as where is the gradient operator and is the divergence operator. Gradient operator ( ) Thus   9  

10. Boundary conditions Boundary conditions Dirichlet boundary condition: Neumann boundary condition: where is the specified solution field on the Dirichlet boundary is the specified flux on the Neumann boundary .   10

11. Finite Element Method 11

12. Finite Element Method (FEM) is a numerical method for solving problems modelled by partial differential equations (PDEs). FEM is not just for solid mechanics. FEM has become the de-facto numerical method for many problems in science and engineering. Simulation software and/or libraries based on FEM are available for Solid mechanics – stress analysis Structural dynamics Wave propagation Heat transfer Fluid dynamics Chemical diffusion Electromagnetics 12 What is FEM? Check my papers for the details.

13. FEM offers a flexible numerical framework with a rich set of elements of different shapes, orders and families. Applicability to complex geometries. Different element shapes – triangle, quadrilateral, tetrahedron, brick, pyramid, wedge. Different element families – Lagrange family, B é zier family, Spectral elements, IGA, H- Div and H-Curl elements. Ease of extensions to higher-order approximations Ease of extension to Multiphysics problems Better accuracies compared to FVM for unstructured meshes 13 Why FEM?

14. Dates back to 1940s Originally used for problems in structural mechanics 1960-1990 was the golden age of FEM Since 1990, an explosion of industrial applications of FEM 14 History of FEM https://link.springer.com/article/10.1007/s11831-022-09740-9

15. The points are called nodes . The cells are called elements . The variables are usually stored at nodes. Variable per cell is also possible. The solution variables at nodes are called as degrees of freedom (DOFs) . The value at nodes (DOFs) can be a Scalar: temperature, pressure, density Vector: displacement, velocity Tensor: deformation gradient, strain, stress 15 FEM terminology

16. Weak formulation 16

17. Weak form Poisson’s equation Weighted residual statement Splitting the integral After using the integrating by parts and divergence theorem Using boundary conditions   Residual:   Weak form 17

18. Weak form Weak form in compact form Weak form in component form   1D 2D 3D 18

19. Finite Element Method In the FEM, we split the domain into smaller parts, called elements . As the size of the elements is reduced, the approximation error is reduced. Elements are described by their connections to nodes. Discretisation   ● is a node. 19  

20. Discretised problem Weak form: Using the discretization, the system of equations can be written as where,   20

21. FEM ingredients We need to evaluate Information needed: Mesh – nodal coordinates and connectivity Element – basis functions Stiffness matrix and force vector - numerical quadrature Linear system – matrix solver (LU, CG, etc..)   21

22. FreeFEM library 22

23. Introduction to FreeFEM 23 An opensource library for FEM. Not an application/solver library like OpenFOAM . Library that provides all the basic FEM ingredients. To develop custom solvers. Well supported. Easy to install and use. Up-to-date tutorials. Supports triangular and tetrahedral elements only. Input through scripts File extension “ . edp ” Similar to C++ syntax No need to compile FreeFem++.exe < inp -script-file>

24. Tutorial #1 24

25. Tutorial #1 25 File : Week1-Tutorial1.edp

26. Tutorial #1 // The triangulated domain Th mesh Th = square (10,10); // The finite element space defined over Th is called here Vh fespace Vh ( Th , P1 ); // Define u and w in the function space Vh u , w ; // Define functions func uexact = ( cosh (pi*y)- sinh (pi*y)/tanh(pi))*sin(pi*x); func f = ; Mesh over a square or mapped square Finite element space. Element type, order etc. P1 is 3-noded triangle. Solution and test function Functions for use 26

27. Tutorial #1 (cont’d) // Define the weak formulation solve Poisson ( u , w, solver = LU ) = int2d ( Th )( // The bilinear part dx (w) * dx (u) + dy (w) * dy (u) ) - int2d ( Th )( // The right hand side w * f ) // The Dirichlet boundary condition + on ( 1,2,3,4 , u = uexact ) ; // Plot the result plot ( u ); load " iovtk " int [ int ] Order = [1]; string filename = "Poisson2D-square-Dirichlet.vtu" ; savevtk ( filename, Th, u, dataname = "u" , order=Order); Weak form LHS Weak form RHS Boundary conditions Plot the solution Export solution in VTK format. Available only as first-order elements . 27  

28. Tutorial #1 (cont’d) // Error real error = int2d(Th)( ( uexact -u)^2 ); error = sqrt(error); 28 Error L2 norm of error  

29. Exercise #1 29 Different number of elements // The triangulated domain Th mesh Th = square (10,10); mesh Th = square (20,20); mesh Th = square (80,80); Different orders of element // The finite element space defined over Th is called here Vh fespace Vh ( Th , P1 ); fespace Vh ( Th , P2 ); load "Element_P3"; fespace Vh ( Th , P3) ; Solve the problem in Tutorial 1 using different number of elements and different orders of elements and observe the error.

30. Neumann BCs 30

31. Tutorial #2 31  

32. Tutorial #2 (cont’d) // Define the weak form solve Poisson(u, w, solver=LU)     = int2d(Th)(    // The bilinear part           dx(w)*dx(u) + dy (w)* dy (u)     )     - int2d(Th)( // Source term           w*f     )     - int1d(Th,3)( // Neumann BC term           w* tbar     )     + on(1,2,4, u=20)  // Dirichlet BC     ; 32  

33. Complex geometries 33

34. “square” function // The triangulated domain Th // unit square , starting at (0,0) mesh Th = square(10, 1 0); // rectangle of size 10x5, starting at (0,0) mesh Th = square(10, 1 0, [10*x, 5*y]); // rectangle of size 5x10, starting at (-4,-2) mesh Th = square(10, 1 0, [-4+5*x, -2+10*y]); 34 //display mesh plot(Th);

35. Complex geometries – inbuilt functions int boundary = 1; border C01(t=0, 1){x=t; y=0.0; label=boundary;} border C02(t=0, 1){x=1; y=0.5*t; label=boundary;} border C03(t=0, 1){x=1.0-0.5*t; y=0.5; label=boundary;} border C04(t=0, 1){x=0.5; y=0.5+0.5*t; label=boundary;} border C05(t=0, 1){x=0.5-0.5*t; y=1.0; label=boundary;} border C06(t=0, 1){x=0.0; y=1-t; label=boundary;} int n = 20; mesh Th = buildmesh (C01(n) + C02(n/2) + C03(n/2) + C04(n/2) + C05(n/2) + C06(n)); plot(Th); C01 C02 C03 C04 C05 C06 35 https://doc.freefem.org/examples/mesh-generation.html

36. Complex geometries – using GMSH load " gmsh " // The triangulated domain Th mesh Th = gmshload (" Lshaped-gmsh-mesh.msh "); int boundary = 1; mesh Th = gmshload ("cylinder- gmsh - mesh.msh "); int inlet = 1; int outlet = 2; int topwall = 3; int bottomwall = 4; int cylinder = 5; 36 //display mesh plot(Th);

37. Exercises 37

38. Exercise #2 Solve the Poisson equation over the L-shaped domain. Use the mesh from “ Lshaped-gmsh-mesh.msh ”. on all boundaries.     38

39. Solution for exercise #2 39

40. Exercise #3 Solve the Poisson equation for the following problem. Use the mesh from “ cylinder- gmsh - mesh.msh ”.             40

41. Solution for exercise #3 41

42. 42

43. Access to matrix and vector // Define the weak (variational) formulation varf PoissonLHS ( u , w ) = int2d ( Th )( // The bilinear part dx ( w) * dx ( u ) + dy ( w ) * dy ( u ) ) // The Dirichlet boundary condition + on ( 1,2,3,4 , u = uexact ) ; matrix A = PoissonLHS ( Vh , Vh ); //stiffness matrix cout << A << endl ; // Define the weak (variational) formulation varf PoissonRHS ( unused , w ) = int2d ( Th )( f * w ) // The Dirichlet boundary condition + on ( 1,2,3,4 , u = uexact ) ; Vh F ; //F[] is the vector associated to the function F F [] = PoissonRHS (0, Vh ); cout << F[] << endl ; //u[] is the vector associated to the function u u [] = A ^ -1 * F []; cout << u[] << endl ; 43

44. More practice with FreeFEM FreeFEM - Getting started https://doc.freefem.org/tutorials/poisson.html FreeFEM – Additional examples https://doc.freefem.org/tutorials/index.html 44

45. 45

46. 46