An Overview Of CFD Tools And Comparison Of A 2D Testcase

2. The Finite Volume and Finite Element Methodhence the velocity field is a solenoidal (divergence-free) vector field. As a consequence, thepressure is no longer coupled to density and temperature. As a result, and with the additionalassumption of constant viscosity, the energy equation and the continuity and momentum equations can be decoupled. Because the temperature of the fluid is not of concern here, the energyequation can be neglected.The assumption that the density is constant is reasonable because the fluid modeled by theNavier Stokes equations is a slow-moving liquid.The momentum equation can be rewritten (see e.g. [6]) as followsρ U · (ρUU) p µ 2 U ρ g F s, t(2.2) where F s represents the surface tension force F s σκnwith the surface tension coefficient σ, the unit vector normal to the interface n and theinterface curvature κ, the latter two are defined asn α α κ · n.In summary, isothermal and incompressible flow is governed by (2.1) and (2.2). They haveto be discretized to numerically calculate their solution. Two different ways to discretize theequations are presented in Sections 2.2 and 2.3.2.2. The Finite-Volume MethodThe finite-volume method (FVM) [7, 8] is a method for the numerical solution of partialdifferential equations. In the field of computational fluid dynamics, these equations describe8

2. The Finite Volume and Finite Element Methodconservation laws, i.e. the conservation of mass, momentum and energy. In general, a conservation law states that a certain property in a closed, isolated system remains constant overtime. An example is the generic equation φ · F (φ) s t(2.3)with a scalar quantity φ, a flow vector F (φ) and a source term s. If the source term is equalto zero and no flow occurs, φ does not change.The FVM requires two steps. The equation of interest is integrated over a control volume Ωthat is constant over time, and is then transformed by the divergence theorem into a volumeand a surface integral, which can then be discretized independently. Applying these steps to(2.3) yieldsZ tZ tZ · F (φ)dV φ dV Zs dV,(2.4)s dV.(2.5)ΩΩΩwhich is equivalent toIφ dV F (φ) · dAΩ ΩΩZNow the balance over a finite volume Vj can be taken, resulting inX sj Vj .(φj Vj ) F (φi ) · dA tcellf aces(2.6)For the discretization, a grid has to be generated. For the sake of simplicity, a two-dimensionalgrid with rectangular cells (cp. Figure 2.1) is assumed. A value φkj is assigned to the centerof each cell Vj for every point in time tk . Now the first part of (2.5) can be discretized usingfor example the forward difference quotient:ZVjφk 1 φkj (φ)jdV x y t tThe discretization of the flux is performed separately for the four cell faces. To do this, thevalue of F (φ) on the cell boundaries has to be approximated.9

2. The Finite Volume and Finite Element MethodFigure 2.1.: Discrete cells on a cartesian grid [9]The discretization method used should be suitable for the problem or equation. If the flux isgoverned by diffusion, central differences are the method of choice, because diffusion describesa balancing process that does not have a specific direction.If, on the other hand, the flux is dominated by convection, the flow direction plays an important role. In this case, an upwind discretization scheme should be applied [7, 10].Often, convection and diffusion occur simultaneously and have to be accounted for. In thecase of the convection-diffusion equation, the flow vector has the formF (φ) Uφ D φwith the diffusion coefficient D. The discretization of the flux will be demonstrated exemplarilyfor horizontal flow over the right (eastern) cell face. F (φ) reduces tof (φ) uφ D φ x(2.7)where u is the component of U in the direction of the x-axis. The value on the cell boundary isdesignated with the subscript e, the values at the cell centers to the left and right are denotedwith the subscripts P and E, respectively. Discretizing the diffusive part of (2.7) with centraldifferences yields10

2. The Finite Volume and Finite Element Methodfe ue φe DφE φP. xEPThe value of φe is still undefined. Because it is part of the convective part of the equation,instead of taking an average of the formφe φE φP,2φe is set to the value upstream of the boundary. This full upwind discretization renders thefollowing one sided formulation:(ue φ e ue φPue φEue 0ue 0Once fe , and analogous fw , fn and fs (the flow over the left, top and bottom cell boundaries)have been fully discretized, the surface integral in (2.5) can be approximated:I (fe fw ) y (gn gs ) xF (φi ) · dA ViIn practice, there are many other and more advanced discretization schemes. The ability tochoose different schemes that are suited to different equations is one advantage of the FVM.Another strength is its ability to work on unstructured meshes [8].The main steps of the FVM will now be applied to the continuity and momentum equations,because they are the equations of interest. For the continuity equation, this yields:Z ρ dV t ΩZ ρ dV t Ω ρ · (ρU) 0Z tΩ · (ρU)dV 0I 0ρU · dA ΩHere, ρU is called the mass flux. It can be seen that the change of mass within the volume Vover time is equal to the flow over the surface or boundary. Again, the balance over a finitevolume Vj is applied, giving11

2. The Finite Volume and Finite Element MethodX 0.(ρj Vj ) ρi Ui · dA tcellf acesApplying the same steps to the momentum equation results inX 0.(ρj Uj Vj ) ρi Ui · dA tcellf aces2.3. The Finite-Element MethodThe finite-element method rests upon the discrete representation of a weak integral form ofthe partial differential equation to be solved [10]. To derive the weak form, two classes offunctions are needed, namely the test functions and trial solutions.The test functions are all functions which are square integrable, have square integrable firstderivatives over the computational domain Ω, and vanish on the Dirichlet part ΓD of theboundary. They are defined byV {w H 1 (Ω) w 0 on ΓD },where H 1 (Ω) is a Sobolev space. The trial solutions are similar, except in that they arerequired to satisfy the Dirichlet conditions on ΓD :S {u H 1 (Ω) u uD on ΓD }.For simplicity, the main steps of the finite-element method will be demonstrated using thePoisson equation 2 u sin Ω,(2.8)with u uD on ΓD and n · u h on ΓN (the Neumann part of the boundary).The first step is to multiply the equation of interest with a test function w and to integrateover Ω:ZZ2 w u dΩ w s dΩ(2.9)ΩΩ12

2. The Finite Volume and Finite Element MethodThen the Green-Gauss divergence theorem is applied to the left hand side of (2.9):ZZ2 w u dΩ ( · (w u) w · u) dΩZ ΩZ w · u dΩ w(n · u)dΓΩΩΓThis changes the requirements on u, it now only has to be once differentiable, as has w.Inserting the Neumann boundary condition gives the weak form of (2.8):ZZZ w · u dΩ Ωw s dΩ Ωw h dΓ(2.10)ΓNObviously, any solution to the original problem (2.8) is a solution to the weak problem. It canbe shown by means of the Lax-Milgram lemma [10] that a solution to the weak problem isunique and is therefore identical to the strong solution, if the latter exists.Now (2.10) can be discretized by partitioning the computational domain Ω into a finite numberof subdomains, or elements.Ω [eΩΩe Ωf for e 6 fandSimilarly, finite dimensional subspaces of the test function and trial solution spaces are defined,V h : {w H 1 (Ω) w Ωe Pm (Ωe ) e and w 0 on ΓD }S h : {u H 1 (Ω) u Ωe Pm (Ωe ) e and u uD on ΓD },where Pm is the finite element interpolating space.The finite-element formulation is then obtained by restricting (2.10) to V h and S h :ZhZhh w · u dΩ ΩZw s dΩ Ωwh h dΓ(2.11)ΓNIn the (Bubnow-) Galerkin method, the test functions are taken equal to the interpolation orshape functions Nj .NXuh ui Ni ,wh Nji 113

2. The Finite Volume and Finite Element MethodFigure 2.2.: Piecewise linear shape functions [11]The key is now to choose these interpolation functions in such a way that their contributionto the solution is non-zero on only a few elements, see Figure 2.2. This will then result in asparse matrix system which can be solved efficiently.It should be noted that similarly to Section 2.2, the current Section serves only to introduce thegeneral idea of the FEM, and that there are of course a variety of variations and sub-classes,notably the Petrov-Galerkin method [8].14

3. Comparison of a 2D TestcaseTo compare the OpenFOAM and COMSOL solvers, the Karman vortex street was chosen asa simple, 2D laminar test case. The Karman vortex street is a well known phenomenon andfrequently used as an example or test case. Indeed, the test case presented here was besed ona COMSOL tutorial case [12].3.1. Case SetupThe test case that was chosen consists of a two-dimensional channel with an obstacle. Fluidflows in from the left and out to the right. The channel has a length of 2.4m and a height2of 0.4m. The inflow velocity was set to 1.2 ms . The kinematic viscosity ν was set to 0.001 ms ,resulting in a Reynolds number of Re Uνl 480 (with the characteristic length l set to thechannel height). The flow is thus in the laminar regime, as the critical Reynolds number forchannel flows is in the range 2300-4000.Placing the obstacle in the middle of the channel results in a symmetric, stationary flow profile.When introducing an irregularity by moving the obstacle a little closer to the upper channelwall, however, a Karman vortex street forms.The test case was run with two variations, once with a circular obstacle and 8104 grid cells orelements, and once with a square obstacle and 32816 cells. In each instance, the grid was firstgenerated in COMSOL and then, lacking the option to convert to a format OpenFOAM canread, defined for OpenFOAM by writing the respective mesh dictionaries by hand using thesame parameters. This ensured that identical meshes were used by COMSOL and OpenFOAM.They are shown in Figure 3.1. In both cases, 6 real-time seconds of fluid flow were simulated.15

3. Comparison of a 2D TestcaseFigure 3.1.: The testcase grids3.1.1. OpenFOAMFor the OpenFOAM simulations, the icoFoam solver was selected. OpenFOAM version 2.4was used. All dictionary files needed to rerun the cases are included in appendix A.3.1.2. COMSOLCOMSOL Multiphysics version 5.0 was used. A laminar single phase flow physics interfacewas utilized, selecting incompressible flow and leaving all other settings at their defaults. Thevelocity field was initialized to 0 ms , and the inlet velocity had to be ramped up from 0 ms , asthe solver would not converge when starting with 1.2 ms right away.3.2. ResultsFigures 3.2 and 3.3 show the velocity magnitude of the fluid. It can be seen that the flowfields generally look similar, as they should. Changing the shape of the obstacle or the numberof cells does not change the result qualitatively, as the vortex street forms in both cases.However, the velocity field obtained by the simulation with OpenFOAM appears to be a bitwashed out, indicating a higher influence of diffusion.16

3. Comparison of a 2D TestcaseFigure 3.2.: The velocity magnitude after 4s. The top image corresponds to the OpenFOAMsimulation, the bottom one to COMSOL.Figure 3.3.: The velocity magnitude after 4s. The top image corresponds to the OpenFOAMsimulation, the bottom one to COMSOL.17

3. Comparison of a 2D Testcase3.2.1. Runtime ComparisonTo quantify the performance of the two solvers, the test cases were run on the RWTH cluster(on the linux frontend cluster.rz.rwth-aachen.de) using one, four and eight cores with a clockrate of 2GHz. The runtimes are shown in Table 3.2.OpenFOAM computes the solutions a lot faster than COMSOL, despite the fact that thenumber of degrees of freedom is nearly identical, see Table 3.1. This measure does of coursenot allow any statement about the solution quality.COMSOLOpenFOAMCircle (8140 Cells)2513124312Square (32816 Cells)10005698448Table 3.1.: Number of degrees of freedomFor finite volumes, the number of degrees of freedom is equal to three (the two velocity components and pressure) times the number of cells. For finite elements, the number is equal tothree times the number of nodes. Because the same (trapezoidal) grids and only linear finiteelements were used, the number of nodes is only slightly larger than the number of cells.COMSOLCircle:(8104 Cells)Square:(32816 Cells)OpenFOAMCoresShared memoryDistributed memoryDistributed �—–———————12 min1487.1 min7.1 min14833.6 min30.3 min12.1 min11.5 min2.1 min1 min4.2 min42.8 min38.5 min17.3 min4.1 min14.2 min1h 17minTable 3.2.: Runtime comparisonCOMSOL supports both shared and distributed memory parallelism, whereas OpenFOAM onlysupports the latter.18

3. Comparison of a 2D TestcaseIt can be seen that the computation runs faster when using shared memory, as expected. Furthermore, the test cases seem to be too small to gain a meaningful speedup from computationon more than four cores. In fact, OpenFOAM takes longer to compute the solution when usingeight cores, indicating that most time is utilized for the intercommunication of the cores.One other interesting thing to note is that OpenFOAM shows superlinear speedup when switching from one to four cores. This could be due to to the fact that dividing the problem leadsto fewer cache misses, avoiding reloading data from RAM.19

4. ConclusionOpenFOAM and COMSOL Multiphysics are two very different software packages. OpenFOAMuses the FVM while COMSOL uses the FEM. For general cases, it is impossible to recommend one method over the other. Historically, the FEM was predominantly used in structuralmechanics, while the FVM was preferred for fluid dynamics. Nowadays, this distinction canbe made to lesser and lesser extent, as both methods are actively researched and undergoingdevelopment, and are now more or less universally applicable. The merits of one method overthe other is the subject of ongoing debate.As for the software packages themselves, they have different strengths and foci. The maindifferences are summarized in Table 4.1.Discretization methodCostSupportGUIEase of useCustomizabilityRuntimeCOMSOLFE -OpenFOAMFV Table 4.1.: CFD-tool quick comparisonOpenFOAM emphasizes openness and the possibility to use and modify the software in anyimaginable way, at the expense of requiring a deep understanding of its inner workings.COMSOL Multiphysics, in contrast, is easy to use, but users cannot extend the functionalityand have to pay for a license.Thus, the decision which of the software packages to use, if any, can vary for every application.20

Bibliography[1] Features of OpenFOAM. http://www.openfoam.org/features/.[2] C. Greenshields.User Guide.OpenFOAM Foundation Ltd., de.pdf, 2.4.0 edition, 2015.[3] C. Greenshields. Programmer’s Guide. OpenFOAM Foundation Ltd., mersGuide.pdf, 2.4.0 edition, 2015.[4] OpenFOAM software support. http://www.openfoam.com/support/.[5] W. Schröder. Fluidmechanik, volume 7 of Aachener Beiträge zur Strömungsmechanik.Mainz, 3 edition, 2010.[6] S. M. Damian. Description and utilization of interfoam multiphase solver. Technical report, Universidad Nacional del Litoral, 477b4b948c4caa7523c8ea52.pdf.[7] B. Binninger. Numerische Simulation von Strömungsvorgängen. Vorlesungsskript, RWTHAachen, 2012.[8] Hirsch C. Numerical Computation of Internal and External Flows, volume 1 of Fundamentals of Computational Fluid Dynamics. Butterworth-Heinemann, 2 edition, 2007.[9] Finite volume methods. otes/MIT2 29F11lect 18.pdf, 2011. Lecture Notes, Massachusetts Institute of Technology.[10] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. Wiley, 2003.[11] Joaquim Peiró and Spencer Sherwin. Finite difference, finite element and finite volumemethods for partial differential equations. In Sidney Yip, editor, Handbook of MaterialsModeling, pages 2415–2446. Springer Netherlands, 2005.21

Bibliography[12] Flow past a cylinder. .mph.cylinder flow.pdf.22

A. OpenFOAM case filesA.1. Circular Obstacle/ C \ \\/ F ield OpenFOAM : The Open S o u r c e CFD T o o l b o x \\/O peration Version :2.4.0 \\ /A nd Web :www . OpenFOAM . o r g \\/M anipulation \ /FoamFile{version2.0;formatascii ;classvolScalarField ;objectp;}// //dimensions[ 0 2 2 0 0 0 0 ] ;internalFielduniform 0 lls{type}frontAndBack{type}zeroGradient ;fixedValue ;uniform 0 ;zeroGradient ;empty ;}// /// C \23

A. OpenFOAM case files \\/ F ield OpenFOAM : The Open S o u r c e CFD T o o l b o x \\/O peration Version :2.4.0 \\ /A nd Web :www . OpenFOAM . o r g \\/M anipulation \ /FoamFile{version2.0;formatascii ;classvolVectorField ;objectU;}// //dimensions[ 0 1 1 0 0 0 0 ] ;internalFielduniform ( 0 0 0 ) lls{typevalue}frontAndBack{type}fixedValue ;uniform ( 1 . 2 0 0 ) ;zeroGradient ;fixedValue ;uniform ( 0 0 0 ) ;empty ;}// ///

1. CFD Tools This chapter will give an overview over the main features of OpenFOAM and COMSOL, as well as their di erences. 1.1. OpenFOAM OpenFOAM (Open Field Operation and Manipulation) [1,2,3] is a free and open source computational uid dynamics (CFD) toolbox. It is developed by OpenCFD Ltd. and distributed