NS2DVD Navier Stokes 2D Variable Density User guide
Contents
1. t 0 completed by the specification of boundary conditions on u t 4 Finally the new density field p is computed by solving on the time interval n 1 At n 2 At the transport equation Op i divx p t u D 0 with suitable boundary conditions on p Then we go back to the first step using n 2 instead of n to compute the solution at the following time steps References 1 2 3 4 C H Bruneau and P Rasetarinera A finite volume method with efficient limiters for solving conservation laws Internal report 96024 MAB Laboratory Bordeaux 1 University France 1996 IR C Calgaro E Chane E Creus and T Goudon L1 stability of vertex based muscl finite volume schemes on unstructured grids simulation of incompressible flows with high density ratios J Comput Phys 229 17 C Calgaro E Creus and T Goudon An hybrid finite volume finite element method for variable density incompressible flows J Comput Phys 227 4671 4696 2008 Y Fraigneau J L Guermond and L Quartapelle Approximation of variable density incompress ible flows by means of finite elements and finite volumes Commun Numer Methods Eng 17 893 2001 TA J L Guermond and A Salgado A splitting method for incompressible flows with variable density based on a pressure poisson equation J Comput Phys 228 2834 2846 2009 T Schneider N Botta K J Geratz and R Klein Extension o2. 0 exp with the folowing preset parameters e zo 0 5 yo 0 the gaussian initial coordinnates e a 0 15 V2 the gaussian radius e o the gaussian amplitude NB See in file compute_sol exacte m for changing parameters Check error between two same positions of the density gaussian on the domain for comparing not the maximum error 7 1 4 Density Gaussian Rotation RRHO This benchmark is quite equivalent to the preceding one Here the density gaussian rotates counter clockwise due to the speed field 2y cos t u 2x cos t The preset domain geometry is x 0 5 12 1 0 yy 1 0 yo 0 5 and boundary conditions are Dirichlet homogeneous ones on every side for density As in the Density Gaussian Translation case due to the imposed velocity field FE part is not used Check error between two same positions of the density gaussian on the domain not the maximum error 13 Density exact solution is 2 2 x cos 2sin t y sin 2sin t zo y cos 2sin t x sin 2sin t yo 2 Pex 1 Po exp a with the folowing preset parameters e zo 0 yo 0 5 the gaussian initial coordinates 2 ea 3 0 15 ES the gaussian radius e Yo 1 the gaussian amplitude See in file compute sol exacte m for changing parameters 7 2 Other benchmarks 7 2 1 Rayleigh Taylor Instability RTIN This benchmark presented in 4 and 3 represents the evolution of two differen
3. Slope 1 H Slope 2 s Velocity A Pressure 2 10 3 Figure 7 Convergence rate 7 1 2 Lid Driven Cavity DCAV This benchmark permits to validate the finite elements scheme by observing that the intial constant density is preserved during the computation and that the results for the velocity and the pressure perfectly coincides with computations that use a standard Navier Stokes code Therefore we can know if this coupling introduces or not spurious variation of density and degradation of the accuracy On the top boundary horizontal velocity is preset to 1m s For denstity and other velocity terms boundary conditions are homogeneous Dirichlet ones 12 Isobars t 29 Density contour t 29 ll ee Figure 8 The lid driven cavity test with a homogeneous density at initial time Pressure Density Vorticity and Streamlines 7 1 3 Density Gaussian Translation TRHO This benchmark pemits to validate the finite volumes scheme It represents the translation of a density gaussian exposed to a constant horizontal speed u Boundary conditions are periodic on vertical sides and homogeneous Dirichlet on horizontal sides for density As velocity field is imposed and not computed the Finite Element part is not solved The domain geometry is preset to the square 0 25 0 25 and the horizontal speed to 1 m s Density exact solution is x zo 0 y a a2 Pex 1
4. be created in OUTPUTS directory_name ANIMATION See section 6 for more details e Save assembling matrix assembling matrix can be saved for working on linear algebra after computation In this case enter the directory name and the saving step time Matrix will be saved in OUTPUTS directory_name MATRIX Backups files are also saved for each simulation see section for more details The man ual_setup m file is created in the main workspace and a copy of this file is also created in the the directory_name eventualy created see in section 5 more details NB Those parameters are grouped in the OUTPUT structure 5 Advanced settings The GUI permits to set several basics parameters but few other ones can also be modified When user starts a simulation from the GUI some default parameters are declared in the same time see part AA of the AUXILIARY FUNCTIONS GUI computation settings m file Those parameters permit to acces to other functionalities schemes and methods Each time that a simulation is started from the GUI a summarization of all parameters is done including the hidden ones set by default This file called manual_setup m is useful for running computation faster and also for accessing to all parameters NB This file is writen thanks to the function manual_setup_writting located in AUXILIARY FUNCTIONS 5 1 Pressure constraint As seen in section 2
5. of those 3 benchmark cases in the GUI check 6 2 Animation If several iteration solutions have been saved see OUTPUTS then user can plot density islolines solutions in a figure by running plot_benchmark m User can also record a video of this simulation by running video_benchmark m 7 Benchmark cases The GUI is preset for reproducing some classical simulations some of them permits to validates schemes New benchmark cases can be implemented by following the procedure presented in subsection 7 1 Validation benchmarks 7 1 1 Analytical Benchmark EXAC This benchmark presented in 3 permits to check convergence rates by evaluating the ability of a scheme to recover analytical solution 11 The analytical solution is Pex t y pi r 0 sint _ y cost talpad x cost Dell ul sinxsinysint Where p1 r 0 2 r cos and where r 0 are the usual polar coordinates The fields pex t x y and Uex t y satisfy the mass conservation equation identically and uex t x y is solenoidal The momentum equation is satisfied with the body force defined by ysint x cos t pex t x y cos z sin y sin t x sin t y cos t pex t 2 y sin z cos y sin t fex t y If computation is performed on a disk domain preset boundary conditions are Dirichlet ones for velocity and pressure Else if domain is a square boudary values are computed thanks to the exact solution
6. step time the FV step time is tooken as Cry j el lulls dr min en Where Cry 0 1 by default and Amar is the space parameter define in remark 1 More details about the FV scheme implemented in NS2DVD are given in DI and 2 5 3 Advanced FE parameters The Finite Element part can be solved by different ways than the ones seen in section 4 4 Figure 6 shows settings posibilities The two methods which are available with the GUI are the ones in red on the scheme Other methods are accessible by modifying the parameter STOKES DIRECT INVER SION from yes to no in the FE PART of the manual_setup m file Then set the iterative methods parameters in part B 6 of the same file Distinction between different methods is done in COMPUTATION FE resol lit ns P1P2P1 m Each ending branch represents a different function some informations about algorithms used are given inside those files 10 5 4 Rerun computation In order to rerun a computation some backup files are saved in OUTPUTS BUFFER DIRECTORY FOR DAC or in OUTPUTS DIRECTORY_NAME RERUN_BACKUPS depending if a directory has been created for the simulation or not ie if we save command window informations if we save matrix or if we save solution for an animation a new directory is automatically created Those files contain variables needed for reruning a simulation from a given iteration For example solutions time t
7. t and t for numerical schemes For reruning a computation in RERUN computation settings part of the manual setup file switch RERUN variable from no to yes Then run the file a message box will open select the backup file from which one you want to restart computation The computation will stop again at the same FINAL TIME than the first computation so for a longer computation modify this parameter in PHYSICAL SETTINGS part The number of iterations between two backups can also be modify in part B 8 of the manual_setup m or in AUXILIARY_FUNCTIONS GUI computation settings m 6 Post processing Some post processing functions are implemented depending of the benchmark case Change directory to POST PROCESSING and run the corresponding function Then pick the directory name of the simulation you want to post process 6 1 Convergence rates When analytical solution is known the computed one can be compared to it Then if few simulations have been run with different edge sizes convergence rate of error depending of the mesh size can be checked e For the analytical benchmark case 7 1 1 errors between speed density and pressure analytical solution and the computed one is plotted e For the gaussian translation 7 1 3 or rotation 7 1 4 only density error depending of the mesh size is plotted because FE part is not computed For checking convergence rates of one
8. 2 an external constraint needs to be set to define pressure solution By default we set zero pressure to the nearest point of the center of the geometry But in case of a high pressure gradient in this region the zero pressure point is placed in an other part of the geometry a boundary for example Instead of imposing a zero pressure point we can also set a zero pressure average condition in order to find pressure solution This parameter is available in part B 6 of manual setup m Constrained optimization rojection methods Direct resolution Direct resolution No Yes Yes No Iterative methods stokes_direct stokes_projection_direct Bloc matrix Uzawa SDP Non SDP SDP Non SDP Conjugated Gradient GMRES Conjugated Gradient GMRES Figure 6 FE schemes tree 5 2 Advanced FV parameters 5 2 1 Flux limiter As presented in subsection 4 3 2 by default a MinMod flux limiter ensures stable simulation and avoids spurious oscillations in the vicinity of the discontinuities see 1 for more details The flowing flux limiters are also implemented e Van Leer e Van Albada e Superbee Those flux limiters can be activated in function function limiter of COMPUTATION FV 5 2 2 FV step time As seen in section 4 5 the main step time is computed as dt C hi With preset parameters C 1 and a 3 2 But the FV scheme is stable under the condition dt lt maz lullo So in order to use a suitable
9. SING This directory contains functions used for plotting data or making a film after the computation The NS2DVD directory also contains the following files e GULlauncher m Run this file for starting Graphical User Interface see section Hl e manual_setting m This file regroups all paramaters and permits to run a computation without using the GUI e main_P1P2P1 m This is the main function it can be called from the GUI or from the man ual_setup m file e user guide pdf The present documentation 4 Graphic User Interface To run GUI change directory to NS2DVD and run the GUl launcher function The first window which permits to set mesh parameters corresponds to the part A of the manual setup m file The second window which permits to set computation and outputs parameters corresponds to the part B 4 1 Mesh settings Two kinds of mesh can be built or a structured one built by hand either a non structured one built thanks to the Matlab Partial Differential Equations Toolbox PDE Toolbox An unstructured mesh can also be used without using the PDE Toolbox then mesh data files needs to be previously gen erated with a mesh editor Some mesh data files have already been generated for for unit disk and for a 0 5 0 5 x 2 2 rectangle domain Those mesh files are placed in MESH_BUILDING MESH _FILES Create a file which contains P1 nodes coordinates and a second one which contains no
10. T I VW Universit 277 V Lil let INVENTEURS DU MONDE NUM RIQUE chnologies NS2DVD Navier Stokes 2D Variable Density T 0 65 2 0 5 04 0 3 02 01 D 0 1 0 2 0 3 0 4 0 5 User guide Laboratoire Paul Painlev Contents I Introduction 2 1 General explanations 2 2 Description of the numerical scheme 2 2 1 Solving the density equation by a Finite Volume method 2 2 2 Solving the velocity equation by a Finite Element method 2 IT Settings 2 3 3 1 Two Possiblities 3 3 2 Directories and files 3 3 pA BG a GE Oe AER And ee EE ee de ad 4 4 1 1 Domain geometry 4 fh oe ee ete es Gene ee ee eee ee a eee ee 4 gee yee Ged War Pee e GE Se eh ee ee ge ee eee 5 Mak a eae a eae ben oe asm Ae ee oe ee 5 A ee teed Gece Di D at PT hse ee ee 6 43 2 PV SCHEME isis ee ee ala eel A ee ee ee oe ee ey 7 eh Mt tah he eet eed SOS EE ogee eae oe te F 7 4 4 1 Constrained optimisation 8 4 4 2 Projection method 8 4 5 Computation parameters 8 A Gee Ee A ees Se ee we ee a a 9 9 5 1 Pressure constraint 9 ne De ms Dam DS Ra DAME Ne Se y 10 DAT Flux limiter score PG a SR ee ae tn E d une ee 10 REENEN 10 SR ne Gre
11. Two Possiblities Several parameters can be set according to the simulation we want All of them are regrouped in the file called manual setup m which is preset for simulating lid driven cavity benchmark see subsection 7 1 2 It can so directly be run to get a first example of a simulation In order to simplify the setup and to exclude incompatible parameters a Graphic User Interface GUI has been developped more details about this interface are given in section 4 3 2 Directories and files Once the archive extracted move to NS2DVD directory Each m file contains an eponymous function most of them are ordered in the folowing sub directories e ASSEMBLING This directory regroups all fonctions used for assembling matrix for the FE part e AUXILIARY_FUNCTIONS Here are the functions which do not compute anything For ex ample the ones used for displaying headers or for building the Graphic User Interface e COMPUTATION Here are the computation functions grouped depending on their use Finite Elements FE Finite Volume FV etc e INITIALISATION This directory regroups functions used for setting boundary conditions and defining initial interface s equation s for bi fluid simulations e MESH_BUILDING Function used for building the mesh or loading an external mesh are grouped in this directory e OUTPUTS Files and data generated during the computation will be stored in this directory e POST_PROCES
12. des numbers of each triangle Mesh design possibilities are summarized in Figure land detailed below DOMAIN_GEOMETRY MESH_DESIGN MESH_DESIGN structured unstructured structured unstructured nbseg_c PDE Toolbox nbseg_x PDE Toolbox nbseg_y Y Ne triangles_orientation Y Ne ho hod ho ho Figure 1 Mesh parameters tree 4 1 1 Domain geometry e For a rectangle geometry set coordinates of the points taking into account that x1 um is the one on the bottom left and x2 y2 the one on the top right e For a circle geometry set the center coordinates and the radius value NB Those parameters are grouped in the domain structure 4 1 2 Mesh design e For an unstructured mesh set the maximum edge length ho e For a structured mesh hand built set the number of subsegments on the radius for a circle or on x and y axis for a rectangle geometry For a rectangle geometry choose between 2 following triangles orientations gt diagonal each quadrant is subdivised with parallels diagonals as presented in figure b gt cross each 2 x 2 rectangle cells are split by and X as presented in figure a NB Those parameters are grouped in the mesh structure Remark 1 Depending on the parameters set in the mesh setting window the edge size parameter hmaz is computed by different ways e If the domain geometry is a rectangle and the mesh is structured LQ T1 gt hmaz maz gt nbseg_x
13. e eee Ge er ar a 10 ee EE ENEE EE 11 6 Post processing 11 11 6 2 Animation cos 24 a OR A AE A EE 11 7 Benchmark cases 11 7 1 Validation benchmarkd ee 11 7 1 1 Analytical Benchmark EXAC o ae 11 7 1 2 Lid Driven Cavity DCAV 2 eked oad aa o 12 7 1 3 Density Gaussian Translation TRHO 13 7 1 4 Density Gaussian Rotation RRHO 2 22403 see eee a eae a wea 13 1 2 Other benchmarks ee 14 7 2 1 Rayleigh Taylor Instability RTIN 14 7 2 2 Falling Droplet DROP e E E ENN EEN AR HE E 15 7 3 New benchmark case 16 II Appendix 16 16 A Euler splitting 4 24 4 4 de na L a a ee eee a Ra 16 A2 n 1 velovity linear second order extrapolation 17 A 3 n 1 2 velovity linear second order extrapolati0M 17 A 4 Strang splitting oeoa aa a ee 17 Part I Introduction 1 General explanations This code collectively developped the Paul Painlev Mathematics Laboratory UMR Lille 1 University CNRS N 8524 and INRIA Lille Nord Europe is devoted to the numerical simulation of the variable density incompressible Navier Stokes system given on a domain 2 C R by Op divx pu 0 la p dru u Vx u Vxp pAxu f 1b div u 0 1c Here p t x gt 0 stands for the density of a viscous fluid whose velocity field is u t x R and pressure p t x R The descripti
14. ey are called star control volumes For a cross structured mesh see section 4 1 2 cells are octog onals and for a diagonal structured mesh cells are hexagonals as presented in figure Hai gt Or cells are built in order to have square control volumes as presented in figure 4 b a b Figure 4 Control volume on structured diagonal mesh a Star b Square Figure 5 Node identification for computing gradient 4 3 2 FV scheme e Beta is the parameter of the so called P scheme More precisely for the MUSCL strategy the gradients V py and V py are defined as V Pi BVpa 1 8 V pa V pa BVpa 1 8 Vi pa with nt SGP j l nt A j l e In case of star volume control see subsection 4 3 1 the gradient computing parameter per mits to choose how the density gradient Vp is computed Vpa Mou Mo gt MEPA M gt Either we take into account the density in the actual and in the upstream cell then we need to identify this upstream cell see Figure 5 gt Or it is computed by taking by considering the value in the actual cell and the mean value in the surounding ones e A flux limiter ensures stable simulation and avoids spurious oscillations in the vicinity of the discontinuities see 1 By default MinMod limiter is implemented but some others are disponible see section 5 2 for more details e The epsilon cri
15. f finite volume compressible flow solvers to multidimensional variable density zero Mach number flows J Comput Phys 155 248 286 1999 T5 G Strang On the construction and comparison of difference schemes SIAM J Numer Anal 5 506 517 1968 18
16. nbseg_y _ NANANAN A LIN NAN a ANINANAON NAN AN ANS a le de Le Un KARARAN LNNNN a b CH Figure 2 Structured meshes a cross b diagonal e If the domain geometry is a circle and the mesh is structured To max nbseg_c e For an unstructured mesh maximal edge size has already been declared so hmar ho 4 2 Physical parameters e Reynolds number e Final simulation time e Gravity value e Density e Maximum density in case of bi fluid flow e Boundary Conditions BC gt Velocity first composant BC gt Velocity second composant BC gt Density BC Remark 2 For most of benchmark cases boundary conditions can t be modified because possi bilities are not implemented or in some cases it does not make sense NB Those parameters are grouped in the PHYSICAL structure 4 3 Finite Volumes parameters The Finite Volume method implemented in NS2DVD is presented in 3 and 2 5 Figure 3 Control cell for unstructured mesh 4 3 1 Control cells design Control cells of finite volume scheme are vertex centered e for an unstructured mesh control cells are built by linking each barycenter to the center of its neighboring edges as presented in figure e in case of a structured mesh for both designs presented in section 4 1 2 two kinds of control cells can be built gt Either cells are built as seen previously for an unstructured mesh then th
17. o back to the first step using n 1 instead of n to compute the solution at the following time steps AA n 1 2 velovity linear second order extrapolation 1 The new density field p t is computed by solving on the time interval nAt n 1 At the transport equation p div lp 0 With u a which corresponds to the linear extrapolation of the velocity at time t 2 2 The new velocity and pressure fields u and p are computed by the resolution on the time interval nAt n 1 At of the system pr Tun 1 u 1 Mani zk Vxprtt uAxu t fort div u 0 Then we go back to the first step using n 1 instead of n to compute the solution at the following time steps A 4 Strang splitting 1 The new density field p is computed by solving on the time interval nAt n 1 At the transport equation apt divx p ttu 0 with suitable boundary conditions on p 17 2 The new velocity and pressure fields u and p are computed by the resolution on the time interval nAt n 1 At of the system p l du 1 u 1 Vu de Vip uAxu t 2 pert div u t 0 completed by the specification of boundary conditions on ui 3 The following velocity and pressure fields u and p are computed by the resolution on the time interval n 1 At n 2 At of the system pr 1 du 2 u 2 Vai dE Vxprt pA xu pota divxu
18. of the heavy fluid in a cavity The computational domain is 0 d x 0 2d where d 1 and at t 0 the fluid is at rest with density 100 f0 lt y lt 1 or 0 lt r lt 0 2 POW 1 if1i lt y lt 2 or 0 2 lt r where r y x 0 5 y 1 75 2 The volumic force in the momentum equation is pxg Boundary condition are Dirichlet ones on every boundary for speed and density On the following figure one can see the evolution of the interface plotted thanks to the plot_rtin_and_drop function e nbseg_x 6k et nbseg y 10k9 k1 k2 EN e T 0 875 15 T 0 1 T 0 5 0 5 0 5 Figure 10 Falling droplet evolution of the interface Re 3132 density ratio 100 structured cross mesh 80 x 160 density contours 1 4 lt p lt 1 6 7 3 New benchmark cases New benchmark cases can be implemented by following this brief procedure Introduce specific parameters in benchmark_initialization m if needed Identify the boundary points in set_cl_P2P1 m If needed identify the boundary points number with Dirichlet conditions on density see in set cl rho m Define the pressure constraint matrix Assemble constan Set initial values in set_init_P1P2P1 m If it is known implement the exact solution in function sol_exact m and compute_sol exact_P1P2P1 m If needed set Dirichlet boundary conditions on velocity in u_diri_P2 m Adapt the post processing Part III Appendix A Splitting time type
19. on of the external force is embodied into the right hand side f t x of 1b and y gt 0 is the viscosity of the fluid The unknowns depend on time t gt 0 and position xENCR 2 Description of the numerical scheme The mass conservation relation is solved thanks to a Finite Volume FV method and the mo mentum equation is solved with a Finite Element FE method taking into acount the divergence free constraint Lc Four different time splittings can be used for solving the system La 1b 1c FV part and FE part are solved separately successively or those part are computed alternately FV FE then FE FV Those different time splittings give different accuracy orders Time semi discreti sation is detailed in appendix A 2 1 Solving the density equation by a Finite Volume method As presented in 3 and 2 the transport equation is solved with a vertex based Finite Volume method combined with a MUSCL scheme and the use of of a flux limiter FV parameters are detailed in section 2 2 Solving the velocity equation by a Finite Element method In the numerical simulation of the Navier Stokes eq 1b a major difficulty is that the velocity and the pressure are coupled by the incompressibility constrain Lc In order to solve this system two types of schemes are implemented in the code e Constrained optimization methods e Projection methods FE parameters are detailed in part Part II Settings 3 Setup introduction 3 1
20. onverges to stationary solution This is Nn _ n l1 n_ fil n 1_ pn pur ang JET A ang LE lull Ip Ile This parameter is also used for defining if solution exploses in finite time This is described as lull gt 1 A or pl gt 1 A or pl gt 1 described as lt Remark 4 Here we consider L R norm e User can choose between the four following time splittings gt Euler splitting gt n 1 velocity linear second order extrapolation gt n 1 2 velocity linear second order extrapolation gt Strang splitting 7 NB Time discretisation of those splittings are presented in appendix A e Computation step time dt will be computed as dt C h amp ap with preset parameters maz gt C 1 which refers to the mesh regularity gt a 3 2 which permits to obtain the good ratio between the order 3 time accuracy and the order 2 space accuracy 4 6 Outputs parameters e Display for displaying figures during the computation Simulation will be slower because of the graphical treatment e Print results in a file errors and norms which are displayed in the Matlab command win dow can be printed in a text file then enter a directory name Results will be printed in OUTPUTS directory_name iterations_informations txt e Keep data for animation data can be saved for ploting an animation after the computation In this case enter the saving step time and a directory name Data files will
21. s Let us denote At the time step and t n At n gt 0 and let us assume that the numerical solution at time t namely p u p is known on the computational domain Four different splitting time type are implement in the NS2DDV code A l Euler splitting This splitting time gives a time first order precision 1 The new density field p is computed by solving on the time interval nAt n 1 At the transport equation dp div pa 0 with suitable boundary conditions on p 16 2 The new velocity and pressure fields u and p are computed by the resolution on the time interval nAt n 1 At of the system n 1 du 1 u 1 Vert a Vip Azur fort p p div u 0 completed by the specification of boundary conditions on ui Then we go back to the first step using n 1 instead of n to compute the solution at the following time steps A 2 n 1 velovity linear second order extrapolation 1 The new density field p t is computed by solving on the time interval nAt n 1 At the transport equation a div p 1u 0 With ur 2u ul which corresponds to the linear extrapolation of the velocity at time pel 2 The new velocity and pressure fields u and p are computed by the resolution on the time interval nAt n 1 At of the system n 1 du 1 u 1 Vaart n Vip _ A u t fort p u div u t 0 Then we g
22. t density fluids The heavier located on the top leaks in the light one located on the bottom The interface is slightly smoothed since we set at time t 0 polz y P Pe p PE tanh y n a 0 01d with pum gt Pm gt 0 and 7 gt 0 the amplitude of the initial perturbation The preset domain geometry is Q d 2 d 2 x 2d 2d where d 1 But for reducing compu tation time simulation can be done on Q 0 d 2 x 2d 2d Boundary condition on the top and bottom are Dirichlet ones and Neumann ones on vertical sides For reproducing this benchmark case set the following parameters e nbseg_x 4k et nbseg y 32k9 k1 k EN e Re 1000 e T 2 pmin 1 Pmax 3 On the following figure one can see the evolution of the interface plotted thanks to the plot rtin and drop function 14 T 1 T 1 5 T 2 T 2 5 T 2 75 T 3 2 2 2 2 2 2 15 7 15 q 15 q 15 q L t q 15 vr 4 ab S 4 1 au a St ab 0 5F 7 ot al 0 57 a BU 1 5 4 15 4 151 1 15 2 2 2 2 2 2 0 0 5 0 05 0 05 0 0 5 0 05 0 0 5 Figure 9 Rayleigh Taylor instability evolution of the interface Re 5000 density ratio 3 initial amplitude y 0 1 structured cross mesh 40 x 320 density contours 1 4 lt p lt 1 6 7 2 2 Falling Droplet DROP This benchmark is inspired from 6 A heavy droplet falls through a light fluid and impacts into the plane surface
23. terion defines when the gradient is too hight between 2 cells so when the flux limiter is used 4 4 Finite Elements parameters In the GUI two FE schemes are implemented in order to solve system 1b 1c 4 4 1 Constrained optimisation It can be first considered as a constrained optimization problem then a Gear scheme is implemented and the time discretization is Su 4u ul j 2At div u t 0 2 sl y a Vx ut uAxu t Van gr The matrix expression of this sytem is so eo b Leo B 0 p f p Where B represents the divergence term matrix and A the Non linear and the mass matrix The pressure term is considered as a Lagrange multiplier and the linear system is solved thanks to the Matlab direct inversion backslash operator 4 4 2 Projection method It can also be solved by using a penalization type projection method with rotational incremental scheme as presented in 5 Main steps are the followings 3t Z Ap ag 2At eg 3y 1 4u ER u P 2At 2x u uy pa 0 4 1 erte iw Art Y r Doa er et 3X A n l A y nl g JAN 0 ptt p al par uV S ut Remark 3 Some other FE schemes are implemented and also other solvers as other projection methods and iterative methods Conjugated Gradient and GMRES algorithms see section 5 4 for more details 4 5 Computation parameters e Stop criterion is used for defining if solution c
Download Pdf Manuals
Related Search