Home

LS-DYNA Incompressible Flow Solver User's Manual

Contents

1. 50 100 150 Time Figure 10 5 Kinetic energy time history 200 Incompressible Flow Solver 10 1 POISEUILLE FLOW SIMPLE DUCT ENTRANCE REGION Time 209 66 Fringe Levels Contours of X Velocity 1 500e 00 min 0 at nodes 1 S max 1 5 at node 6 1 364e 00 Vector of XY velocity min 0 at nodes 1 1 227e 00 max 1 5 at node 6 1 091e 00 9 545e 01 8 182e 01 6 818e 01 5 455e 01 4 091e 01 2 727e 01 1 364e 01 0 000e 00 Figure 10 6 X velocity field after 209 66 time units SIMPLE DUCT ENTRANCE REGION Time 209 66 in 0 0520345 at nodo 227 8 400701 max 0 594762 at nodes 1 5 436e 01 Fringe Levels 4 9246 01 4 413e 01 3 901e 01 3 390e 01 2 878e 01 2 367e 01 1 855e 01 1 343e 01 8 319e 02 3 203e 02 Figure 10 7 Pressure contours after 209 66 time units Incompressible Flow Solver 109 CHAPTER 10 EXAMPLE PROBLEMS 0 0 0 1 02 03 04 0 5 0 6 0 7 0 8 0 9 10 L1 12 13 14 L5 16 X velocity Figure 10 8 Steady state x velocity profile at outlet plane 110 Incompressible Flow Solver 10 2 NATURAL CONVECTION IN A SQUARE CAVITY 10 2 Natural Convection in a Square Cavity The thermal cavity benchmark introduced by De Vahl Davis 99 100 is used here to demonstrate the application of the explicit semi implicit and fully implicit algorithms to buoyancy driven flow Figure 10 9 shows the computational domain for the differentially he
2. The conservation of energy is expressed in terms of temperature T as oT oc uvr V 4 0 2 40 where C is the specific heat at constant pressure q is the diffusional heat flux rate and Q represents volumetric heat sources and sinks e g due to exothermic endothermic chemical reactions Fourier s law relates the heat flux rate to the temperature gradient and thermal conductivity q kVT 2 41 where k is the thermal conductivity tensor In many cases the fluid properties are only available as scalar quantities i e the thermal conductivity is considered to be isotropic 2 2 4 Boundary and Initial Conditions The prescription of boundary conditions is based on a flow domain with boundaries that are either physical or implied for the purposes of performing a simulation A simple flow domain is shown in Figure 2 1 where the boundary of the domain is I Ty Ul The momentum equations Eq 2 28 are subject to boundary conditions that consist of specified velocity on I as in Eq 2 42 or traction boundary conditions on I as in Eq 2 44 u x t x 1 on T 2 42 In the case of a no slip and no penetration boundary u 0 is the prescribed velocity boundary condition The prescribed traction boundary conditions are o n f x t on T 2 43 where n is the outward normal for the domain boundary and f are the components of the prescribed traction In terms of the pressure and strain rate the traction boundary c
3. 998 kg m and dynamic viscosity u 1073 kg m s The boundary conditions for the sheath chamber consist of no slip and no penetration conditions on the surface of the chamber The pressure at the axial inlet is prescribed to be 10 kPa while the side inlet ports have a prescribed pressure of 7 5 kPa For each inlet port the prescribed tangential velocity components are zero At the axial outlet from the chamber natural do nothing boundary conditions are used The initial conditions for the problem consist of the initial pressure difference and velocities u v w 0 The sheath p2 k keyword input file is shown below In this problem a fixed maximum CFL number is used with the time step being updated every two time steps as indicated by ickdt 2 Because a steady state solution is anticipated 10 A conjugate base vectors are used in the pressure solution strategy which yields pressure solutions in an average of 4 iterations per time step Although a steady state laminar flow solution is expected in this problem a time accurate calculation is performed to test the steady state hypothesis Figure 10 26 shows the kinetic energy time history plot for the sheath chamber and illustrates the trend towards steady conditions for t gt 0 005 s Similar results are observed for the velocity time history plots in Figure 10 27 and 10 28 For t gt 0 005 s the peak x velocity at the outlet is approximately 4 times larger than along the center line of the ax
4. In comparison to Germano s dynamic procedure which can exhibit pathalogical behavior requiring ad ditional spatial averaging in homogeneous directions the evaluation of C in Eq 7 36 is well defined and non zero 6 0 0 In addition in this dynamic procedure there is no requirement for the model constant C to be filter independent i e the evaluation is strictly local Thus the model is referred to as the local dynamic k 8 model LDKM In order to obatin the dissipation model constant Cz similarity between e and e is assumed Thus wee Ey 7 38 which permits evaulation of the dissipation constant as WEEDS e QUON qu Om 7 39 Cete p ktest 3 2 Ox Ox Ox Ox 64 Incompressible Flow Solver Chapter 8 Flow Statistics The purpose for this chapter is to provide sufficient information for an experienced analyst to use LS DYNA to generate flow statistics from a transient flow simulation The goal of this chapter is to present the built in tools that are available in LS DYNA and LS POST to analyze time dependent flow data specifically turbulent flow data It is through the use of these tools that the user will develop an under standing of the simulated flow field which may lead to design changes or implementation of turbulence controls LS DYNA provides the capability to accumulate time averaged data that may be later post processed with LS POST to derive flow statistics For large eddy simulations turbul
5. 75 BOUNDARY OUTFLOW CFD 76 BOUNDARY PRESCRIBED CFD 77 BOUNDARY PRESSURE CFD SET 79 CONTROL_CFD_AUTO 81 CONTROL_CFD_GENERAL 83 CONTROL CFD MOMENTUM 85 CONTROL CFD PRESSURE 88 CONTROL CFD TRANSPORT 90 CONTROL CFD TURBULENCE 94 INITIAL_CFD 95 MAT CFD CONSTANT 97 shared keywords 99 CONTROL OUTPUT 99 CONTROL SOLUTION 99 CONTROL TERMINATION 99 DATABASE BINARY D3DUMP 99 DATABASE BINARY D3MEAN 99 DATABASE BINARY D3PLOT 99 DATABASE BINARY D3THDT 99 DATABASE BINARY RUNRSF 99 DATABASE HISTORY NODE 99 DEFINE CURVE 99 SECTION SHELL 99 SECTION SOLID 99 LS POST statistics levels 73 higher order statistics level 3 73 mean statistics level 1 73 second moment statistics level 2 73 Material models fluid properties 97 Output ASCII 99 binary 99 d3dump 99 d3mean 65 99 d3plot 99 d3thdt 99 runrsf 99 state database 99 time history database 99 Parallelism communication costs 25 domain decomposition 22 domain decomposition message passing 21 parallel assembly via message passing 22 parallel explicit algorithm 24 Pressure solution methods 41 PPE solvers 88 pressure stabilization 42 pressure Poisson solvers 88 projection CG method 44 updating 46 saddle point problems 41 solver performance 47 Pressure Poisson equation pressure levels 40 Pressure Poisson equation PPE 13 28 31 Projection method 27 fully implicit projection method 3
6. Both the semi implicit and fully implicit methods used consistent mass matrices imass 2 for the momentum and energy equations in conjunction with the default second order time weighting parameters for the viscous diffusion and advective terms Node Set ID 2 a Th Shell Set ID 100 112 Node Set ID 4 E L Node Set ID 3 a L Tc Node Set ID 1 Figure 10 9 21 x 21 thermal cavity mesh Incompressible Flow Solver 10 2 NATURAL CONVECTION IN A SQUARE CAVITY Example box p2 k KEYWORD TITLE Thermal Cavity Ra 1000 9 5 9995599 5559955525592559555595952525955552552555555525255555555555595555599555 Setup the CFD problem soln 4 CONTROL SOLUTION soln 4 CONTROL CFD AUTO iauto epsdt dtsf dtmax 2 0 001 1 50 0 500 CONTROL CFD GENERAL insol dtinit cfl ickdt istats tstart lavg 3 1 0e 4 05 5 10 CONTROL CFD MOMENTUM imass iadvec Ifet divu thetak thetaa thetaf 2 10 1 msol maxit ichkit iwrt ihist eps ihg ehg 20 100 2 0 0 CONTROL CFD TRANSPORT itemp iden nspec ifrac 1 0 0 0 imass iadvec ifct thetak thetaa thetaf 2 10 1 itsol maxit ichkit iwrt ihist eps ihg ehg 100 2 0 0 CONTROL CFD PRESSURE ipsol maxit ickit iwrt ihist eps 11 nvec stab beta sid plev lcid 100 0 0 2 CONTROL TERMINATION endtim endcyc dtmin endeng
7. and the base vectors in the A norm Here is a set of A conjugate vectors derived from N prior solutions to Ax b where i 1 N As will be demonstrated in 6 4 N 5 to N 10 is a reasonable choice Let the initial guess for a solution at time step n be given by N xe S 00 6 7 i l Let the difference between the solution at time step n and the initial guess x be defined as lella x lla 6 8 where lla V C AC With this definition pe AA any Ay x Ax x Ax x Ax 6 9 and for A being symmetric positive definite Ix Al 4 Ag AMY Ax 4 x Ax 6 10 Incompressible Flow Solver 45 CHAPTER 6 PRESSURE SOLUTION METHODS Now substituting Eq 6 7 and using the fact that b is a set of A conjugate vectors N N lela Ge Ax 2 V ao Ax M 07 6 11 i l i Minimization of e 4 with respect to 01 yields o 6 Ax 9 b 6 12 Thus given a set of A conjugate vectors O the best approximation to x that minimizes the error in the A norm is obtained by projecting the right hand side b at time step n onto the set of base vectors b This suggests the following solution procedure Algorithm 1 A conjugate projection CG method 1 oj 0j b 2 xem ye 00 where N is the number of A conjugate base vectors 3 Solve AAx r using the conjugate gradient method where r b Ax 4 Update the solution x x Ax
8. and the referential coordinates 5 7 in Eq 3 23 Incompressible Flow Solver 17 CHAPTER 3 EXPLICIT TIME INTEGRATION 1 EX EL ye E ge JO z nix ny nz 3 23 Cixe CT ye C zt The element volume is simply the determinant of the Jacobian in Eq 3 24 Note that the element volume in three dimensions may only be computed exactly with single point integration for elements that are bricks or parallelepipeds V det J 0 3 24 The computation of the element level gradient operators proceeds by first evaluating the co factors of the Jacobian the inverse Jacobian in Eq 3 25 and the gradient operators as shown in Eq 3 26 Again only the unique gradient operators must be stored 1 e 12 words of storage per element are required for Cy Cy and C7 1 D Dj J 0 3 25 1 C g Dust Din Dist 1 G 3 D216 Don D234 3 26 1 QE z P316 Daan Das The computation of the element level mass matrix using one point quadrature and row sum lumping yields the operator for 2 D in Eq 3 27 and the operator for 3 D in Eq 3 28 Af ab Saby 3 27 Mg Ba 3 28 The direct evaluation of the advection operator in Eq 3 18 requires an integral of triple products that is very computationally intensive Therefore the advection operator is approximated using an ad hoc modification known as the centroid advection velocity This modification assumes that u in Eq 3 18 may be approximated by Nnpe
9. as in Eq 2 15 or traction boundary conditions on I as in Eq 2 17 uj xi t i x t on Ty 2 15 In the case of a no slip and no penetration boundary u 0 is the prescribed velocity boundary condition The prescribed traction boundary conditions are ojjnj fixit on I5 2 16 Incompressible Flow Solver 5 CHAPTER 2 GOVERNING EQUATIONS r 2 Outflow gt Boundary i I Figure 2 1 Flow domain for conservation equations where n is the outward normal for the domain boundary and fi are the components of the prescribed traction In terms of the pressure and strain rate the traction boundary conditions are p j 2uSij nj fi xi t on Pa 2 17 The traction and velocity boundary conditions can be mixed In a two dimensional sense mixed boundary conditions can consist of a prescribed normal traction and a tangential velocity For example at the outflow boundary in Figure 2 1 a homogeneous normal traction and vertical velocity on I constitutes a valid set of mixed boundary conditions A detailed discussion of boundary conditions for the incompressible Navier Stokes equations may be found in Gresho and Sani 5 Turning attention to the species transport equations boundary conditions for Eq 2 9 may consist of either a prescribed concentration or a mass flux rate In the binary mixture example the prescribed con centration is Zi xit Zi xit 2 18 where Z is the known va
10. boundaries on the channel and cylinder walls Natural boundary conditions were used at the channel outflow The prescription of ujn 0 on the channel walls results in a two dimensional checkerboard pressure mode that exists in each plane of three dimensional elements normal to the flow direction In this situation pressure stabilization is necessary to filter the checkerboard modes For this reason all computations were performed using the pressure stabilization outlined in 6 2 Tables 6 3 and 6 4 report the grind times and iteration counts for the channel flow problem using the explicit and semi implicit projection methods respectively The results in Table 6 3 show that the cost of the PPE solution dominates the solution time per time step In contrast the solution of the stabilized PPE in the projection algorithm is less than one third of the computational cost per time step when the ESSOR PCG solver is used Again the best overall grind times are found using the ESSOR PCG solver combined with 5 to 10 base vectors where a speedup of 3 to 4 is observed relative to the basic JPCG solver The influence of the large half bandwidth of the PPE is seen on the grind times resulting in normalized grind times approaching unity for the ESSOR PCG solver Based on the results for a range of flow problems the use of a stabilized consistent pressure Poisson operator an A conjugate projection technique and SSOR preconditioning with the conjugate gradient m
11. t u 4 38 2 Solve the momentum equations treating the simplified form of the linearized equations and perform the projection step 1 e M At OKK 4 94A u M At 1 0x K D RS u u Ar OrF 1 0g F 4 MM CP 4 39 where une ees 4 40 2 3 Solve for the Lagrange multiplier CM Che Cr tao 4 41 4 Perform the velocity update ug MDICA 4 42 After the velocity update an updated pressure at time level n 4 1 1s obtained via A pedet TP 4 43 5 Update the acceleration as 2 u y pies c em 4 44 Arn 2 ey 6 Update the time step based on the user specified local time truncation error according to 3 Apt At l 4 45 Lem P where Hyp 1 Ndof 1 umm suy d u di ur 4 46 Y M Nnp ki max max ut AD Ndof is the number of degrees of freedom DOF per node e g u v w T Zl and Nnp is the number of unknowns per DOF Here d is the error between the corrected and predicted values of the field variables utt Qrt aure E AM M 4 47 3 1 LAM In LS DYNA the time step prediction is limited by a maximum time step scale factor 1 l P DTSF 4 48 LL lt dp where 1 5 lt DTSF lt 2 0 In addition a user specified upper limit Atmax is placed on the time step size 34 Incompressible Flow Solver Chapter 5 Boundary Conditions and Source Terms This chapter summarizes the boundary cond
12. 01 shell 434 min 7 53e 01 shell 806 Surface middle 695e 01 Surface middle phi 1 7 04e 01 Coorse mesh vortex shedding a t 2003006402 b mox 3 51e 01 shell 571 phi 2 mox 1 70e 01 shell 1441 phi 3 min 2 48e 01 shell 670 min 5 71e 01 shell 631 Surfoce middle 351e 01 Surfoce middle 170 01 300e 01 200e 01 149e 08 c d mox 2 28e 01 shell 1441 phi 4 mox 4 88e 01 shell 590 min 4 72e 01 shell 610 min 4 52e 01 shell 589 Surface middle 228e 01 Surface middle phi 5 1 00e 01 e f Figure 6 4 a Snapshot of the pressure field during vortex shedding and five fields b f based on pressure solution at the five prior time steps 48 Incompressible Flow Solver 6 4 PPE SOLVER PERFORMANCE a Temperature b Vorticity c Pressure Figure 6 5 Snapshot of a temperature field b z vorticity field and c pressure for the 2 D momentum driven jet Re 3561 Fr 326 The temperature vorticity and pressure fields have been reflected about the vertical centerline for presentation purposes The first series of computations used to evaluate the PPE solution methods were carried out for a momen tum driven slot jet with Re 3561 based on a 15 mm slot width and Fr 326 Fr v gBAT Le where g 9 81 m s BAT 1 3 Le 15 mm and v 4 m s In this computation an energy equation was solved in conjunction with the Navier Stokes equations using
13. 10 um 2 12 2 1 3 Energy Conservation Conservation of energy is expressed in terms of temperature T as oT oT Ogi ocr Sr euer SH 49 2 13 where Cp is the specific heat at constant pressure q is the diffusional heat flux rate and Q represents volumetric heat sources and sinks e g due to exothermic endothermic chemical reactions Fourier s law relates the heat flux rate to the temperature gradient and thermal conductivity oT UD 2 14 qi k where k is the thermal conductivity tensor In many cases the fluid properties are only available as scalar quantities i e the thermal conductivity is considered to be isotropic Remark In LS DYNA the viscosity thermal conductivity and mass diffusivities are treated as second rank tensors even though these properties may only be available as scalar quantities for some fluids In the limiting case of a scalar material property such as viscosity the internal code representation assumes that the viscosity is u j 6 f where fiis the user input scalar viscosity 2 1 4 Boundary and Initial Conditions The prescription of boundary conditions is based on a flow domain with boundaries that are either physical or implied for the purposes of performing a simulation A simple flow domain is shown in Figure 2 1 where the boundary of the domain is I Ty UI The momentum equations Eq 2 1 are subject to boundary conditions that consist of specified velocity on I
14. 6 1a Here I7 represents the shared inter element boundary is the jump operator B is a non dimensional parameter and wy is the pressure basis function for element In two dimensions Izy represents the length of the 42 Incompressible Flow Solver 6 2 0100 STABILIZATION element edge shared by element and J and in three dimensions it represents the area of the shared face The inclusion of the PPE term C M C in Eq 6 5 yields proper dimensionality of the stabilization matrix accounts for scaling due to irregular elements and still preserves the symmetry of the original PPE For alternative scaling procedures see Chan and Sugiyama 62 The diagonal contributions are calculated as CTM C Su py EM Chul i wil volar 6 6 J IJ Ty where the sum on J indicates a summation over all shared faces of element This formulation preserves the row sum zero property of the original PPE The global jump stabilization is simple to implement from a linear algebra point of view since the contributions to S are simply augmented row entries from the original PPE Pressure DOF Velocity DOF Figure 6 1 Element configuration for pressure stabilization a global jump b local jump In contrast to the global method the local stabilization procedure relies on the construction of macro elements that contain one velocity node per edge of the macro element in two dimensions and one velocity node per face in thr
15. 9732E 01 div u 8 u_i mn mx 367E 01 0 CHAPTER 10 EXAMPLE PROBLEMS 0000E 00 pmn pmx 115E 02 0 115E 02 000E 00 0 000E 00 0 000E 00 100E 01 dt 1 00E 04 write d3plot file 4427E 11 pmn pmx 126E 02 0 126E 02 434E 02 927E 02 0 927E 02 100E 01 1 5000 old dt 1 0000E 04 new dt 1 5000E 04 at node 311 Y Velocity at node 149 Y Velocity 4505E 11 pmn pmx 147E 02 0 147E 02 907E 02 192E 01 0 192E 01 100E 01 1 5000 old dt at node at node 1 5000E 04 new dt 15 Y Velocity 428 Y Velocity 2 2500E 04 5514E 11 pmn pmx 168E 02 0 168E 02 166E 01 347E 01 0 347E 01 100E 01 1 4996 old dt 2 2500E 04 new dt at node 407 Y Velocity at node 38 Y Velocity 3 3742E 04 3505E 10 pmn pmx 103E 03 0 875E 02 367E 01 367E 01 0 367E 01 tmn mx 0 000E 00 0 100E 01 LTE 1 0034E 03 DTSF Min diff 6 7904E 03 Max diff 6 7904E 03 t 1 9876E 01 div u 8 u_i mn mx 367E 01 0 tmn mx 0 000E 00 0 LTE 1 0034E 03 DTSF Min diff 6 7883E 03 Max diff 6 7883E 03 0 9989 old dt 1 4461E 01 new dt at node 139 Y Velocity at node 322 Y Velocity 1 4445E 01 3491E 10 pmn pmx 103E 03 0 875E 02 367E 01 367E 01 0 367E 01 100E 01 0 9989 old dt 1 4445E 01 new dt at node 322 Y Velocity at node 139 Y Velocity 1 4429E 01 Figure 10 10 Screen output of RMS divergence minimum maximum pressure velocity temperature and time step control monitors using the Adams Bashforth
16. Dynamics Review 97 chapter Problems and solutions generalized and FEM related to rapid and impulsive changes for incompressible flow John Wiley and Sons New York New York 1997 11 Mark A Christon and Thomas E Spelce Visualization of high resolution three dimensional nonlin ear finite element analyses In Visualization 92 pages 324 331 Boston MA 1992 IEEE LLNL UCRL JC 110110 141 BIBLIOGRAPHY 12 Mark A Christon Visualization methods for high resolution transient 3 D finite element simu lations In International Workshop on Visualization Paderborn Germany January 1994 invited paper LLNL UCRL JC 119879 13 Gerald L Goudreau and John O Halquist Recent developments in large scale finite element lagrangian hydrocode technology Computer Methods in Applied Mechanics and Engineering 33 725 1982 14 T Belytschko J S Ong W K Liu and J M Kennedy Hourglass control in linear and nonlinear problems Computer Methods in Applied Mechanics and Engineering 43 251 276 1984 15 W K Liu and T Belytschko Efficient linear and nonlinear heat conduction with a quadrilateral element International Journal for Numerical Methods in Engineering 20 931 948 1984 16 W K Liu J S Ong and R A Uras Finite element stabilization matrices Computer Methods in Applied Mechanics and Engineering 53 13 46 1985 17 Robert G Whirley and Bruce E Engelmann Dyna3d A nonlinear explicit t
17. GENERAL 36 ee E ER EMER ES 83 9 1 6 CONTROL_CFDMOMENTUM 0040 85 9 1 7 CONTROL CFD _PRESSURE 6 26 os 6 oe Ve 9 EA 88 9 1 8 CONTROL_CFD_TRANSPORT llle 90 9 1 9 CONTROL_CFD_TURBULENCE 4 94 91 10 INITIAL ED Saito oe xt aote PAM PER A AAA AAA 95 J LIT FMATECED OPTION sito Spic tegi At eh A Pa RE EIER RES 97 2 Shared Keywords s pi op hth tats p ade pter ees oer pie stage Seat dendo de estes 99 10 Example Problems 101 TOST Poiseuille FIOW dcm ge qae IE eI y EET TEE REEL STRE dC Y qe Re oer 101 10 2 Natural Convection in a Square Cavity 2 2 111 10 3 The Momentum Driven Jet zzz 54 6 RE ER UE 121 10 4 The Sheath Flow Chamber ooo emo eye an em eee RR GR s 131 A Sense Switch Controls 139 vi Preface The incompressible flow solver in LS DYNA originated from research conducted by the author at Lawrence Livermore National Laboratory and Sandia National Laboratories during the 1990s This re search was done in collaboration with Phil Gresho and Stevens Chan at Lawrence Livermore National Laboratory The author is indebted to Phil Gresho for his support and his willingness to share his insights into the behavior of time dependent incompressible flows Many of the algorithmic aspects of the explicit flow solution algorithm derive directly from John Hallquist s pioneering work with low order finite ele ments and reduced integration in DYNA3D I would like to ack
18. Hp pur AN Node Set ID 20 Node Set ID 10 Figure 10 17 Momentum driven jet mesh with 11250 elements and 11466 nodes Incompressible Flow Solver 123 CHAPTER 10 EXAMPLE PROBLEMS o Kinetic Energy Temperature MOMENTUM DRIVEN JET 0 12 J 0 1 0 08 G D 2 0 06 0 04 0 02 0 1 0 0 2 0 4 0 6 Time Figure 10 18 Time history plot of kinetic energy MOMENTUM DRIVEN JET A d Node No A 11 27 380 51 i l 340 J 320 J 300 i 0 0 005 0 01 0 015 0 02 Time Figure 10 19 Nodal time history plot of temperature along the jet axis for O lt t lt 0 02 s at node 11 x y 7 50e 4 6 3 le 3 node 27 x y 7 50e 4 1 07e 2 node 51 x y 7 50e 4 1 52e 2 and node 83 x y 7 50e 4 1 98e 2 126 Incompressible Flow Solver 10 3 THE MOMENTUM DRIVEN JET MOMENTUM DRIVEN JET X velocity eo 1 2 osi uedwog i Time EO 3 MOMENTUM DRIVEN JET 2 5 Y velocity 2 9 0 0 2 0 4 0 6 b Figure 10 20 Nodal time history plot of a x velocity and b y velocit
19. IADVEC 40 2 For IAUTO 3 the default maximum time step DTMAX is set 10 times larger than the starting time step DTINIT 82 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 5 CONTROL CFD GENERAL Purpose Set the solver parameters for the Navier Stokes flow solver CONTROL CFD MOMENTUM CONTROL CFD TRANSPORT and CONTROL_CFD_PRESSURE are used in conjunction with this keyword to control the flow solver options Material models may be specified with the MAT CFD OPTION keyword input and turbulence models are activated with the CON TROL_CFD_TURBULENCE keyword input Card Format Mitre rata ESEA epee oe fama II Ss Sa ie eyed Ea ae VARIABLE DESCRIPTION INSOL Set the solver type EQ 0 INSOL 3 default EQ 1 Explicit transient incompressible Navier Stokes EQ 3 Semi implicit fully implicit transient incompressible Navier Stokes using staggered velocity pressure DTINIT Set the initial time step for the Navier Stokes and all auxiliary transport equations The time step is computed based on either the prescribed CFL number INSOL 3 or stability INSOL 1 unless ICKDT lt 0 or IAUTO 3 on the CONTROL_CFD_AUTO keyword CFL Set the maximum grid CFL number to be maintained during the computation EQ 0 CFL 0 9 default for INSOL 1 CFL 2 default for INSOL 3 Incompressible Flow Solver 83 CHAPTER 9 KEYWORD INPUT VARIABLE ICKDT IACURC Remarks 1 There are multiple solver optio
20. R Leland A multilevel algorithm for partitioning graphs Technical Report SAND93 1301 Sandia National Laboratories Report October 1993 26 B Hendrickson and R Leland The chaco user s guide version 1 0 Technical Report SAND93 2339 Sandia National Laboratories Report October 1993 27 Livermore Software Technology Corporation LS DYNA Keyword User s Manual Nonlinear Dy namic Analysis of Structures Livermore Software Technology Corporation 7374 Las Positas Road Livermore California 94550 USA version 960 edition 2001 28 Alexandre Joel Chorin Numerical solution of the navier stokes equations Mathematics of Com putations 22 745 762 1968 142 Incompressible Flow Solver BIBLIOGRAPHY 29 J Van Kan A second order accurate pressure correction scheme for viscous incompressible flow SIAM Journal for Scientific and Statistical Computing 7 870 891 1986 30 John B Bell Philip Colella and Harland M Glaz A second order projection method for the incompressible navier stokes equations Journal of Computational Physics 85 257 283 1989 31 Philip M Gresho Stevens T Chan Mark A Christon and Allen C Hindmarsh A little more on stabilized q1q1 for transient viscous incompressible flow International Journal for Numerical Methods in Fluids 21 837 856 1995 32 Philip M Gresho and Stevens T Chan Projection 2 goes turbulent and fully implicit preprint International Journal for Computa
21. The mass conservation principle in divergence form is Piy pVu 0 2 33 In the incompressible limit the velocity field is solenoidal V u 0 2 34 which implies a mass density transport equation op tu Vp 0 2 35 For constant density Eq 2 35 is neglected with Eq 2 34 remaining as a constraint on the velocity field LS DYNA provides the capability to transport up to 10 species represented by the mass concentration Z1 Z2 Z10 In order to simplify the presentation a single mass fraction is presented representing a binary mixture In order to account for the change in mass concentration mass conservation applied to the individual species yields for Zi OZ pS Pu VZ VJi r 2 36 where J is the diffusional mass flux rate and r1 is a volumetric mass source The mass source may include the injection of mass concentration from a boundary or the source sink terms from chemical reac tions The diffusional mass flux rate is based on Fick s law of diffusion Ji pD VZ 2 37 8 Incompressible Flow Solver 2 2 CONSERVATION EQUATIONS VECTOR NOTATION where D is a tensorial mass diffusivity Typically mass diffusivities are only available as scalars so that Ji p 2 VZ 2 38 In the most general form the species concentration transport equations are oZ P pu VZ 2 V Jp rn 2 39 where indicates the species mass concentration i e Z 1 2 10 2 2 3 Energy Conservation
22. Trapezoid rule time integrator Incompressible Flow Solver 10 2 NATURAL CONVECTION IN A SQUARE CAVITY Algorithm Atmin Amar Umar Vmax De Vahl Davis 100 40 x 40 grid Explicit Algorithm 2 00e 4 2 00e 4 3 619 Semi Implicit Algorithm 1 00e 4 6 92e 3 3 629 Fully Implicit Algorithm Table 10 1 Maximum velocities for the 21 x 21 grid using a the explicit method b semi implicit projec tion and c the fully implicit projection method with the AB2 TR time integrator a THERMAL CAVITY RA 1000 2 5 2 gt o o 15 c ul 2 9 Xx 1 0 57 0 1 9 0 5 10 15 20 El 3 Time 0002 82 11 60 pziz Figure 10 11 Thermal cavity kinetic energy time history Incompressible Flow Solver 117 CHAPTER 10 EXAMPLE PROBLEMS THERMAL CAVITY RA 1000 Time 19 152 Fringe Levels Contours of Stream Function min 1 16467 at node 121 oe max 1 45795e 14 at node 332 1294e 01 2 588e 01 3 882e 01 5 176e 01 6 470e 01 7 764e 01 9 059e 01 1 035e 00 1 165e 00 Figure 10 12 Thermal cavity stream function contours THERMAL CAVITY RA 1000 Time 19 152 Fringe Levels Conto f Te tul imn 006 0 max 1 at nodes 1 9 000e 01 8 000e 01 7 000e 01 6 000e 01 5 000e 01 4 000e 01 3 000e 01 2 000e 01 1 000e 01 0 000e 00 Figure 10 13 Thermal cavity temperature contours 118 Incompressible Flow Sol
23. Velocity Average 0 50 100 150 200 250 300 Time Figure 8 1 Time averaged and instantaneous velocity where Q is the volume of the domain Of course homogeneous flow is an idealization because boundary conditions such as no slip no penetration surfaces introduce flow inhomogeneities For flows that are both stationary and homogeneous all three forms of averaging are equivalent i e the ergodic hypothesis applies The current version of LS POST calculates Reynolds averaged statistics of stationary flows Spatial averaging of homogeneous data ensemble averaging of transient flow data and calculation of probability densities have been neglected in the current version of the software 8 1 2 Mean and Fluctuating Quantities Turning now to the calculation of mean and fluctuating components the definition of the Reynolds average in Eq 8 1 permits any variable can be decomposed into its mean and fluctuating component 0 x 1 lt 6 x gt 4 xt 8 5 where lt gt is the mean and 9 is the fluctuating part Given the definition of the Reynolds average application of the average to the fluctuating component yields lt q gt 0 8 6 Repeated application of the averaging procedure does not change the average lt lt gt gt lt 0 gt 8 7 Averaging a product of variables i e and y yields lt 0y gt lt 0 gt lt w gt lt 0w gt 8
24. and suggest several remedies In LS DYNA outflow boundary conditions that similar to the pressure boundary condition in Eq 5 4 are used However the outflow boundary conditions use a known pressure distribution that has been computed during the solution procedure This boundary condition results in an equilibrating force applied on the outflow boundary F Na p nj 5 5 D where p is the pressure field computed at time t The outflow boundary condition is activated using the BOUNDARY OUTFLOW CFD SET keyword which applies the outflow boundary conditions on a segment set basis 5 5 Body Forces The presence of body forces in the momentum equations results in element level force contributions At the element level the forces are Fe Na pfi dQ 5 6 Q where f represents the body force per unit volume In thermal convection problems where a Boussinesq fluid is appropriate the body force due to buoyancy is Fe Jj Na pgi B T Tref dQ 5 7 where g is the acceleration due to gravity D is the coefficient of thermal expansion and T y is a reference temperature The specification of the gravity vector reference temperature and fluid properties is achieved with the MAT_CFD_OPTION keyword 38 Incompressible Flow Solver 5 6 FLUX BOUNDARY CONDITIONS Remark By default the fluid temperature is initialized to the reference temperature T ef specified in the MAT_CFD_OPTION keyword Each fluid is i
25. constant support for pressure The methods for obtaining the weak form of the conservation equations are well known and will not be repeated here see for exam ple Gresho et al 7 Hughes 8 and Zienkiewicz and Taylor 9 The spatially discrete form of Eq 2 1 and 2 7 are Mu A uju Ku Cp F 3 1 and Clu g 3 2 where M is the mass matrix A u and K are the the advection and the viscous diffusion operators respec tively F is the body force and g accounts for the presence of prescribed velocity boundary conditions C is the gradient operator and C is the divergence operator Here u and p are understood to be dis crete approximations to the continuous velocity and pressure fields Equations 3 1 and 3 2 constitute a differential algebraic system of equations that precludes the direct application of time marching algo rithms due to the presence of the discrete incompressibility constraint Following Gresho et al 1 a consistent discrete pressure Poisson equation PPE is constructed using a row sum lumped mass matrix Mz ic M c p CM F Ku A u ju 3 3 Here g accounts for time dependent velocity boundary conditions The PPE constitutes an algebraic system of equations that is solved for the element centered pressure during the time marching procedure Figure 3 1 shows the dual staggered grid associated with the pressure variables The PPE in Eq 3 3 incorporates the effect of the essential velocity bound
26. create each d3mean time average data window Hence this derived data is also a subset of the input data from the LS DYNA flow simulation The level 1 derived mean quantities include the magnitude of the vorticity lt gt enstrophy 1 2 lt 0 gt lt gt and helicity lt gt lt uj gt At the intermediate level of derived statistics level 2 second moments are computed in addition to the mean quantities just mentioned In this case the turbulent stress tensor and scalar flux vectors are com puted That is lt un gt lt uip gt lt uT gt lt uz gt and the turbulent kinetic energy are added to the output time average database At the highest level of statistics skewness and flatness of the velocity are computed in addition to the mean and derived flux quantities Currently the scalar dissipation is not treated in LS POST Thus the time averaged statistics derived by LS POST from the d3mean database are level 1 the magnitude of the vorticity lt gt enstrophy 1 2 lt 6j gt lt 0 gt and helicity lt gt lt u gt Table 8 4 shows these mean variables and the associated LS POST labels level 2 In addition to the mean quantities for level 1 the turbulent stress tensor turbulent kinetic energy and scalar flux vectors are computed That is lt un Sod lt wp gt lt uT gt lt uz gt are added to the time average database Table 8 5 shows th
27. element mesh during the initialization phase This requires the nodal as sembly of the partial nodal mass at nodes on sub domain boundaries as shown in Figure 3 4 3 20 4 Communication Costs This section outlines the communications costs associated with the explicit time integration algorithm The communications costs may be broken down into the cost per time step for the momentum equations and the cost per time step for solving the PPE To begin Nr is the number of nodes on the boundary of a sub domain In two dimensions the average number of nodes communicated to an adjacent processor is Nr 4 and in three dimensions the average is Nr 6 Assuming that there are 8 bytes per floating point word then the total number of bytes per processor to be communicated via a send is N 8NrNpor 3 46 where Npor is the number of degrees of freedom per node The communication cost for a send operation may be broken into three parts the time to initiate the message passing fsrartup the cost per packet tpacket and the transmission time frransmir Thus the time to send a message 1s Ne t packet N cltransmit 5 3 AT tsend startup N packet where N packet is the number of bytes per packet From Eq 3 46 and 3 47 it is clear that the number of adjacent sub domains determines the total startup time for message passing Next the communication cost per time step is estimated for the explicit algorithm For the momentum equat
28. eta grid CFL number Element TOs seared aie avt re RU ADU idad Element Lengths xa vns RA eye eA nn Grid GEL Number weed usa e tos A re ede xi 101 8 7298E 02 9 2771E 01 2 0000E 00 1 1877E 01 200 5 3391E 01 1 0060E 01 97 1 0000E 01 3 5299E 05 101 8 7298E 02 2 0000E 00 97 1 0000E 01 4 1926E 06 Figure 10 2 Screen report of grid parameters showing the minimum and maximum grid Reynolds and CFL numbers and and associated element numbers 106 Incompressible Flow Solver 10 1 POISEUILLE FLOW initialization completed 0 t 0 0000E 00 div u 0 0000E 00 pmn pmx 0 320E 01 0 595E 00 u_i mn mx 0 000E 00 0 000E 00 0 000E 00 0 000E 00 0000E 00 dt 1 00E 02 write d3plot file 0000E 02 div u 8 4941E 21 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 120E 02 481E 09 0 426E 09 0000E 02 div u 1 6672E 20 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 240E 02 954E 09 0 848E 09 0000E 02 div u 2 9852E 20 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 360E 02 142E 08 0 127E 08 0000E 02 div u 4 5121E 20 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 480E 02 188E 08 0 168E 08 0000E 02 div u 6 5703E 20 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 600E 02 234E 08 0 209E 08 ho N A WwW cr amp tat ate cr ct u H U H d Be WP BH oO T 1995 2 0907E 02 div u 3 7079E 16 pmn pmx 0 320E 01 0 595E 00 mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 2 0919E 02 div u 3 7104E 16 pmn
29. gt uaT gt lt uZ gt lt UZ gt lt uaZ gt lt uZ gt lt u gt lt uz gt uZ1o gt lt u Z40 gt lt u3Z10 gt lt ulul gt lt ulu2 gt lt ulu3 gt u2u2 u2u3 u3u3 ulp lt u2p gt lt u3p gt lt ulT gt lt u2T gt lt u3T gt lt ulZ1 gt lt u2Z1 gt lt u3Z1 gt lt ulZ2 gt u222 lt u3Z2 gt lt ulZ10 gt lt u2Z10 gt lt u3Z10 gt Table 8 2 LS DYNA ISTATS 2 output variables and LS POST labels 8 2 3 Higher Order Statistics The computation of the higher order statistics includes the variables listed in Table 8 3 This data along with the variables listed in Table 8 2 must be present in the input time average database in order for the statistics associated with ISTATS 3 to be computed Variable Name LS POST Label 12 Third Moments Fourth Moments Table 8 3 LS DYNA ISTATS 3 output variables and LS POST labels Incompressible Flow Solver 8 3 LS POST STATISTICS LEVELS Variable Name LS POST Label Mag vorticity 0 vorticity Enstrophy 1 2 lt 0 gt lt 0 gt lt enstrophy gt Helicity lt u gt lt gt lt helicity gt Table 8 4 LS POST level 1 derived variables and LS POST labels 8 3 LS POST Statistics Levels The derived statistics generated by LS POST from the d3mean database are also provided in three lev els The derived statistics are based on the same user specified time average window used to
30. information and in this sense the procedure may be viewed as an approximate means of deflating the eigenvalue spectrum for the PPE The combination of the A conjugate projection method PPE stabilization and SSOR preconditioning has proved to be the most computationally efficient method for solving the PPE 6 4 PPE Solver Performance This section summarizes the performance of the PPE solvers via several computational experiments per formed using the A conjugate projection strategy with several variants of the preconditioned conjugate gradient method In addition the effect of pressure stabilization on the iterative solver is considered for both local and global stabilization The CG variants include a diagonally scaled conjugate gradient JPCG solver an SSOR preconditioned CG method and an SSOR preconditioned CG solver implemented using the Eisenstat transformation the ESSOR CG solver Each CG solver used a row compressed storage scheme for the PPE The SSOR preconditioners were implemented in a stand alone form SSOR PCG that requires separate matrix vector preconditioning functions and in the ESSOR PCG form using the transformation first suggested by Eisenstat 68 to reduce the operations count in the preconditioned con jugate gradient method See also Ortega 69 and Eisenstat et al 70 Incompressible Flow Solver 47 CHAPTER 6 PRESSURE SOLUTION METHODS mox 6 95e 01 shell 120 Pressure mox 7 04e 01 shell 100 min 6 87e
31. lead to singularity in the assembled global operator In two dimensions there is only one improper singular mode while in three dimensions there are four singular modes These modes are commonly referred to as hourglass modes and when excited in a numerical solution they can remain undamped and pollute the field solution A detailed discussion of the hourglass stabilization methods used in LS DYNA is beyond the scope of this chapter However for the sake of completeness a brief overview of the so called h stabilization is presented The term h stabilization derives from the fact that the outer product of the element hourglass vectors is used to form the stabilization operator In two dimensions the single hourglass mode is Incompressible Flow Solver 19 CHAPTER 3 EXPLICIT TIME INTEGRATION r 1 1 1 1 3 35 This mode when excited may be detected visually as a w mode in isometric plots of the field variable from a numerical solution Following Goudreau and Hallquist 13 and Gresho et al 1 the element level stabilization operator is formed using the outer product of the hourglass vector H Engt ATIT 3 36 a Eng is a non dimensional parameter that in practice is unity by default for outer product stabilization see Gresho et al 1 In three dimensions the four hourglass vectors are re bed 1 1 BS hahaa V 1 T dT bb ci d T 3 37 r4 1 1 1 1 1 1 1 1 The resulting 3 D hourg
32. predicted velocity at t At using the known initial velocity and acceleration u Aou 4 32 2 Solve the momentum equations for M At OKK 94A u M Ar 1 0 amp K 1 04 A u u Ar OF 1 6 F MM Cp 4 33 where 0 lt 0 amp 07 04 lt 1 with 04 0x 0r 1 2 for second order accuracy and l 0 w E 4 34 Remark Although the fully discrete momentum equations are presented using a generalized trapezoidal rule in Eq 4 33 and below in Eq 4 39 the second order Adams Bashforth predictor re quires that 04 Oy Or 1 2 be used This is the default for the fully implicit projec tion method For additional information see the CONTROL CFD MOMENTUM and CON TROL_CFD_TRANSPORT keywords in Chapter 9 3 Solve Eq 4 19 for the Lagrange multiplier and project the approximate velocity to obtain u M C 4 35 After the velocity update an updated pressure at time level n 4 1 1s obtained via MH pipe 4 36 p power 4 36 4 Update the acceleration as wv il 4 37 4 2 3 General Time Step Given a div free velocity u and pressure p with the acceleration vectors and u the fully implicit algorithm proceeds as follows Incompressible Flow Solver 33 CHAPTER 4 THE PROJECTION METHOD 1 Compute a predicted velocity at t using the known velocity and acceleration Ar Ar Ar 2 i y t cpu aes
33. procedure is presented The parallel assembly procedure then sets the stage for a discussion of the parallel iterative solution methods applied to the PPE Incompressible Flow Solver 21 CHAPTER 3 EXPLICIT TIME INTEGRATION 3 2 1 Domain Decomposition Domain decomposition is the process of sub dividing the spatial domain into sub domains that can be assigned to individual processors for the subsequent solution process There has been a great deal of work on the problem of spatial decomposition for unstructured grids over the last several years 19 20 21 22 23 24 25 26 In general the requirements for mesh decomposition is that the sub domains be generated in such a way that the computational load is uniformly distributed across the available processors and that the inter processor communication is minimized LS DYNA uses recursive coordinate bisection for the decomposition task Details on the use of decomposition tool and decomposition options may be found in the LS DYNA User s Manual 27 In order to exploit the finite element assembly process 8 for parallelization the dual graph of the finite ele ment mesh i e the connectivity of the dual grid shown in Figure 3 1 is used to perform a non overlapping element based domain decomposition Note that the dual grid corresponds to the grid associated with the element centered pressure variables in the Q1Q0 element Implicit in this choice of a domain decomposi tion strategy 1s the idea t
34. row compressed matrix vector multiply There is one final modification to the finite element formulation that derives from the explicit treatment of the advective terms For advection dominated flows it is well known that the use of a backward Euler treatment of the advective terms introduces excessive diffusion Similarly Gresho et al 1 have shown 20 Incompressible Flow Solver 3 2 THE DOMAIN DECOMPOSITION MESSAGE PASSING PARADIGM Table 3 1 Memory requirements and operations counts for a matrix vector multiply for various el ement integration rules stabilization operators and storage schemes Nel number of elements Nnp number of nodes EBE element by element Dimension Element Formulation Total Op 2 D 1 pt h stabilization EBE 1 pt y stabilization EBE 2x2 Quadrature EBE 2x2 ITPACK V 18 1 pt EBE 1 pt h stabilization EBE 1 pt y stabilization EBE 2x2x2 Quadrature EBE 2x2x2 ITPACKV 18 that forward Euler treatment of the advective terms results in negative diffusivity or an under diffusive scheme In order to remedy this problem balancing tensor diffusivity BTD derived from a Taylor series analysis to exactly balance the diffusivity deficit is adopted In the one point quadrature element the BTD term is simply added to the kinematic viscosity in Eq 3 39 to form the tensorial diffusivity used in Eq 3 34 i At fij Mij PZ ui 3 39 In summary the modifications made
35. sense the time average is approximated by N ni n A lt 03 gt lim St Go A 8 2 N 300 pe At 9 where Ar is the time step increment for the nth time integration step and 1 7 Ar is the accumulated time A flow is considered to be non stationary if there is no time invariant mean that can be defined In order to define mean quantities for a non stationary flow an ensemble average is required Here the mean values are computed by averaging over a significant number of flow trials until the mean quantities are independent of the number of trials Formally the ensemble average of x f is Ntrial 2 On xi t t 8 3 lt 0 x T gt lim o E Niras Nial where T t t is the elapsed time from the start of the flow experiment fo and Ntrial is the number of computational experiments In a spatially homogeneous flow the statistical properties of the flow do not vary in space All the mean values are independent of location and can by described by a single value at an instant in time instead of a value at each spatial location The mean values for homogeneous flows are determined from a volume average 1 e or lim Sl xj t d 2 8 4 Q gt Q 66 Incompressible Flow Solver 8 1 DERIVED FLOW STATISTICS
36. test scale and SGS stresses a Lij p uju uju 7 32 qe where L T ij ij Tj j The test scale kinetic energy is defined as the trace of Eq 7 32 i e 1 Ke zk 7 33 Incompressible Flow Solver 63 CHAPTER 7 TURBULENCE MODELS The test scale kinetic energy is analogous to k but it is produced at the large scales and dissipated at the small scales by H un Ou Ou Oii OU 7 34 p ox j ox j ox j ox j Similarly the dissipation is based on the molecular and turbulent viscosity since K represents the re solved kinetic energy at the test scale The SGS stresses Uii and the Leonard stresses L have been observed to exhibit similarities in experi mental data as reported by by Liu et al 97 for the fully developed region of a circular jet This suggests that L be modeled in a manner similar to the SGS stresses that is a 1 Lij 2C AV Re Sij z ijler 7 35 Here C is the only unknown since L may be obtained by application of the test scale filter to the com puted velocities Thus Eq 7 35 provides an explicit representation for C in terms of known quantities However this equation is overdetermined because there are five equations and only a single unknown Cy Using the least squares method introduced by Lilly 96 Lijoij C 7 36 OijOij where Oj Ay kest Sij 7 37 which can be determined entirely from data available at the test scale
37. the mean flow As shown in Figure 7 1 T is the time scale associated with the turbulent fluctuations while 75 is the time scale for the periodic mean flow For this case the decomposition becomes 0 x 1 P x t 0 xit 7 3 accounting for time dependent mean flow variables x t A homogeneous turbulent flow is one for which the flow field is uniform in all directions on the average For homogeneous turbulent flows a spatial average is appropriate where lim gt A xp t dV 7 4 Vw Y For flows that are not stationary an ensemble average is appropriate In this case the average is computed over a significant number of flow trials until the mean quantities are independent of the number of trials Formally the ensemble average of is Ntrial 2 0n x 1 1 7 5 Gi t lim Ntrial gt Ntrial where T t fo is the elapsed time from the start of the flow experiment fo and Ntrial is the number of experimental trials 58 Incompressible Flow Solver 7 1 AVERAGING FILTERING AND SCALE SEPARATION t Figure 7 1 Time averaging for non stationary flow For flows that are both stationary and homogeneous time space and ensemble averaging are all equivalent 1 e the ergodic hypothesis applies However in a practical sense homogeneous turbulent flow is strictly an idealization because boundary conditions such as at no slip no penetration surfaces introduce flow inhomogeneities 7
38. to have a normalized grind time using the iterative solvers that approaches or beats the re solve cost for the PVS direct solver Note that the CPU time associated with the PVS re solve scales roughly as Incompressible Flow Solver 49 CHAPTER 6 PRESSURE SOLUTION METHODS Nel x Nb where Nb is the 1 2 bandwidth of the PPE As Tables 6 1 and 6 2 show there is a small but noticeable effect of the stabilization on the iteration count and grind time for the 2 D jet for both the explicit and semi implicit projection algorithms This is due to the fact that for all computations a convergence criteria of 107 was used As shown in Figure 6 3 this level of convergence limits makes the effect of the stabilization appear reduced However the SSOR preconditioning reduces the iteration count by nearly a factor of 3 without the A conjugate projection The use of the Eisenstat formulation further reduces the grind time by a factor of approximately 1 6 For the explicit algorithm the use of 5 base vectors for the A conjugate projection algorithm reduces the iteration count by nearly a factor of 3 and the grind time by nearly a factor of 3 regardless of the preconditioner For the semi implicit projection algorithm the effect is reduced somewhat resulting in a factor of about 2 25 reduction in the iteration count and better than a factor of 2 reduction in the grind time The combination of the ESSOR PCG algorithm and 10 projection vectors results in grin
39. u Y N 0 u 3 29 a 1 where N 0 indicates evaluation of the shape functions at the origin of the referential coordinate system 18 Incompressible Flow Solver 3 1 MODIFIED FINITE ELEMENT FORMULATION The application of single point integration further simplifies the advection operator In two dimensions 044 and oo in Eq 3 30 are used to form the advection operators in Eq 3 31 where Uab Ug up Qj uC VCy1 0 UCh YCy2 3 30 A u u Af 1 1 1 1 lous 02u24 A A u v 1 1 1 1 louis O2v24 3 31 For the evaluation of the three dimensional advection operators D B4 are defined as Bi ula T YCy WCz1 Bo IC VCy2 WCoa B3 0034 VC3 WCz 3 32 Ba Wy4 VCy4t WCa and used in the element level advection operators P ve A u u L 1 1 1 1 1 1 1 uus Bats Bauzs Bara A u v ye 1 l 1 l L l 1 I Bii B2v28 zn B3v35 Bavae 3 33 A w Ve 1 1 1 1 1 1 1 I gre Bows Bawss Bawas The single point diffusion operator may be stated simply as Kip Ci MijCj Ve 3 34 where i represents a tensorial viscosity Here i and j range from 1 to the number of spatial dimensions while a b range from 1 to the number of nodes per element Generation of the diffusion operator in Eq 3 34 using single point integration leads to rank deficiency of the element level operator The presence of an improper singular mode in the element level operator may also
40. with the continuous projection operators P and Q Eq 4 19 is the consistent discrete form of the elliptic operator for the projection algorithm It represents an algebraic system of equations that is solved for the element centered Lagrange multiplier during the time marching procedure Figure 3 1 shows the dual staggered grid where A and P are centered The pressure Poisson operator in Eq 4 19 incorporates the effect of the essential velocity boundary condi tions and automatically builds in the boundary conditions from Eq 2 15 and 2 17 see Gresho et al 4 The projection algorithm proceeds as follows Given a div free velocity u and its corresponding pressure field p solve the momentum equations for an approximate velocity field at n 1 1 Calculate the approximate velocity field M A O K M Ar 1 O K u Ar Oc F 1 07 F A u u MM Cp 4 20 Here Og and Or control the time weighting applied to the viscous and body force terms with 0 lt 0k 0r lt 1 For Or 1 the viscous terms are treated explicitly while for Og 0 the viscous terms are treated implicitly i e via backward Euler and for Og 1 2 the viscous treatment corresponds to a Crank Nicolson integrator 2 Given the approximate velocity solve Eq 4 19 for A 3 Project the approximate velocity to a div free subspace u Sa M7 CK 4 21 4 After the velocity update
41. 0 wall temperature BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 20 104 1 300 00000 5 outflow boundaries right side top BOUNDARY OUTFLOW CFD SET ssid 35 BOUNDARY OUTFLOW CFD SET ssid 45 BOUNDARY OUTFLOW CFD SET ssid 55 CEEEEPEEREEREEEEEPEREREREEEEEEEEREEEEPEREREEREREEEPEEEREEREREEEPEEEREEEEREEEEEEEEERERS 5 Setup fluid properties 5 MAT CFD CONSTANT mid ro mu k Cp 1 1 1770 1 9832e 5 2 6240e 2 1 0060e 3 gx gy gz diffi diff2 0 29 81 0 diff6 diff7 diff8 diff9 diff10 5 Section amp part 5 SECTION SHELL sid elform shrf nip propt 1 31 tl t2 t3 t4 nloc beta 3 33e 3 diff3 qr irid 1 0000000 1 0000000 1 0000000 1 0000000 0 0000000 0 0000000 Incompressible Flow Solver Tref 300 0 diff4 icomp diff5 123 CHAPTER 10 EXAMPLE PROBLEMS PART Fluid Domain pid secid mid eosid hgid grav adpopt tmid 1 1 1 0 0 0 0 0 5 O O O 5 Handle output of state time history and restart data 5 DATABASE BINARY D3PLOT dt cycl ledt beam npltc 1 0e 2 DATABASE BINARY D3THDT dt cycl ledt beam npltc 1 0e 6 DATABASE HISTORY NODE ndi nd2 nd3 nd4 nd5 nd6 nd7 nd8 11 24 51 83 289 5 DDD Input the mesh node amp segment sets INCLUDE momjet dat k END Example momjet p2 k 124 Incompressible Flow Solver 10 3 THE MOMENTUM DRIVEN JET Segment Set ID 55 Segment Set ID 45 Node Set ID 60 Symmetry Plane Segment Set ID 35
42. 1 4 Spatial Filtering An alternative approach to turbulence modeling is large eddy simulation LES which relies on the direct computation of the largest flow structures with a model for the small scale turbulence The formulation of the LES problem relies on a separation of the flow field into resolved and sub grid fields xi f xit 0 xit 7 6 where 6 is the resolved field and 0 is the sub grid field The resolved field is defined in terms of a spatial filter which for homogeneous turbulence may be defined as unt ENEE de n where G is the filter kernel and Q is the flow domain In order to reproduce constant fields the filter must satisfy the condition that f 66040 1 7 8 Qa where 2 corresponds to the filter support Examples of the Gaussian sharp cut off and top hat filters may be found in 75 83 The filtering procedures used for LES are different than the time averaging discussed above Contrary to the traditional Reynolds time averaging o 0 7 9 Incompressible Flow Solver 59 CHAPTER 7 TURBULENCE MODELS and in general 97 7 10 1 e repeated filtering operations remove additional short wavelength information from the field The support for the filter indicated above by 92 introduces a length scale that corresponds to the smallest turbulence scale A that may be considered to be resolved Thus the spatial filtering procedure provides the mechanism to separate the flow
43. 2 first time step 33 general time step 33 start up 32 projection properties 28 semi implicit projection method 29 Properties fluid properties 97 Restart d3dump 99 d3mean 99 runrsf 99 Sense switches 139 flow solver sense switches 139 Sets node 35 Incompressible Flow Solver INDEX segment 35 shell 35 solid 35 Source terms 35 body forces 38 buoyancy forces 38 Start up procedure explicit time integration 14 fully implicit projection method 32 semi implicit projection method 31 Time integration explicit time integration 13 domain decomposition message passing 21 modified finite element formulation 16 communication costs 25 parallel explicit algorithm 24 reduced integration operators 17 stability 15 start up 14 projection method 27 fully implicit projection method 32 projection properties 28 semi implicit projection method 29 start up 31 32 termination 99 time step control 81 time step selection 81 Turbulence models 57 94 k 8 model 62 averaging 58 dynamic k transport model 63 ensemble averaging 58 66 filtered equations 61 filtering 58 k transport model 62 LES filtering 59 scale separation 58 Smagorinsky subgrid scale model 62 space averaging 58 spatial filtering 59 time averaging 58 66 time averaged equations 60 149
44. 5 However for meshes graded to resolve boundary layers the advective diffusive stability limit usually dictates the time step At lt 3 15 3 1 Modified Finite Element Formulation Several ad hoc modifications are made to the standard Galerkin finite element formulation for the explicit time integration algorithm These modifications include the use of a row sum lumped mass matrix single point Gaussian quadrature balancing tensor diffusivity BTD and hourglass stabilization to damp the spurious zero energy modes known as keystone or hourglass modes A detailed numerical analysis of these modifications is discussed in Gresho et al 1 Before discussing the reduced integration operators a brief overview of the element matrices associated with Eq 3 2 and 3 1 is presented The element level gradient mass advection and diffusion operators are presented in Eq 3 16 3 19 Here N is the element shape function the subscripts a and b range from 1 to the number of nodes per element Nnpe and 1 lt i j lt Ndim Cf 9Na ay 3 16 a Qe Ox 5 f_ NaMaV 3 17 16 Incompressible Flow Solver 3 1 MODIFIED FINITE ELEMENT FORMULATION aN e u N utay 3 18 ab U de u A i ON ON KS E H i dV 3 19 ab Qe x Ox 3 1 1 Reduced Integration Operators From an examination of the explicit algorithm it is clear that the bulk of the computational effort for a discrete time step consists o
45. 5 Update the base vectors to include new information from the last solution As an aside the solution at a given time level may be written as N x A V 005 6 13 i l Updating For the initial solution or when the number of existing base vectors exceeds N the basis is started by normalizing the solution as xl d1 ben 6 14 For all other cases a solution vector is a candidate for addition to only when it contains non trivial information not already present in Here the basic idea is that Ax is A conjugate to x as well as to the individual base vectors Therefore the addition of a new base vector should be based on criteria that guarantees that new information is being added to the existing base vectors Addition of a solution vector 46 Incompressible Flow Solver 6 4 PPE SOLVER PERFORMANCE or a part of a solution vector proceeds by first computing the part of the solution that is not contained in the basis as l Yi x M adi 6 15 i l Now each 0 is A conjugate to y 1 by construction and is added to the basis as Wi 1 1 1 6 16 lyla The idea of seeding the A conjugate vectors in with several initial vectors has been suggested by T J R Hughes 67 as a mechanism for reducing cost of the initial linear solves i e the first PPE solve when is empty This may be achieved by selecting several global polynomials and even random vector as seeds for Construction o
46. 5 95 559952522999559555599522599555599595299599299 O Setup fluid properties MAT CFD CONSTANT mid ro mu k Cp beta Tref 1 1000 0 1 0e 3 gx gy gz diff1 diff2 diff3 diff4 diff5 diff6 diff7 diff8 diff9 diff10 Section amp part SECTION SOLID 1 31 PART Fluid domain 1 1 1 0 0 0 0 0 5 O O O O 5 Handle output of state time history and restart data 5 DATABASE BINARY D3PLOT 5 0e 4 DATABASE BINARY D3THDT 1 0e 6 DATABASE HISTORY NODE 134 Incompressible Flow Solver 10 4 THE SHEATH FLOW CHAMBER ndi nd2 nd3 nd4 nd5 nd6 nd7 nd8 17400 17427 17397 18989 35693 6842 25474 9 DATABASE BINARY RUNRSF 2000 5 PAPA 5 Input the mesh amp node sets S INCLUDE sheath dat k END Example sheath p2 k z PROTOTYPE SHEATH FLOW CHAMBER 0 025 O 0 02 2 e a o o A Kinetic Energy E 06 0 005 0 0 002 0 004 0 006 0 008 0 01 0 012 Time Figure 10 26 Time history plot of kinetic energy Incompressible Flow Solver 135 CHAPTER 10 EXAMPLE PROBLEMS PROTOTYPE SHEATH FLOW CHAMBER Node No F B B 9 17400 X velocity a 0 0 002 0 004 0 006 0 008 0 01 0 012 hued Time 0002 62 11 00 Sp z4 Figure 10 27 Time history plot of x velocity
47. 5255 INCLUDE duct dat k END Example duct p2 k The grid CFL and Reynolds numbers were reported every 100 time steps by setting ickdt 100 on the CONTROL CFD GENERAL keyword Figure 10 2 shows a sample of the screen output with the CFL limited time step and associated grid parameters When the grid parameters are reported they are pre sented in terms of the xi grid eta grid and zeta grid directions which correspond to the element local n and G coordinate directions as described in Chapter 3 For each local coordinate direction a minimum and maximum grid CFL and Reynolds number is reported along with the corresponding ele ment id In this example the time step is controlled so that the maximum grid CFL number never exceeds 2 0 During the calculation the print sense switch toggles the output of the minimum maximum pres sure and velocities with the RMS divergence at each time step Sample screen output of the reported minimum maximum pressure and velocities with the RMS divergence is shown in Figure 10 3 The fully implicit treatment of the viscous terms is based on a first order Forward Euler time integrator that in combination with a lumped mass matrix yields a rudimentary marching scheme for obtaining a steady state solution to the channel flow problem In order to verify that a steady state has been achieved nodal time histories of the x velocity along the centerline of the channel at the inlet and outlet boundaries are shown
48. 8 Incompressible Flow Solver 67 CHAPTER 8 FLOW STATISTICS Similarly the averaging procedure is commutative and averaging with a fluctuating component yields Qd lt y gt gt lt y lt o gt gt 0 8 9 Let represent a scalar quantity and let u i 1 2 3 represent the ith Cartesian component of the instantaneous velocity For a scalar quantity the variance is o lt osos 8 10 and the root mean square RMS value is c The triple correlation between three independent scalars 0 y and is QA lt owAa gt lt o gt lt y gt lt A gt 8 11 lt o gt lt WV gt 4 lt y gt lt o0N gt 4 lt A gt lt y s Correlations of this type appear in second order moment closure models such as the Reynolds stress turbulence model Higher order correlations can be defined in a similar manner The Reynolds stress tensor is defined in a similar way Rij lt UU gt lt ujuj gt uj gt lt uj 8 12 where u are the velocity components The turbulent kinetic energy is 1 1 T 5 ui FoU ukuk gt lt uk gt lt uy gt 8 13 The turbulent scalar flux vector is Qu gt lt du gt lt 0 gt lt u gt 8 14 8 1 3 Higher Order Statistics Higher order moments of velocity differences and velocity derivatives are frequently used to analyze isotropic turbulence The skewness is defined as uj Soi m 8 15 and the flatness factor is uj Foi 5 8 16 NEU uj The skew
49. 9 1 7 CONTROL CFD PRESSURE Purpose Set the pressure solver parameters to be used for the incompressible Navier Stokes equations Card Format Naxiwisqereraltlw fm ef oe paella pa pete 7 Ae EP EP PPP eee fw foe few me E pel EA VARIABLE DESCRIPTION IPSOL Set the pressure solver type EQ 0 IPSOL 22 for serial IPSOL 21 for MPP default EQ 10 Sparse direct solver EQ 11 PVS direct solver EQ 20 Jacobi preconditioned conjugate gradient method EQ 21 SSOR preconditioned conjugate gradient method EQ 22 SSOR preconditioned conjugate gradient using the Eisenstat transformation 88 Incompressible Flow Solver VARIABLE MAXIT ICHKIT IWRT IHIST EPS NVEC ISTAB BETA SID PLEV LCID 9 1 INCOMPRESSIBLE FLOW KEYWORDS DESCRIPTION Set the maximum number of iterations for the pressure solver EQ 0 MAXIT 200 default Set the interval to check the convergence criteria for the pressure solver EQ 0 ICHKIT 5 default Activate the output of diagnostic information from the pressure solver EQ 0 Diagnostic information is off default EQ 1 Diagnostic information is on NOTE during execution sense switch print can be used to toggle this flag on or off Activate the generation of a convergence history file for the pressure solver The ASCII history file is ppe his EQ 0 Convergence history is off default EQ 1 Convergence history is on Set the convergence crit
50. DESCRIPTION typeID Node ID NID or node set ID NSID see SET_NODE DOF Applicable degrees of freedom EQ 101 x velocity EQ 102 y velocity EQ 103 z velocity EQ 104 temperature EQ 107 turbulent kinetic energy EQ 110 turbulent eddy viscosity EQ 121 species mass fraction 1 EQ 122 species mass fraction 2 EQ 123 species mass fraction 3 EQ 124 species mass fraction 4 EQ 125 species mass fraction 5 EQ 126 species mass fraction 6 EQ 127 species mass fraction 7 EQ 128 species mass fraction 8 EQ 129 species mass fraction 9 EQ 130 species mass fraction 10 EQ 201 x y z velocity EQ 202 x y velocity EQ 203 y z velocity EQ 301 all species LCID Load curve ID to describe motion value versus time see DEFINE CURVE SF Load curve scale factor default 1 0 Example Input 993399999599999999 99 99 99 ooo SoS So Sooo SSS SS SSS SS SS SSS SS SSS SSSSEPUP OT USUS SS BOUNDARY PRESCRIBED CFD SET 5 AAA AAA 5 A set of nodes is given a prescribed velocity in the x direction according to a specified vel time curve which is scaled 5 BOUNDARY PRESCRIBED CFD SET 5 A o O nsid dof lcid sf 4 101 8 2 0 5 nsid 4 nodal set ID number requires a SET_NODE option dof 101 x velocity is prescribed leid 8 velocity follows load curve 8 requires a DEFINE CURVE sf 2 0 velocity specified by load curve is scaled by 2 0 5 O O O 78 Incompressibl
51. EHG 92 DESCRIPTION Time weighting for the diffusion terms Valid values are 0 Ox 1 with Ox 1 2 for second order accuracy in time EQ 0 THETAK 0 5 default Time weighting for advection terms Time weighting for body forces and boundary conditions Valid values are 0 lt 0r lt 1 with Or 1 2 for with for second order accuracy in time EQ 0 THETAF 0 5 default Set the equation solver type for the momentum equations EQ 0 MSOL 20 default EQ 20 Jacobi preconditioned conjugate gradient method EQ 30 Jacobi preconditioned conjugate gradient squared method default when IADVEC 40 Set the maximum number of iterations for the iterative equation solver EQ 0 MAXIT 100 default Set the interval to check the convergence criteria for the iterative equation solver EQ 0 ICHKIT 2 default Activate the output of diagnostic information from the equation solver EQ 0 Diagnostic information is off default EQ 1 Diagnostic information is on Activate the generation of a convergence history file from the equation solver EQ 0 Convergence history is off default EQ 1 Convergence history is on Set the convergence criteria for the iterative equation solver EQ 0 EPS 1 0e 5 default Set the type of hourglass stabilization to be used with the momentum equations This only applies to the explicit treatment of the momentum equations INSOL 1 EQ 0 IHG 1 default EQ 1 LS DYNA CFD viscous h
52. Global Global Local 1 E E 4 2 1 2 3 4 3 3 1 3 a Sequential Mesh 4 5 8 4 Sub domain P1 5 6 9 lt 9 13191 1400 hisp J16 121 623 LL 6 7 9 A 7 f 8 10 12 8 905 1006 11171 1218 NN iden d P H 14 10 9 7 fos Juo 7181 8 A 16 12 rd z Y 5 4 6 5 7 6 3 1 4 2 c Parallel Assembly Mapping Y 111 2 2 3 3 Sub domain PO b 2 Processor Sub domain Mapping Figure 3 4 Parallel assembly procedure a Sequential Mesh b 2 Processor sub domain mapping c Parallel assembly mapping where A is the assembly operator Here the diffusive and advective contributions are computed at the element level using un rolled right to left matrix vector multiplies In order to exploit both register to register vector and cache based processors data independent groups of elements are identified This not only permits the vectorization of the assembly process but also maximizes the number of element level operations with respect to the number of load store operations For cache based architectures this permits all of the element data in a data independent group to be loaded into cache once for the element level operations e g right to left matrix vector multiplies Thus the assembly process using the vector cache blocks proceeds block by block with all element level operations for a block being performed before completing the add scatter portion of the gather add s
53. LS DYNA s Incompressible Flow Solver User s Manual MARK A CHRISTON Ph D and GRANT O COOK JR Ph D March 2001 Version 960 Development version Copyright 2001 LIVERMORE SOFTWARE TECHNOLOGY CORPORATION All Rights Reserved Mailing address Livermore Software Technology Corporation 2876 Waverley Way Livermore California 94550 1740 Support Address Livermore Software Technology Corporation 7374 Las Positas Road Livermore California 94550 FAX 925 449 2507 TEL 925 449 2500 EMAIL sales lstc com Copyright 2001 by Livermore Software Technology Corporation All Rights Reserved Contents 1 Introduction 1 1 Guide to the User s Manual a us a see ee ee m ee Se eS SR de 2 Governing Equations 2 1 Conservation Equations a xe eo ee Ll Ro ele eh he a A Bo ee BP oe 2 1 1 Momentum Conservation 2a ae Gig rece edo audi go ROSE So a eed 2 1 2 Mass Conservation 4 uui due Madd RA Be Be ee GA aera AER 2 1 3 Energy Conservation oe 6 59 v ee ee ioe be he ese hs 2 1 4 Boundary and Initial Conditions llle 2 2 Conservation Equations Vector Notation o e e 2 2 1 Momentum Conservation uw oe en dR eee ee ewe Ry RR 222 Mass Conservation v ug kone ah ey Se ch DA Bk 2 2 3 Energy Conservation s 23 ose QR Pe ue Ex Se be a as 22 4 Boundary and Initial Conditions less 3 Explicit Time Integration 3 1 Modified Finite Element Formulation ss RR m Ree 3 1 1 Red
54. OBLEMS Prescribed inlet BC BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 3 102 2 1 0000000 Setup the no slip BC s BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 101 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 102 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 2 101 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 2 102 2 1 0000000 59 9 5 9 559952259995595555995252299555559525299555299959599959 92999 9 9 9 9 9 9 99OS Setup fluid properties 5 MAT CFD CONSTANT mid ro mu k Cp beta Tref 1 1 0000 1 0e 2 gx gy gz diff1 diff2 diff3 diff4 diff5 diff6 diff7 diff8 diff9 diff10 5 Section amp part 5 SECTION SHELL sid elform shrf nip propt qr irid icomp 1 31 tl t2 t3 t4 nloc 1 0000000 1 0000000 1 0000000 1 0000000 PART Fluid domain pid secid mid eosid hgid grav adpopt tmid 1 1 1 0 0 0 0 0 5 AAA AAA 5 Handle output of state time history and restart data DATABASE BINARY D3PLOT dt eyel ledt beam npltc 104 Incompressible Flow Solver 10 1 POISEUILLE FLOW 10 0 0 DATABASE BINARY D3THDT dt cycl ledt beam npltc 1 0e 6 DATABASE HISTORY NODE ndl nd2 nd3 nd4 nd5 nd6 nd7 nd8 6 28 39 50 215 226 DATABASE BINARY RUNRSF dt oyel ledt beam npltc 1000 DATABASE BINARY D3DUMP dt cycl ledt beam npltc 0 99 552555355252525225552225323525252525252255522252555252525252253225552225225522555
55. RMANCE Table 6 4 Effect of the A Conjugate projection on the PPE solution time for the 3 D square cylinder vortex shedding problem The grind times for the projection algorithm have been normalized with respect to the solution time for 1000 time steps using a direct solver 71 72 Nr is the average number of iterations required per time step to solve the PPE problem Nel 22380 Nnp 25568 N 772 Projection 2 Algorithm Global Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG RT I ROO No of Grind Grind Grind Vectors Pur ame x LI Nir T a Nur Time 2 Algorithm Local Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG No Ns Grind Grind Grind Vectors Mir ame ud HE DIT ame T 25 Ar in i s Incompressible Flow Solver 55 CHAPTER 6 PRESSURE SOLUTION METHODS 56 Incompressible Flow Solver Chapter 7 Turbulence Models Turbulent flows occur in a wide variety of engineering applications and they are characterized by their inherent unpredictability enhanced mixing properties and broad range of length and time scales The length scales in a turbulent flow are bounded by the bulk or integral length scale e g a pipe diameter or body dimension and by the dissipation scale where physical viscosity acts Similarly the time scales are bounded by the time scale associated with the turn over time of the largest eddies in the flow and the time scale for viscous dissipation Resolution of all the length and t
56. Reynolds stress tensor the SGS stress tensor is symmetric and introduces new unknowns which must be represented via a model The definition of a model for the subgrid scale stress tensor permits the solution of the filtered Navier Stokes equations for the resolved scale velocity 7 and pressure p However like the time averaged equations closure of the spatially filtered LES equations can vary in the degree of complexity and success in simulating turbulent flow The most simple and most common prescription for the SGS stresses relies on a parameterization in terms of the strain rate and an eddy viscosity 1 e the Boussinesq hypothesis The Boussinesq hypothesis assumes that the principle axes of the SGS stress tensor 1 are aligned with the principle axes of the filtered strain rate tensor Thus the relationship between the deviatoric SGS stress and strain rate is A 1 e ti i GT 2urSij 7 21 where u is a turbulent eddy viscosity The parameterization of the deviatoric SGS stresses in terms of the strain rate permits the filtered momen tum equations to be written as Qu Ou Q Pp Q U 4 32 Si Ff 7 22 py pu Ox 2 aed u u Sij pfi where 1 p pt E 7 23 Incompressible Flow Solver 61 CHAPTER 7 TURBULENCE MODELS 7 3 1 Baseline Smagorinsky Subgrid Scale Model One of the most common models for the LES subgrid scale stresses was first proposed by Smagorinsky 85 and uses the Boussinesq hypoth
57. S 7 2 The Time Averaged Equations 2 aie UR EUR OUR ee Re x oe 7 3 The Spatially Filtered Equations 2 a 7 3 1 Baseline Smagorinsky Subgrid Scale Model lll 7 4 Subgrid Kinetic Energy Models oe eu qe qe e qe ee pe fa npe Roe Ter apt 7 4 1 Local Dynamic JE model us Ium vem ivan SUR SCR MES Flow Statistics 8 1 Derived Flow Statistics 5 ss o oe de M et edo ih a ih das de 8 1 1 Reynolds Averaged Statistics loko Re REOR ex OR bate Paes 8 1 2 Mean and Fluctuating Quantities 220 dei ee eek ec De a ete Tena 8 1 3 Higher Order Statistics ae Led od bn dole deg do Deb ed Se 4 8 8 1 4 Anisotropic Stress Tensor LAA Al e ius as RA fan AR ge o at rs Ge n 82 LS DYNA Statistics Levels gt s sce vue X ees Rea Bad C ATO us eR as SS S D Level Mean Statistics e uou m ue oh k ee Pe AA CUR ORE WO 8 2 2 Second Moment Statistics le be Sb ee ee Bene Sod Eus 8 2 3 Higher Order Statistics o in te ano se eee ee ne ee BE eue 8 3 LS POST Statistics Levels uu dox oC ees eos AS ang RA Keyword Input 9 1 Incompressible Flow Keywords 2o uo be Se UR RUMOR EUR MEE OS x 9 1 1 BOUNDARY OUTFLOW CFD OPTION lee 9 1 2 BOUNDARY PRESCRIBED CFD OPTION 9 3 BOUNDARY PRESSURE CFDSET 57 58 59 60 61 62 62 63 65 65 66 67 68 69 70 71 71 12 73 75 9 1 4 CONTROL CFD_AUTO 5 0 LS EE RRS RAE RRS URS 81 9 1 5 CONTROL CFD
58. Silvester Stabilised vs stable mixed methods for incompressible flow Technical Report NA 307 University of Manchester UMIST Manchester Centre for Computational Mathematics September 1997 144 Incompressible Flow Solver BIBLIOGRAPHY 62 Stevens T Chan and Gayle Sugiyama User s manual for mc_wind A new mass consistent wind model for arac 3 Technical Report UCRL MA 129067 Lawrence Livermore National Laboratory Livermore California 94550 November 1997 63 Mark A Christon A domain decomposition message passing approach to transient viscous in compressible flow using explicit time integration Computer Methods in Applied Mechanics and Engineering 148 329 352 1997 Sandia SAND96 2566J 64 Paul F Fischer Projection techniques for iterative solution of ax b with successive right hand sides Technical Report 93 90 ICASE 1993 Also NASA CR 191571 65 Youcef Saad On the lanczos method for solving symmetric linear systems with several right hand sides Mathematics of Computation 48 178 65 1 662 April 1987 66 Tony F Chan and W L Wan Analysis of projection methods for solving linear systems with multiple right hand sides SIAM Journal on Scientific Computing 18 6 1698 1721 November 1997 67 T J R Hughes personal communication October 1998 68 S C Eisenstat Efficient implementation of a class of preconditioned conjugate gradient methods SIAM Journal for Scientific and Statistical Comput
59. YNA control file for the mesh with 441 nodes is shown below for the semi implicit projection method box p2 k Similar to the duct p2 k example the box p2 k uses a fixed CFL 2 condition with the automatic time step selection and an initial time step of dtinit 1 0e 4 In order to solve an energy equation in conjunction with the momentum equations itemp 1 on the CONTROL CFD TRANSPORT keyword In addition boundary conditions specifying 0 1 on the left wall node set 2 and 0 0 on the right wall node set 3 are included in the keyword input In this example the fluid properties have been translated into non dimensional parameters i e u Pr D RaPr with the inclusion of the gravity vector gy 1 and a reference temperature 0 27 0 5 As a basis for comparison the explicit semi implicit and fully implicit methods have been applied to this problem For the explicit method a fixed time step was used with At 2 0e 4 This requires iauto 1 on the CONTROL CFD AUTO keyword The semi implicit method started with a time step At 1 0e 4 and increased the time step according to the maximum CFL 0 5 condition iauto 2 on the CONTROL_CFD_AUTO keyword The fully implicit method used a second order Adams Bashforth Incompressible Flow Solver 111 CHAPTER 10 EXAMPLE PROBLEMS predictor with a trapezoid rule corrector with the local time truncation error specified as epsdt 0 001 and iauto 3 on the CONTROL CFD AUTO keyword
60. additional solvability constraint as shown in Eq 2 27 Ou Fr 0 2 25 nju x 0 niu x 2 26 nu dr 0 2 27 E 2 2 Conservation Equations Vector Notation In the ensuing discussion the invariant bold face vector notation of Gibbs is used with boldface symbols representing vector tensor quantities The reader may refer to Gresho and Sani 5 pp 357 359 for an outline of notation for the Navier Stokes equations 2 2 1 Momentum Conservation To begin the conservation of linear momentum is O V o pf 2 28 where u u v w the velocity o is the stress tensor p the mass density and f the body force The body force contribution pf typically accounts for buoyancy forces with f representing the acceleration due to gravity The stress may be written in terms of the fluid pressure and the deviatoric stress tensor as 6 pl t 2 29 Incompressible Flow Solver y CHAPTER 2 GOVERNING EQUATIONS where p is the pressure I is the identity tensor and T is the deviatoric stress tensor A constitutive equation relates the deviatoric stress and the strain rate e g T 2S 2 30 Here the dynamic viscosity u is a second rank tensor Frequently the fluid viscosity is only available as an isotropic viscosity in which case the constitutive relation becomes T 2p8 2 31 The strain rate tensor is written in terms of the velocity gradients as l T S Vu Vu 2 32 2 20 0 Mass Conservation
61. algorithm INSOL 3 84 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 6 CONTROL CFD MOMENTUM Purpose Set the solver parameters to be used for the momentum equations in the Navier Stokes solver Card 1 is used to control the time integrator and advective transport options Card 2 is used to set the linear solver options such as the maximum iteration count and interval to check the convergence criteria Card Format pwr a a fos fae fos fie or es Variable IMASS IADVEC IFCT DIVU THETAK THETAA THETAF EM eet pe pete fe fe fe Remarks pew fe fate fs Pe ffs Pe sm oe oe VARIABLE DESCRIPTION IMASS Select the mass matrix formulation to use EQ 0 IMASS 1 default EQ 1 Lumped mass matrix EQ 2 Consistent mass matrix EQ 3 Higher order mass matrix Incompressible Flow Solver 85 CHAPTER 9 KEYWORD INPUT VARIABLE IADVEC IFCT DIVU THETAK THETAA THETAF MSOL MAXIT ICHKIT IWRT 86 DESCRIPTION Toggle the treatment of advection between explicit with balancing tensor diffusivity or fully implicit EQ 0 IADVEC 10 for forward Euler with BTD default EQ 1 IADVEC O for forward Euler without BTD EQ 10 forward Euler with BTD EQ 40 fully implicit simplified trapezoid rule Toggle the use of the advective flux limiting advection scheme EQ 0 IFCT 1 default EQ 1 Advective flux limiting is on EQ 1 Advective flux limiting is off Set the RMS div
62. an updated pressure at time level n 1 is obtained via atl AG 4 22 p p Re 4 22 30 Incompressible Flow Solver 4 1 SEMI IMPLICIT PROJECTION METHOD Remark In the projection algorithm the start up procedure to obtain the pressure p proceeds as follows 1 Calculate the partial acceleration 1 e acceleration neglecting the pressure gradient at time level n by solving the mass matrix problem M F Ku A u u 4 23 where is the instantaneous acceleration neglecting the pressure gradient 2 Solve the global PPE for the current pressure field C7 M c p C73 gl 4 24 In LS DYNA the prescribed initial conditions and boundary conditions are tested and if necessary a projection to a div free subspace is performed on the initial velocity field u This guarantees that the flow problem is well posed even if the user prescribed initial conditions violate the conditions of Eq 2 25 2 26 In practice the criterion for performing a div free projection is based upon the RMS diver gence error CTu CTu m 4 25 Nel S where Nel is the number of elements and e is a user specified tolerance typically 1071 to 1077 If the RMS divergence error is greater than the specified tolerance for the initial can didate velocity field a then the PPE problem in Eq 4 26 is solved for A and a mass consistent projection performed using Eq 4 27 c M cjA C a 2 4 26
63. ariable according to a prescribed value or function of time All nodal boundary conditions for the incompressible flow solver are prescribed using the BOUND ARY PRESCRIBED CFD NODE and BOUNDARY_PRESCRIBED_CFD_SET keywords As an example of nodal boundary conditions consider the prescription of boundary conditions based on the flow domain shown in Figure 5 1 Inflow boundary velocity conditions that emulate free stream conditions are prescribed as uj U and u2 0 on I No slip and no penetration velocity boundary conditions are prescribed at the cylinder wall u 0 on Ts 35 CHAPTER 5 BOUNDARY CONDITIONS AND SOURCE TERMS Outflow Boundary Inflow Boundar Segment Set I Node Set E U E Y ONT E SOKO PK SSS SS SSRI ae 25 No Slip No Penetration Node Set I Figure 5 1 Flow domain with inflow I outflow I and no slip no penetration I boundaries 36 Incompressible Flow Solver 5 3 TRACTION BOUNDARY CONDITIONS The boundary conditions for the energy equation Eq 2 13 consist of a prescribed temperature or heat flux rate For example a prescribed temperature at the cylinder wall is T x t f on T Turning attention to the species transport equations boundary conditions for Eq 2 9 may consist of either a prescribed concentration or a mass flux rate In the case of a binary mixture with a prescribed free stream concentration Z1 xj t Z on T where is the known value of co
64. ary conditions from Eq 2 15 and automatically builds in the boundary conditions from Eq 2 17 see Gresho et al 6 Equations 3 1 and 3 3 form the basis for a description of the explicit time integration algorithm It is assumed that the explicit algorithm begins with a given divergence free velocity field u that satisfies the essential velocity boundary conditions and an initial pressure p To simplify the description of the 13 CHAPTER 3 EXPLICIT TIME INTEGRATION explicit algorithm velocity boundary conditions that are constant in time are assumed making g 0 see Gresho and Sani 10 for additional algorithmic issues for the situation when g 4 0 The explicit algorithm proceeds as follows 1 Calculate the partial acceleration i e acceleration neglecting the pressure gradient at time level n m 1 a M F 3 4 where 7 F F Ku A u u 3 5 2 Solve the global PPE for the current pressure field Ic M c p C a 3 6 3 Update the nodal velocities u u Ar 8 Mj Cp 3 7 4 Repeat steps 1 3 until a maximum simulation time limit or maximum number of time steps is reached Primary Grid 29 F e F e L o Pressure MS N vasiy DOF ooo i Dual Grid 77 id ual Gri duos dilo Figure 3 1 Velocity mesh with two degrees of freedom DOF per node and the PPE dual grid with one DOF per element Remark In LS DYNA the pre
65. at the centerline of the inlet outlet planes PROTOTYPE SHEATH FLOW CHAMBER Node No _A 6842 0 6 B NL B B _B 25474 0 4 0 2 E 07 o 9 L o 1 gt 0 24 0 4 A A A A A 0 6 Y 0 8 1 1 1 1 o 0 0 002 0 004 0 006 0 008 0 01 0 012 El 3 Time 0002 62 11 oziz Figure 10 28 Time history plot of y velocity at the centerline of the side inlet planes 136 Incompressible Flow Solver 10 4 THE SHEATH FLOW CHAMBER PROTOTYPE SHEATH FLOW CHAMBER Time 0 013068 Fringe Levels Isosurfaces of Pressure min 1422 76 at nodes 17463 max 10943 9 at node 20594 1 094e 04 9 886e 03 8 828e 03 7 770e403 6 712e403 5 654e 03 4 596e 03 3 539e 03 2 481e 03 1 423e403 Figure 10 29 Pressure isosurfaces at t 0 013 s PROTOTYPE SHEATH FLOW CHAMBER Time 0 013068 Isosurfaces of Helicity min 49474 9 at nodes 32356 max 49654 9 at node 30022 Fringe Levels 4 965e 04 3 064e 04 2 763e 04 1 661e 04 5 597e 03 5 417e 03 1 643e 04 2 745e 04 3 846e 04 4 947e 04 Figure 10 30 Helicity isosurfaces at t 0 013 s Incompressible Flow Solver 137 CHAPTER 10 EXAMPLE PROBLEMS 138 Incompressible Flow Solver Appendix A Sense Switch Controls The status of an LS DYNA calculation can be determined by typing C control c This sends an interrupt to LS DYNA which is trapped and the user is prompted to input a
66. ated cavity The non dimensional governing equations for the time dependent thermal convection problem in vector form are the incompressible Navier Stokes equations conservation of mass and the energy equation written in terms of temperature ou t Vu VP PrV u RaPr jO 10 2 V u 0 10 3 and 30 7 u V0 V 9 10 4 where u u v P and 0 are the velocity the deviation from hydrostatic pressure and temperature re spectively and j the unit vector in the y direction These non dimensional equations were obtained using the characteristic length L velocity U o L time scale t Q IE and pressure P pU Here p is the mass density g the gravitational acceleration amp k pC is the thermal diffusivity and v is the kinematic viscosity The Prandtl number is Pr v and fixed at Pr 0 71 and the Rayleigh number is T TeL Ra 8S TD 10 5 va where 7 T the temperature difference between the hot and cold walls and B the coefficient of thermal expansion The non dimensional temperature is defined in terms of the wall temperature difference and a reference temperature as iT Tp TE where T is the prescribed temperature of the hot wall and T is the temperature of the cold wall 10 6 The boundary conditions for this problem consist of no slip and no penetration walls with the top and bottom walls insulated The initial conditions are prescribed with u v 0 and T T T 2 The LS D
67. cal solution International Journal for Numerical Methods in Fluids 3 249 264 1983 Incompressible Flow Solver 147 Index Bibliography 147 Boundary conditions 5 9 35 flux 39 heat flux rate 39 mass flux rate 39 natural 37 nodal 35 77 node sets 35 outflow 38 76 pressure 37 79 segment sets 35 time dependent 99 traction 37 79 Conservation equations energy equation 90 momentum equations 85 species transport 90 vector notation 7 Element formulation 99 modified finite element formulation 16 reduced integration operators 17 Example problems 101 2 D channel flow 101 momentum driven jet 121 natural convection 111 Poiseuille flow 101 sheath flow chamber 131 thermal cavity 111 Flow solver fluid structure interaction 90 stand alone 99 Flow statistics 65 anisotropic stress tensor 69 derived flow statistics 65 fluctuating quantities 67 higher order statistics 68 mean quantities 67 Reynolds averaged statistics 66 statistics levels 70 higher order statistics level 3 72 mean statistics level 1 71 second moment statistics level 2 71 Fluid structure interaction 99 Governing equations 3 energy conservation 5 9 mass conservation 4 8 momentum conservation 3 7 vector notation 7 Grid parameters CFL number 15 105 Reynolds number 15 105 Initial conditions 5 9 95 Introduction 1 guide to the user s manual 2 Keyword input 75 incompressible flow keywords
68. catter assembly procedure The parallel assembly procedure may be viewed as a generalized form of the finite element assembly algorithm However inter processor communication is an inherent part of the parallel assembly As an example of a two processor assembly consider the sequential mesh and the assignment of the global nodes to two processors as shown in Figure 3 4 In Figure 3 4b the local node numbers are enclosed in brackets to the right of the global node number global node 1 in sub domain PO is local node 1 The sub domain assignment of the global nodes is the consequence of the unique assignment of elements to processors and reveals the existence of the nodes at the sub domain boundaries on multiple processors Figure 3 4b shows the global inter processor assembly of the sub domain boundary nodes Figure 3 4c Incompressible Flow Solver 23 CHAPTER 3 EXPLICIT TIME INTEGRATION illustrates the use of the local to global mapping required for the gather add scatter operation and the arrows between global node numbers identify a send receive pair Thus the parallel assembly procedure induces communication in the form of a gather send receive add scatter process The parallel right hand side assembly for Eq 3 5 may be viewed as a generalized assembly with off processor communication where first the on process vector cache blocked assembly is performed according to Eq 3 41 and then the assembly at the sub domain boundaries is perf
69. ce criteria Card Format pow a fos fe fos fe fade Remarks P e ale ee ET Remarks 90 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS pews fe fe fe fs fe fle Variable ITSOL MAXIT ICHKIT IWRT IHIST EPS IHG EHG Remarks VARIABLE DESCRIPTION ITEMP Solve the energy equation in terms of temperature EQ 0 No energy equation default EQ 1 Energy equation is solved in terms of temperature NSPEC Activate the solution of NSPEC species transport equations EQ 0 No species equations are solved default EQ NSPEC Solve for NSPEC species Up to 10 species transport equations may be active 0 lt NSPEC lt 10 IMASS Select the mass matrix formulation to use EQ 0 IMASS 1 default EQ 1 Lumped mass matrix EQ 2 Consistent mass matrix EQ 3 Higher order mass matrix IADVEC Toggle the treatment of advection between explicit with balancing tensor diffusivity or fully implicit EQ 0 IADVEC 10 for forward Euler with BTD default EQ 1 IADVEC 0 for forward Euler without BTD EQ 10 forward Euler with BTD EQ 40 fully implicit with simplified trapezoid rule IFCT Toggle the use of the advective flux limiting advection scheme EQ 0 IFCT 1 default EQ 1 Advective flux limiting is on EQ 1 Advective flux limiting is off Incompressible Flow Solver 91 CHAPTER 9 KEYWORD INPUT VARIABLE THETAK THETAA THETAF ITSOL MAXIT ICHKIT IWRT IHIST EPS IHG
70. d heat flux rate is kVT n Q x t 2 49 where q is the known flux rate through the boundary with normal n The heat flux rate may also be prescribed in terms of a heat transfer coefficient kVT n h T To 2 50 where h is the heat transfer coefficient and To is a reference temperature Initial conditions take on the form of prescribed velocity species and temperature distributions at t 0 i e Zi 2 51 Remark For a well posed incompressible flow problem the prescribed initial velocity field in equation 2 52 must satisfy equations 2 52 and 2 53 see Gresho and Sani 6 If T2 0 the null set i e enclosure flows with u prescribed on all surfaces then global mass conservation enters as an additional solvability constraint as shown in equation 2 54 V u 0 2 52 10 Incompressible Flow Solver 2 2 CONSERVATION EQUATIONS VECTOR NOTATION n u x 0 n u x 2 53 i nud 0 2 54 Tr Incompressible Flow Solver 11 CHAPTER 2 GOVERNING EQUATIONS 12 Incompressible Flow Solver Chapter 3 Explicit Time Integration This chapter presents the spatial discretization and explicit time integration method for the incompressible Navier Stokes equations The spatial discretization is achieved using the Q1Q0 element with bilinear sup port for velocity and piecewise constant support for the pressure in two dimensions In three dimensions the velocity support is trilinear with piecewise
71. d structure interaction problems where the flow is in the low Mach incompressible regime LS DYNA s incompressible flow solver is based in part upon the work of Gresho et al 1 2 3 4 and makes use of advanced solution algorithms for both implicit and explicit time integration The explicit solution algorithm 1 2 sacrifices some phase accuracy but decouples the momentum equations and min imizes the memory requirements While both the diffusive and Courant Freidrichs Levy CFL stability limits must be respected in the explicit algorithm balancing tensor diffusivity meliorates the restrictive diffusive stability limit and raises the order of accuracy of the time integration scheme The explicit al gorithm in combination with single point integration and hourglass stabilization has proven to be both simple and efficient in a computational sense In the second order projection algorithm 3 4 a consistent mass predictor in conjunction with a lumped mass corrector legitimately decouples the velocity and pressure fields thereby reducing both memory and CPU requirements relative to more traditional fully coupled solution strategies for the Navier Stokes equations The consistent mass predictor retains fourth order advective phase accuracy while the lumped mass corrector a projection to a divergence free subspace maintains a divergence free velocity field Both the predictor and the corrector steps are amenable to solution via direct or pr
72. d times that are within a factor of 3 of the grind time using the PVS solver For the projection algorithm the grind times are within a factor about 2 emphasizing the additional cost of the semi implicit treatment of the momentum equations Relative to the basic JPCG solver the ESSOR PCG solver with 10 projection vectors results in a reduction of nearly 9 for the explicit algorithm and nearly 4 for the projection algorithm Interestingly the storage cost of a separate preconditioning matrix is approximately equivalent to 10 projection vectors Although further reductions in the execution times may be had with increasing number of projection vectors the increased storage cost and diminishing reduction in the computational costs makes suggests that 5 to 10 vectors is nearly optimal The second series of computations were carried out for a three dimensional channel flow past a circular cylinder with Rep 100 Figure 6 6 shows snapshots of the computation that was motivated by a suite of benchmark problems used as round robin tests to evaluate incompressible flow solution methods devel oped under the Deutsche Forschungsgemeinschaft DFG Priority Research Program Flow Simulation on High Performance Computers Figure 6 6a shows a snapshot of the 3 D z vorticity field Isosurfaces of the corresponding pressure field are shown in Figure 6 6b In this computation a parabolic velocity profile was prescribed at the inlet to the channel with no slip
73. derived from large eddy simulations Journal of Fluid Mechanics 200 511 562 1989 V K Chakravarthy and S Menon Large eddy simulations of turbulent premixed flames in the flamelet region to appear Combustion Sience and Technology 2001 W W Kim A New Dynamic Subgrid Scale Model for Large Eddy Simulation of Turbulent Flows PhD thesis Georgia Institute of Technology Atlanta Georgia September 1996 C C Nelson Simulations of Spatially Evolving Compressible Turbulence Using a Local Dynamic Subgrid Model PhD thesis Georgia Institute of Technology Atlanta Georgia December 1997 Incompressible Flow Solver BIBLIOGRAPHY 95 M Germano Ugo Piomelli Parviz Moin and William H Cabot A dynamic subgrid scale eddy viscosity model Physics of Fluids A 3 7 1760 1765 July 1991 96 D K Lilly A proposed modification of the germano subgrid scale closure method Physics of Fluids A 4 3 633 635 March 1992 97 S Liu C Meneveau and J Katz On the properties of similarity subgrid scale models as deduced from measurements in a turbulent jet Journal of Fluid Mechanics 275 83 1994 98 G K Batchelor The Theory of Homogeneous Turbulence Cambridge University Press 1953 99 G De Vahl Davis Natural convection in a square cavity a comparison exercise International Journal for Numerical Methods in Fluids 3 227 248 1983 100 G De Vahl Davis Natural convection of air in a square cavity a bench mark numeri
74. du uj y 5 Seed 2 5 Ox Ox i 2 1 2 Mass Conservation The mass conservation principle in divergence form is SURE 277 MINIS 2 6 ot Ox j In the incompressible limit the velocity field is solenoidal Qui 0 2 7 TEA 2 7 which implies a mass density transport equation dp O 2 8 ar ax ce For constant density Eq 2 8 is neglected with Eq 2 7 remaining as a constraint on the velocity field LS DYNA provides the capability to transport up to 10 species with Z Zo Zi9 representing the 10 species mass concentrations In order to simplify the presentation a single mass fraction is presented representing a binary mixture Mass conservation applied to one species yields for Zi OZ 0Z wi Poe Pas d m 2 9 where J is the diffusional mass flux rate and m is a volumetric mass source The mass source may include the injection of mass concentration from a boundary or the source sink terms from chemical reac tions The diffusional mass flux rate is based on Fick s law of diffusion OZ fh pps 2 10 J where 7D is a tensorial mass diffusivity Typically mass diffusivities are only available as scalars so that 0Z Ji p2D 3s 2 11 4 Incompressible Flow Solver 2 1 CONSERVATION EQUATIONS In the most general form the species concentration transport equations are OZ OZ MN Puy en CU where indicates the mass concentration i e I 1 2
75. e 01 1 441e 00 2 038e 00 2 635e 00 3 232e 00 3 829e 00 4 426e 00 5 023e 00 Figure 10 24 Pressure snapshot at t 0 7758 s Incompressible Flow Solver 129 CHAPTER 10 EXAMPLE PROBLEMS 130 Incompressible Flow Solver 10 4 THE SHEATH FLOW CHAMBER 10 4 The Sheath Flow Chamber Recent advancements in micro machining has permitted the development of portable particle analyzers that sort and categorize particles e g blood cells at rates exceeding 5000 particles per second At the core of these devices is a fluid dynamic focusing chamber commonly referred to as a sheath flow chamber The sheath geometry effects the flow patterns which in turn effects the particle trajectory through the sheath and into the capillary tube used to count and categorize the particles In the ideal situation the sheath flow is steady laminar and void of recirculation zones and swirl effects This example presents the analysis of the flow field in a proposed sheath chamber manufactured from etched silicon wafers The geometry for the sheath chamber is shown in Figure 10 25 In this configuration the inlet and outlet axial are square sections with an edge length of 100 um and a hydraulic diameter of 100 um The overall length of the device is 1250 um with a 100 um thickness Axial velocities ranging from 1 m s to 10 m s are typical when water is used as the sheath fluid In this analysis the properties of water at room temperature are density p
76. e Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 5 BOUNDARY PRESSURE CFD SET Purpose Apply a pressure load over each segment in a segment set The pressure convention follows Figure 9 1 This keyword activates what is commonly referred to as the pressure boundary condition according to 5 3 1 Card Format fo EEA ESE EE CC LIII Remarks VARIABLE DESCRIPTION SSID Segment set ID see SET SEGMENT LCID Load curve ID see DEFINE CURVE P Pressure to be applied Remarks 1 The load curve multipliers may be used to increase or decrease the pressure The time value is not scaled Incompressible Flow Solver 79 CHAPTER 9 KEYWORD INPUT a t t t s n4 n3 n y E ny nj No r n5 b Figure 9 1 Node numbering for pressure boundary segments a two dimensional elements b three dimensional elements Positive pressure acts in the negative t direction For two dimensional problems repeat the second node for the third and fourth nodes in the segment definitions 80 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 4 CONTROL_CFD_AUTO Purpose Set the time step control options for the Navier Stokes flow solver CON TROL_CFD_GENERAL is used in conjunction with with keyword to control the flow solver time integration options Card Format al EES SESE emm I Remarks VARIABLE DESCRIPTION IAUTO Set the time step control type EQ 0 IAUTO 1 for fixed time step size EQ 1 Fixed tim
77. e additional output variables associated with level 2 level 3 In addition to the mean and derived flux quantities higher order averages are generated that include velocity skewness and velocity flatness Table 8 6 shows the additional output variables associated with level 3 At this time the scalar dissipation is not included in the time average database Incompressible Flow Solver 73 CHAPTER 8 FLOW STATISTICS Variable Name LS POST Label Reynolds Stresses curas lt ul u2 gt lt ul u3 gt u2 u2 u2 u3 RA a gt Velocity Pressure Correlations Velocity Temperature Correlations Velocity Species 121 lt ul Z1 gt Correlations lt u2 Z1 gt lt u3 Z1 gt lt ul Z2 gt lt u2 Z2 gt lt u3 Z2 gt lt ul Z10 gt lt u2 Z10 gt lt u3 Z10 gt Table 8 5 LS POST level 2 derived variables and LS POST Labels Variable Name LS POST Label 2 Velocity Skewness lt u gt lt u gt lt X skewness gt uj gt u 23 lt Y skewness gt ui gt u lt Z skewness gt Velocity Flatness j lt X flatness gt Y flatness Z flatness Table 8 6 LS POST level 3 derived variables and LS POST labels 74 Incompressible Flow Solver Chapter 9 Keyword Input This chapter outlines the primary keywords used to control the incompressible flow solver and discusses how these keywords interact with LS DYNA s keywords for solid struct
78. e incompressible navier stokes equations Part 1 International Journal for Numerical Methods in Fluids 1 17 43 1981 55 Robert L Sani P M Gresho Robert L Lee D F Griffiths and Michael Engelman The cause and cure of the spurious pressures generated by certain fem solutions of the incompressible navier stokes equations Part 2 International Journal for Numerical Methods in Fluids 1 171 204 1981 56 David Griffiths and David Silvester Unstable modes of the q1 pO element Technical Report NA 257 University of Manchester UMIST Manchester Centre for Computational Mathematics September 1994 57 T J R Hughes and L P Franca A mew finite element formulation for computational fluid dynam ics Vii the stokes problem with various well posed boundary conditions symmetric formulations that converge for all velocity pressure spaces Computer Methods in Applied Mechanics and Engi neering 65 85 96 1987 58 David J Silvester and N Kechkar Stabilised bilinear constant velocity pressure finite elements for the conjugate gradient solution of the stokes problem Computer Methods in Applied Mechanics and Engineering 79 71 86 1990 59 N Kechkar and D Silvester Analysis of locally stabilised mixed finite element methods for the stokes problem Mathematics of Computation 58 1 10 1992 60 David Silvester Stabilised mixed finite element methods draft preprint January 1995 61 Sean Norburn and David
79. e step based on DTINIT default see CONTROL CFD GENERAL EQ 2 Time step based on CFL stability for INSOL 3 1 EQ 3 Automatic time step selection based on local truncation error using a second order Adams Bashforth predictor and trapezoidal rule corrector EPSDT Set the tolerance for the local truncation error in time EQ 0 EPSDT 1 0e 3 default DTSF Set the maximum time step scale factor that may be applied at any given time step This sets the upper limit on the amount that the time step can be increased during any given step EQ 0 DTSF 1 25 default DTMAX Set the upper limit on the time step size This value puts a ceiling on how far the time step may be increased for IAUTO 3 EQ 0 10 0 DTINT default Incompressible Flow Solver 8l CHAPTER 9 KEYWORD INPUT Remarks 1 There are multiple solver options for a variety of flow related physics in LS DYNA The selection of the time step control mechanism is dependent upon the flow solver that is selected IAUTO 1 may be used with any of the solution methods AUTO 2 forces the time step to be based on either stability or the CFL number see CONTROL_CFD_GENERAL for either the explicit INSOL 1 or the semi implicit INSOL 3 methods For IAUTO 2 the ICKDT parameter may be used to control the interval at which the time step is checked and adjusted The use of the second order predictor corrector time step control is restricted to the fully implicit solver 1 e INSOL 3 and
80. econditioned iterative techniques making it possible to tune the algorithm to the computing platform i e parallel vector or super scalar The second order projection algorithm can accurately track shed vortices and is amenable to the incorporation of either simple or complex multi equation turbulence sub models appropriate for a broad spectrum of applications CHAPTER 1 INTRODUCTION 1 1 Guide to the User s Manual The purpose for this document is to provide sufficient information for an experienced analyst to use LS DYNA s flow solver in an effective way The assumption is that the user is somewhat familiar with modern computing practices common CFD practices and to a certain degree the current CFD literature This manual provides sufficient references to the literature to permit the interested reader to pursue the technical details of LS DYNA In Chapters 2 6 an overview of the theoretical background for LS DYNA is presented Chapters 7 and 8 provide an overview of the turbulence modeling and the computation of derived flow statistics Chapter 9 presents information on LS DYNA s keyword input Several sample calculations are presented in Chapter 10 to provide guidance for new users 2 Incompressible Flow Solver Chapter 2 Governing Equations This chapter presents the basic forms of the partial differential equations that LS DYNA s flow solver treats and Chapters 3 4 provide a general description of the methodologies emp
81. ee dimensions For the local stabilization the entries of the stabilization matrix are calculated according to Eq 6 5 and 6 6 but only for the faces shared with elements in the same macro element as shown in Figure 6 1b Since the same scaling is used for both the local and global jump stabilization the only remaining pa rameter to be determined is B From the inherited scaling of the stabilization matrix B 0 provides the limit where no stabilization is applied For D 1 the stabilization matrix will have entries of the same order as the original PPE In the context of Stokes flow Norburn and Silvester 61 bounded B by com puting the extremal eigenvalues for the stabilized Schur complement as a function of B They found that 0 01 lt B lt 0 1 minimizes the distance between the extremal eigenvalues Numerical testing has shown that in general values of p in this range yield acceptable results for the stabilized PPE as well Figure 6 2 shows the variation of the iteration count for SSOR preconditioned conjugate gradient as a function of B for a 2 D vortex shedding problem and 3 D flow past a circular cylinder using both local and global jump stabilization This plot shows that for B gt 0 002 in 2 D B gt 0 004 in 3 D there is very Incompressible Flow Solver 43 CHAPTER 6 PRESSURE SOLUTION METHODS 450 D Unstabilized is E ni i 3 D Global 3 D Local 400 350 300 250 No of Iterations p Fi
82. ence statistics such as the mean and fluctuating velocity Reynolds stresses turbulent kinetic energy velocity temperature and velocity pressure correlations and higher order statistics such as velocity skewness and flatness may be generated LS DYNA generates a time averaged database the d3mean file that may be used with LS POST which provides the capability to select a time averaging window for the generation of all the derived turbulence statistics In this model LS DYNA generates the raw statistical data that is analyzed by LS POST to generate flow statistics In 8 1 a theoretical overview is presented as background while 8 2 summarizes the statistics level generated by LS DYNA and 8 3 presents the derived statistics available in LS POST Note that LS POST support for d3mean and d3thins files will only be available in LS POST 2 0 and later versions d3thins is the new time history database for the incompressible flow solver 8 1 Derived Flow Statistics Some of the most important statistical quantities of interest and the most physically relevant are the mean values of the flow variables sometimes referred to as moments The first two moments are the mean and variance These moments are routinely measured in experiments and form the basis of the Reynolds Averaged Navier Stokes RANS solution approach for turbulent flows There are several techniques used to obtain these moments from raw data Before proceeding with an overvi
83. endmas 20 00 2000 9 5 9 95 599 55559525255925529955559525259555552555555952525255555555525595555559555 Define the load curves DEFINE CURVE lcid sidr sfa sfo offa offo dattyp 1 al o1 0 0 1 0 1000000 0 1 0 DEFINE CURVE Incompressible Flow Solver 113 CHAPTER 10 EXAMPLE PROBLEMS lcid sidr 2 al 0 0 1000000 0 sfa sfo O oOo OF 0 0 offa offo dattyp 99933999959999999999 Sooo SSS SS Sooo SSS SS SSS SS SS SSS SS SSS SSSSEPUP OT USUS SS 5 Setup the boundary conditions Setup the no slip BC s set all velocity components at once BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 202 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 2 202 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 3 202 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 4 202 2 1 0000000 Temperature BC s BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 2 104 1 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 3 104 2 1 0000000 999339995989999999999 99 99 ooo SoS SS Sooo SSS S SS SSS SS SS SSS SS SSS SS SSS SS SSS USUS SS 5 Setup non dimensional Boussinesq fluid properties 5 MAT CFD CONSTANT mid ro T 1 0000 gx gy 0 0 1 0 diff6 diff7 5 Section amp part 5 SECTION SHELL sid elform 1 31 t1 t2 114 mu ST AD gz diff8 shrf t3 k 1 0 diffl diff9 nip t4 C
84. eraging used to ob tain the mean flow fields depends on the flow regime and the type of turbulence model that is applied The three types of averaging typically used in turbulence modeling are time space and ensemble averaging Time averaging is used for flows that are considered stationary i e the mean values are independent of the starting time for the averaging process Stationary flows are also referred to as statistically steady i e the mean values are not a function of time For a generic flow variable x t where 6 is taken to represent flow variables such as velocity pressure etc the time average is 1 to T Jaj i ai T T t x t dt 7 1 where x is the averaged quantity and T is the duration for the averaging procedure Note that in subsequent chapters the average is also identified as lt Q x gt or lt gt Using the average flow variable the instantaneous value may be decomposed into average b and fluctuating components as olit Pi 0 xit 7 2 For flows that have a slowly varying mean that is not turbulent in nature e g flows with an imposed mean periodic behavior the time averaging procedure can be extended to account for the long and short time scales In this situation the time average is computed using Eq 7 1 where Ti lt T lt To i e the time scale for the averaging must be large relative to the turbulent fluctuations and small relative to the time scale associated with
85. ergence tolerance i e V ul rws This tolerance is used for the initial startup procedure to insure that proper initial conditions are prescribed for the momentum equations EQ 0 DIVU 1 0e 5 default Time weighting for viscous diffusion terms Valid values are 0 lt Ox 1 with Ox 1 2 for second order accuracy in time EQ 0 THETAK O 5 default Time weighting for advection terms Time weighting for body forces and boundary conditions Valid values are 0 lt Op lt 1 with Of 1 2 for with for second order accuracy in time EQ 0 THETAF 0 5 default Set the equation solver type for the momentum equations EQ 0 MSOL 20 default EQ 20 Jacobi preconditioned conjugate gradient method EQ 30 Jacobi preconditioned conjugate gradient squared method default when IADVEC 40 Set the maximum number of iterations for the iterative equation solver EQ 0 MAXIT 100 default Set the interval to check the convergence criteria for the iterative equation solver EQ 0 ICHKIT 2 default Activate the output of diagnostic information from the equation solver EQ 0 Diagnostic information is off default EQ 1 Diagnostic information is on Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS VARIABLE DESCRIPTION IHIST Activate the generation of a convergence history file from the equation solver The ASCII history files are velx his vely his and velz his EQ 0 Convergence his
86. eria for the pressure solver EQ 0 EPS 1 0e 5 default Set the number of A conjugate vectors to use during the iterative pressure solve EQ 0 NVEC 5 default LT 0 A conjugate projection is disabled Set the stabilization type EQ 0 ISTAB 2 default EQ 1 Local jump stabilization EQ 2 Global jump stabilization EQ 1 No stabilization is active Stabilization parameter for ISTAB 1 2 Valid values for the stabilization parameter are 0 lt p lt 1 EQ 0 BETA 0 05 default Solid element set or shell element set see SET SOLID SET SHELL OPTION to be used for the prescription of hydrostatic pressure Set the hydrostatic pressure level This value multiplies the values of the load curve specified with the LCID option EQ 0 PLEV 0 0 default Load curve to be used for setting the hydrostatic pressure By default LCID 0 which forces a constant pressure level to be set at the level prescribed by PLEV EQ 0 LCID 0 default Incompressible Flow Solver 89 CHAPTER 9 KEYWORD INPUT 9 1 8 CONTROL CFD TRANSPORT Purpose Activate the calculation of transport variables and associated solver parameters to be used for the auxiliary scalar transport equations Card 1 is used to activate the auxiliary transport equations and Card 2 is used to set the mass matrix advection and time weighting options Card 3 is used to set the linear solver options such as the maximum iteration count and interval to check the convergen
87. ering a relatively cold quiescent fluid and a classical shear instability phenomena known as the Kelvin Helmholtz instabil ity The parameters for the momentum driven jet are prescribed so that the Reynolds number is Re 3561 based on a 15 mm slot width and Fr 326 Fr v gBATL where g 9 81 m s BAT 1 Le 15 mm and v 4 m s Here the relatively large Froude number indicates that the influence of buoyancy forces is small compared to the inertial forces i e a momentum driven jet The grid for the momentum driven jet is shown in Figure 10 17 where a single plane of symmetry along x Q is used with the jet half width H 2 In this computation an energy equation is solved in conjunction with the Navier Stokes equations using a Boussinesq fluid i e air The initial conditions consist of an initial div free velocity field with a free field temperature of 300K and an inlet air jet temperature of 400K The working fluid is air with a near unit Prandtl number Pr 0 71 resulting in a Peclet number Pe 2528 where Pe RePr The keyword file momjet p2 k is shown below In this problem a fixed CFLz1 was prescribed with the automatic time step selection based on the CFL number Due to the effect of the buoyancy terms the vertical jet velocity increased beyond the prescribed 4 m s so the interval to check and adjust the time is set to ickdt 2 A time accurate simulation is of interest here so the default second order time weighting para
88. esis to relate the deviatoric part of the SGS stress tensor to the strain rate via an eddy viscosity The basis for the Smagorinsky model lies in the observation that the turbulent eddy viscosity can be parameterized in terms of a length and a velocity scale This observation permits the eddy viscosity to be written in terms of a length scale and the rate of dissipation as ty CAB 3 7 24 where the length scale has been chosen as the filter scale A Assuming that the dissipation rate is in equilibrium with the rate of production of subgrid kinetic en ergy P Ti i j permits estimation of the eddy viscosity in terms of the filtered strain rate In the Smagorinsky model the turbulent eddy viscosity is ui p CsA 4 SiS 7 25 The Smagorinsky constant Cs is not universal and requires adjustment for each unique set of flow con ditions Various applications have shown that Cs typically is in the range of 0 1 lt Cs lt 0 24 see Rogallo and Moin 86 In LS DYNA the default for the baseline Smagorinsky model is Cs 0 1 see the CON TROL_CFD_TURBULENCE keyword in Chapter 9 7 4 Subgrid Kinetic Energy Models The Smagorinsky model has found broad application and is attractive from a computational perspective due to its simplicity While the Smagorinsky model provides approximately the correct amount of dissi pation to the resolved velocity field the model has several deficiencies For example the Smagorinsky constant is no
89. ethod yields the overall best performance relative to the re solve cost of a direct solver The effect of 50 Incompressible Flow Solver 6 4 PPE SOLVER PERFORMANCE Table 6 1 Effect of the A Conjugate projection on the PPE solution time for the 2 D momentum driven jet The grind times for the explicit algorithm have been normalized with respect to the solution time for 1000 time steps using a direct solver 71 72 Nrr is the average number of iterations required per time step to solve the PPE problem Nel 11250 Nnp 11466 N 146 direct solver grind time 1 2254 x 107 CPU seconds per element per time step Explicit Algorithm No Stabilization JPCG SSOR PCG ESSOR PCG No Sr Grind Grind Grind vectors AUT ius 7 To TT we L Mr Bus a s RE Algorithm Global Jump Stabilization D 0 05 JPCG SSOR PCG ESSOR PCG_ No mum Grind Grind Grind vectors Nir Time 7 2 TIT nus 7 n MT Time 7 TE 92 Explicit Algorithm Local Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR CG No m Grind Grind Grind vectors AUT m 7 T Mir me 7 a Mir me 7 a Incompressible Flow Solver 51 CHAPTER 6 PRESSURE SOLUTION METHODS Table 6 2 Effect of the A Conjugate projection on the PPE solution time for the 2 D momentum driven jet The grind times for the projection algorithm have been normalized with respect to the solution time for 1000 time steps using a direct solver 71 72 Nrr is the average number of iterations
90. ew of the methods used to derive statistics from a time accurate flow simulation several definitions are required 65 CHAPTER 8 FLOW STATISTICS 8 1 1 Reynolds Averaged Statistics Reynolds time averages homogeneous spatial averages and ensemble averages are three types of averag ing commonly used to determine mean quantities or moments of a turbulent flow Reynolds averaging describes the time invariant mean quantities of a stationary flow as a function of position Ensemble av erages are well suited for transient flows and describe the mean values as a function of both space and elapsed time Spatial averages describe the spatial invariance of mean values in a homogeneous flow that may or may not be stationary A stationary flow is one for which the mean values are independent of the starting time for the averaging process Stationary flows are also referred to as statistically steady i e the mean values are not a function of time For a generic flow variable o x t where u v w p etc the Reynolds time average is 1 to T lt o x gt lim x t dt 8 1 To T where lt x gt is the averaged quantity and T is the duration for the averaging process as shown in Figure 8 1 Remark Note that the time average was defined in an equivalent way in Chapter 7 and was referred to as o x Here lt gt is used because this notation is directly translated into the graphics database and LS POST In a discrete
91. f the A conjugate vectors in proceeds according to the algorithm above Experimentation with vectors representing a random solution hydrostatic and simple polynomials have been tested with limited success In the case of simple flows e g Pouiselle flow the linear polynomials contribute directly to the initial solution and reduce the number of CG iterates dramatically However in moderately complex domains the simple polynomials provide little reduction in computational cost It is thought that initial polynomial fields that respect the implied pressure boundary conditions may provide an effective set of seed vectors but this has not yet been implemented To illustrate the use of the A conjugate projection Figure 6 4a shows a snapshot of the pressure field for a Re 100 vortex shedding computation In addition Figures 6 4b 6 4f show snapshots of the vectors based on the previous five time steps It is clear that provides primarily long wavelength information while the other 4 vectors provide detailed information about the wake Thus the vectors 5 5 may be viewed as short wavelength corrections to that yield the best approximation to the current pressure field The A conjugate projection procedure in effect selects the appropriate information from each vector in order to minimize the residual in the A norm before performing any CG iterations The A conjugate projection procedure retains both long and short wavelength
92. f the formation and assembly of the right hand side in Eq 3 5 and the subsequent PPE solution in Eq 3 6 The right hand side vector is formed element by element using a right to left matrix vector multiply with reduced integration operators for the diffusive and advective terms In this context reduced integration is considered synonymous with single point Gaussian quadra ture Using single point integration the element area and gradient operators may be written strictly in terms of element local nodal coordinates 1 Af Pays x24y32 3 20 C x l X 54e 24 31 24 Y31 1 yo E day 21 e zae 42 13 x42 x13 3 21 In Eq 3 20 and 3 21 Xab Xa Xp where the subscripts a and b identify the local node number and range from 1 to 4 for the two dimensional bilinear element The fact that C 3 Cx1 and C4 C 2 permits the storage of only the unique values in the gradient operator at the element level In two dimensions this requires only 4 floating point values per element for Cy and C7 The computation of the element level gradient operators in three dimensions is somewhat more involved To begin the element local nodal coordinates in the referential domain are defined as follows o TL TT o oes Dust n 1 1 1 1 1 1 1 1 3 22 CS Slee Uo DN E Y The Jacobian evaluated at the element centroid or central Gauss point is defined in terms of the element local nodal coordinates Xe ye Ze
93. fusivity 7 At lij Hij py uiu 4 28 4 20 Fully Implicit Projection Method The fully implicit projection method relies upon the implicit treatment of the advection terms as well as the viscous terms in the momentum energy and species transport equations The fully implicit method is implemented with automatic time step selection based on local time truncation error requiring a slightly different start up procedure than the semi implicit projection method The details of the automatic time step selection may be found in Gresho et al 7 The start up procedure for the fully implicit algorithm follows the start up algorithm outlined in 4 1 to obtain u and p but requires the consistent calculation of the acceleration field at t 0 4 2 1 Start Up Given an initial velocity field u the start up procedure for the fully implicit algorithm proceeds as fol lows 1 Calculate the partial acceleration i e acceleration neglecting the pressure gradient by solving the mass matrix problem M F Ku A u u 4 29 where is the instantaneous acceleration neglecting the pressure gradient 2 Solve the global PPE for the current pressure field C7 M c p C73 gl 4 30 3 Calculate the total acceleration at t 0 CM7 p 4 31 32 Incompressible Flow Solver 4 2 FULLY IMPLICIT PROJECTION METHOD 4 2 2 First Time Step The first time step at t O proceeds as follows 1 Compute a
94. gure 6 2 Variation of iteration count vs stabilization parameter for a 2 D Re 100 vortex shedding mesh Nel 11200 and 3 D Re 100 flow past a cylinder Nel 22380 The iteration count for the un stabilized PPE is indicated by an arrow for the 2 D and 3 D problems little variation in the iteration count with respect to D In addition the benefit of stabilization appears to be considerably larger for 3 D a factor of 6 reduction in iteration count for 3 D compared to a factor of about 1 3 for the 2 D case The marked decrease in iteration count in 3 D is undoubtedly a consequence of the relatively larger null space for the PPE in three dimensions Operational experience has shown that 0 01 lt B lt 0 1 works well for most applications while B gt 0 25 can often yield inaccurate velocity fields The effect of pressure stabilization on the convergence rate is shown in Figure 6 3 which shows conver gence histories for coarse medium and fine grids for a vortex shedding problem Here a convergence criteria of e lt 10712 was used to obtain an initial RMS divergence of 1078 using global stabilization Similar results may be found in Gresho and Sani 5 where the local and global jump stabilization formu lations have been applied to a 3 D problem using an ICCG 0 preconditioner with similar improvements in convergence rate 6 3 The Projection CG Method The solution of the time dependent incompressible Navier Stokes requires the repeated
95. hat elements are uniquely assigned to processors while the nodes at the sub domain interfaces are stored redundantly in multiple processors Figure 3 3 illustrates a non overlapping domain decomposition for a simple two dimensional vortex shedding mesh partitioned for four processors Af ter obtaining the domain decomposition the assignment of nodes boundary conditions and materials to individual processors is performed internally before parallel execution begins Figure 3 3 Four processor spatial domain decomposition 3 2 2 Parallel Assembly via Message Passing The finite element assembly procedure is an integral part of any finite element code and consists of a global gather operation of nodal quantities an add operation and a subsequent scatter back to the global memory locations A complete description of the sequential assembly algorithm may be found in Hughes 8 The assembly procedure is used to both form global coefficient matrices right hand side vectors and matrix vector products in an element by element sense In the case of the explicit time integration algorithm the emphasis is upon the assembly of the element level contributions to the global right hand side vector AY f Keu Ae u uz 3 40 22 Incompressible Flow Solver 3 2 THE DOMAIN DECOMPOSITION MESSAGE PASSING PARADIGM 13 14 15 16 9 10 11 12 5 6 7 Sub domain PO Sub domain P1 8 Local
96. hat have been artificially imposed to emulate far field conditions in a large physical domain For the SET option define the following card Card Format CU RIESE EREAERESED For the SEGMENT option define the following card Card Format pos a fos fa fos fe for fa VARIABLE DESCRIPTION SSID Segment set ID N1 N2 Node ID s defining the segment 76 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 2 BOUNDARY PRESCRIBED CFD OPTION Available options include NODE or SET Purpose Define an imposed nodal variable velocity temperature species etc on a node or a set of nodes for the Navier Stokes flow solver Do not use the NODE option in r adaptive problems since the node ID s may change during the adaptive step The prescription of all nodal variables in the incompressible flow solver are defined using this keyword It is similar in function to the BOUNDARY PRESCRIBED MOTION OPTION keyword but permits the specification of boundary conditions for velocities mass concentration species temperature turbulent kinetic energy etc Additionally the keyword optionally provides for prescribing all velocities for exam ple on a no slip and no penetration surface Multiple instances of the keyword permit individual nodal variables to be prescribed using independent load curves and scale factors Card Format oo a fate fos fe for fa Incompressible Flow Solver 77 CHAPTER 9 KEYWORD INPUT VARIABLE
97. hat is referred to as global jump sta bilization and one that is referred to as local jump stabilization The global jump stabilization first proposed by Hughes and Franca 57 and the local jump stabilization techniques of Silvester and his co workers 58 59 60 61 are applied to the PPE problem for both the explicit and projection methods In effect the jump stabilization techniques provide an a priori filter for the weakly unstable pressure modes associated with the Q1Q0 element The stabilized Q1Q0 element yields a regularized saddle point problem for the projection method M C utl Mya bee pan e where S is a symmetric positive semi definite matrix Here pressure stabilization results in an approximate projection method since the associated stabilized PPE is no longer constructed using only the discrete div and grad operators The Schur complement of Eq 6 3 i e the stabilized pressure Poisson equation PPE is jc tc if s a cTan gt 6 4 where S remains to be defined for the global and local jump stabilization techniques The global jump stabilization formulation attempts to control the jump in pressure across element bound aries and results in a PPE that is perturbed by a pressure diffusion operator The off diagonal entries in the global jump stabilization matrix are defined as C M C sy p E Ch S yy i v at 6 5 where J and J identify adjacent elements that share a common face as shown in Figure
98. he computation of the derived turbulence statistics The starting time for the collection of data in LS DYNA is controlled by the TSTART parameter on the DATABASE BINARY D3MEAN keyword For each time window the LS DYNA collects the raw statistics for the active time window Considering a generic variable x t the time average for a given window is simply 1 n lavg x gt gntavg qn gt x t AU 8 22 j n where avg is the number of time samples accumulated in the window 70 Incompressible Flow Solver 8 2 LS DYNA STATISTICS LEVELS LS DYNA records the time associated with the end of each time window the elapsed time associated the window and the number of samples collected in the window The elapsed time pr tlavg _ 1 and the number of time samples during this time interval avg is written into the d3mean database The number of time samples in a window is controlled by the JAVG parameter on the DATABASE BINARY D3MEAN keyword 8 2 1 Level 1 Mean Statistics For level 1 time averaged statistics LS DYNA accumulates the variables shown in Table 8 1 with their associated database labels for the d3mean file In this statistics level derived mean quantities that are based on the user specified time window include the velocity lt u gt pressure lt p gt vorticity lt Qj gt and for 2D problems the vorticity lt gt and the stream function lt y gt For thermal problems the mean quantities a
99. hree dimensional finite element code for solid and structural mechanics user manual Technical Report UCRL MA 107254 Rev 1 Lawrence Livermore National Laboratory 1993 18 D R Kincaid T C Oppe and D M Young Itpackv 2d user s guide Technical report Center for Numerical Analysis The University of Texas at Austin May 1989 19 C Farhat and Michel Lesoinne Automatic partitioning of unstructured meshes for the parallel solution of problems in computational mechanics International Journal for Numerical Methods in Engineering 36 745 764 1993 20 A Pothen H Simon and K Liou Partitioning sparse matrices with eigenvectors of graphs SIAM J Matrix Anal Appl 11 430 452 1990 21 H D Simon Partitioning of unstructured problems for parallel processing Computing Systems in Engineering 2 135 1991 22 S T Barnard and H D Simon A fast multilevel implementation of recursive spectral bisection for partitioning unstructured problems In 6th SIAM Conference on Parallel Processing for Scientific Computing pages 711 718 1993 23 B Hendrickson and R Leland An improved spectral graph partitioning algorithm for mapping parallel computations Technical Report SAND92 1460 Sandia National Laboratories Report September 1992 24 B Hendrickson and R Leland Multidimensional spectral load balancing Technical Report SAND93 0074 Sandia National Laboratories Report January 1993 25 B Hendrickson and
100. ial inlet In addition to approaching a steady state the flow field is also symmetric about the x y plane as shown by the pressure isosurfaces in Figure 10 29 In contrast the y velocity field is skew symmetric about the x y plane as shown by the time history plot in Figure 10 28 The isosurfaces of helicity u shown in Figure 10 30 highlights the rotational structure of the velocity field at the outlet plane Here there are four counter rotating cells each pair originating at the side inlet port Incompressible Flow Solver 131 CHAPTER 10 EXAMPLE PROBLEMS Side Inlet Port Node Segment Set ID 3 Axial Inlet Port Node Segment Set ID 2 N Side Inlet Port Node Segment Set ID 1 Axial Outlet Port Figure 10 25 Sheath chamber mesh with 33088 elements and 37196 nodes Example sheath p2 k KEYWORD TITLE Prototype sheath flow chamber 8 98 99999999999999 99 9 PS SS Sooo ooo oS SS SO OSS SS SSS SS SS SSS SS SSS SS SSS SSSSSS Sssss Setup the CFD problem soln 4 CONTROL SOLUTION soln 4 CONTROL CFD GENERAL insol dtinit cfl ickdt 3 1 0e 6 4 0 2 CONTROL CFD MOMENTUM 132 istats tstart iavg Incompressible Flow Solver imass advec i 10 msol maxit 20 100 CONTROL CFD PRESSURE ipsol maxit 22 250 nvec stab 10 CONTROL TERMINATION endtim endcyc 1500 00000 2000 ifct HH ichkit 2 ickit 2 beta dtmin 10 4 THE SHEATH FLOW CHAMBER divu thetak 1 0e 7 i
101. ical element local node numbering scheme coordinate system and centroid velocity for the 2 D and 3 D elements Figure 3 2 Grid parameters for a Two dimensional element with centroid velocity and characteristic element dimensions hg and hy b Three dimensional element with centroid velocity and characteristic dimensions he hy and he The grid Re and CFL numbers are defined as u h Ross v 3 11 ju h Ar Incompressible Flow Solver 15 CHAPTER 3 EXPLICIT TIME INTEGRATION where i n Gare the element local coordinates The grid Re and CFL numbers rely upon the projection of the centroid velocity onto the element local coordinate directions that are oriented according to the canonical local node numbering scheme for each element type 11 12 In order to use Eq 3 12 to estimate a stable time step a unit vector for each element local coordinate direction is defined as n h where denotes the unit vector for each of the n coordinate directions Using the grid Re and the element size h the advective diffusive stability limit becomes 1 At vitis m 3 14 where a minimum over all elements and all element local coordinate directions establishes a global mini mum time step The advective stability limit is established in a similar manner using CFL hi u hj gt The stable time step is based upon the minimum time step derived from either Eq 3 14 or 3 1
102. ime scales in high Reynolds numbers flows i e direct numerical simulation remains out of reach even for the largest supercomputers available today particularly when complex geometry or complex physics is involved For industrial flows which are typically turbulent with relatively high Reynolds numbers turbulence modeling is the pacing technology for computational fluid dynamics Turbulence models range in complexity from simple mixing length or zero equation models to n equation models that provide second moment closure The n equation models require the solution of additional partial differential equations that describe the transport of turbulence quantities e g turbulent kinetic energy The n equation models range in complexity from one equation models e g the Spalart Allmaras model 73 and the Baldwin Barth model 74 to models that provide second moment closure and require up to seven additional transport equations Two equation models are typified by the standard k model and its multitude of variations Large eddy simulation is a promising alternative to traditional turbulence modeling approaches that rely upon Reynolds averaging and closure models based upon ad hoc transport equations for the rate of dis sipation Unlike traditional Reynolds averaged models LES provides a high degree of accuracy with minimal empiricism However the computational cost associated with performing pure LES i e where a significant fraction
103. in Figure 10 4 In addition a time history of the global kinetic energy 1 2u7 Mu is shown in Figure 10 5 As shown by these plots a steady state solution is achieved after approximately 100 time units A plot of the velocity profile and the corresponding linear pressure field are presented in Figures 10 6 and 10 7 In Figure 10 6 velocity vectors are plotted at the outlet boundary to show the parabolic profile The outlet velocity profile is compared to the exact solution to the Poiseuille problem in Figure 10 8 where it is seen that the computed velocities interpolate the exact solution at the nodal points Incompressible Flow Solver 105 CHAPTER 10 EXAMPLE PROBLEMS grid parameters CFL based time step Local coordinate direction lesse Element n mber sj v4 ue REIHE El ment lengthlia iia ti Oats oe EA OS Grid Reynolds noia ia ERU SEE ES GEO CEBOnumbes cec OMe Moo Coed Ade E ecd CFL based time step llle eee Maximum xi grid Reynolds numbers Element NO asc caer e AA Rae ea Element length RR ee nente Grid Reynolds number lvl ae eye ans Maximum eta grid Reynolds numbers Elemente nOu su dtu de va eoa Cut tog mele ou Sea ya Ns Element lengths dos sb eae dae ote A ae dorus Grid Reynolds number 0 020002 ee eee Maximum xi grid CFL number Element HO sie ghee ave ae eae NE ERE Y aye e lee ROS BRR Element Length reke oia tete e Pu dr ers Grid CTL UMD d ee OER Maximum
104. ing 2 1 1 4 1981 69 J M Ortega Efficient implementations of certain iterative methods SIAM Journal for Scientific and Statistical Computing 9 5 882 891 September 1988 70 S C Eisenstat J M Ortega and C T Vaughan Efficient implementation of a class of pre conditioned conjugate gradient methods SIAM Journal for Scientific and Statistical Computing 11 5 859 872 1990 71 O O Storaasli D T Nguyen and T K Agarwal A parallel vector algorithm for rapid structural analysis on high performance computers Technical Report 102614 NASA April 1990 72 T K Agarwal O O Storaasli and D T Nguyen A parallel vector algorithm for rapid structural analysis on high performance computers AZAA 90 1149 1990 73 P R Spalart and S R Allmaras A one equation turbulence model for aerodynamic flows In AIAA 92 0439 Reno Nevada January 1992 AIAA 30th Aerospace Science Meeting and Exhibit 74 Barrett S Baldwin and Timothy J Barth A one equation turbulence transport model for high reynolds number wall bounded flows In A AA 91 0610 Reno Nevada January 1991 AIAA 29th Aerospace Science Meeting and Exhibit 75 W C Reynolds The potential and limitations of direct and large eddy simulations In J L Lumley editor Whither Turbulence Turbulence at the Crossroads pages 313 343 Berlin March 1989 Springer Verlag 76 David C Wilcox Turbulence Modeling for CFD DCW Industries Inc La Canada Ca
105. ions there are two primary messaging steps The first occurs during the distributed assembly of the right hand side in Eq 3 41 The second occurs during the assembly of the pressure gradient in Eq 3 45 Thus there are 2 vector valued messaging steps For the solution Incompressible Flow Solver 25 CHAPTER 3 EXPLICIT TIME INTEGRATION of the PPE there is 1 vector valued messaging step per iteration where rr iterations are required in the conjugate gradient solution From this it is possible to estimate the communication cost per time step for the explicit algorithm on a per processor basis as Cstep 8 2 Nir NpoFNrtsend 3 48 26 Incompressible Flow Solver Chapter 4 The Projection Method The solution of the time dependent incompressible Navier Stokes equations poses several algorithmic is sues due to the div free constraint and the concomitant spatial and temporal resolution required to perform time accurate solutions particularly where complex geometry is involved Although fully coupled solution strategies are available the cost of such methods is generally considered prohibitive for time dependent simulations where high resolution grids are required The application of projection methods provides a computationally efficient alternative to fully coupled solution methods A detailed review of projection methods is beyond the scope of this chapter but a partial list of relevant work is provided Projection methods also commo
106. it projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix part 1 Theory International Journal for Numerical Methods in Fluids 11 587 620 1990 Philip M Gresho and Stevens T Chan On the theory of semi implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix part 2 Implementation International Journal for Numerical Methods in Fluids 11 621 659 1990 P M Gresho and R L Sani Incompressible flow and the finite element method Advection diffusion and isothermal laminar flow John Wiley amp Sons Chicester England 1998 Philip M Gresho and Robert L Sani On pressure boundary conditions for the incompressible navier stokes equations International Journal for Numerical Methods in Fluids 7 1111 1145 1987 Philip M Gresho Robert L Lee and Robert Sani On the time dependent solution of the Navier Stokes equations in two and three dimensions volume 1 page 27 Pineridge Press Swansea U K 1980 Thomas J R Hughes The Finite Element Method Prentice Hall Inc Englewood Cliffs New Jersey 1987 9 O C Zienkiewicz and Robert L Taylor The Finite Element Method volume 2 McGraw Hill Book Company Limited Berkshire England fourth edition 1991 10 Philip M Gresho and Robert L Sani Computational Fluid
107. ithms available in LS DYNA that may be applied to this problem 1 an explicit time integration method with a reduced integration element 2 a semi implicit projection method with explicit advection and 3 the fully implicit projection method with a smart time integrator In this example a the semi implicit projection method is applied with the selection of a backward Euler treatment of the viscous terms and a row sum lumped mass matrix The LS DYNA keyword file for this problem is shown below in a typewriter font The keyword in put is segregated into 5 primary parts 1 the control keywords 2 load curves 3 boundary condi tions 4 fluid properties and 5 I O controls In this example the second order projection algorithm insol 3 on the CONTROL CFD GENERAL keyword is used with the viscous terms treated in a fully implicit mode a lumped mass matrix and explicit advection terms This is obtained by setting thetak 1 0 imass 1 and iadvec 10 on the CONTROL CFD MOMENTUM keyword The CON TROL CFD AUTO keyword is used to enable automatic time step selection based on the maximum CFL number specified on the CONTROL CFD GENERAL keyword i e c 1 2 0 In addition the initial time step is set to dtinit 0 01 A direct solver ipsol 11 via the CONTROL CFD PRESSURE key word is used with a termination of endcyc 2000 time steps Constant fluid properties are set using the MAT CFD CONSTANT keyword which in the non dimensional context re
108. itions for the incompressible flow solver in LS DYNA Key words that are relevant to the boundary conditions are summarized in Chapter 9 and fully documented in the LS DYNA Keyword User s Manual 27 5 1 Node and Segment Sets Sets provide a generalized concept for grouping nodes elements and materials In LS DYNA sets provide the basic entity for the prescription of boundary conditions and are used for all keyword input Node sets consist of an arbitrary list of nodes that are treated as one entity for the application of nodal boundary conditions Segment sets consist of an arbitrary list of quadrilateral edges in two dimensions and a list of hexahedral faces in three dimensions Figure 5 1 shows a node set associated with inflow conditions and a segment set associated with outflow boundary conditions for an external flow problem In addition to node and segment sets the incompressible flow solver uses shell and solid sets to identify elements with specific attributes Additional information on sets and the SET keywords may be found in the LS DYNA keyword manual 27 5 2 Nodal Boundary Conditions The prescription of nodal boundary conditions for the incompressible flow solver encompasses all nodal degrees of freedom e g velocity temperature mass concentration species and turbulent kinetic en ergy These boundary conditions are typically referred to as essential or Dirichlet boundary condi tions and fix the nodal values of the field v
109. ive mass transfer coefficient as Ji ni ha Zi Z1 5 13 where ho is the mass transfer coefficient and Z is a reference species concentration Incompressible Flow Solver 39 CHAPTER 5 BOUNDARY CONDITIONS AND SOURCE TERMS 5 7 Pressure Levels In incompressible flow the pressure is typically only known up to an additive constant the hydrostatic pressure level The prescription of a specific hydrostatic pressure level is achieved through the CON TROL CFD PRESSURE keyword For the Q1Q0 element technology setting the hydrostatic pressure level requires a shell set in 2 D or a solid set in 3 D to identify those elements for which the pressure level will be fixed 40 Incompressible Flow Solver Chapter 6 Pressure Solution Methods This chapter presents a survey of solution strategies stabilization techniques and linear algebra procedures used in LS DYNA for calculating the pressure in the time dependent Navier Stokes equations The pri mary focus is on the treatment of the pressure Poisson equation deriving from the index 1 formulations of the Navier Stokes equations presented in Chapters 3 4 Operational experience has shown that de pending on the solution strategy the pressure solution can consume 80 90 of the CPU time per time step in semi implicit and fully implicit projection methods even more when explicit time integration is used Based on extensive operational experience with a variety of solution s
110. l ickdt istats tstart lavg 3 1 0e 4 1 0 2 CONTROL CFD MOMENTUM imass iadvec ifct divu thetak thetaa thetaf 3 10 1 1 0e 7 msol maxit ichkit iwrt ihist eps ihg ehg 20 100 2 0 0 CONTROL CFD PRESSURE ipsol maxit ickit iwrt ihist eps 22 300 2 0 0 1 0e 5 nvec stab beta 5 CONTROL CFD TRANSPORT itemp iden nspec ifrac 1 0 0 0 imass ladvec ifct thetak thetaa thetaf 3 10 1 itsol maxit ichkit iwrt ihist eps ihg ehg 0 CONTROL TERMINATION endtim endcyc dtmin endeng endmas 5 00000 1500 O O O 99999 99 99999 Define the load curves DEFINE CURVE 1 0 0 1 0 9999 0 1 0 DEFINE CURVE 2 0 0 0 0 9999 0 0 0 DEFINE CURVE 3 0 0 400 0 9999 0 400 0 CEEEEPEEREEREEEEEEEEEEEREEEEEEEEREEEEEEREREEREREEEPEEEEEREREEEEEEEEEEEEEEEEEEREERERS 5 Setup the boundary conditions 5 no slip no penetration BC s on nodeset 20 BOUNDARY PRESCRIBED MOTION SET nsid dof lcid sf 10 101 2 0 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 122 Incompressible Flow Solver 10 3 THE MOMENTUM DRIVEN JET 20 102 2 0 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 20 101 2 1 0000000 5 no x velocity on symmetry plane using nodeset 60 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 60 101 2 1 0000000 5 jet inlet velocity on nodeset 20 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 10 102 1 4 0000000 jet inlet temperature BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 10 104 3 1 000000
111. lass stabilization operator is Ci Ii PENIS C r H ej u F1 T2 Ts Pa 2 C rs 3 38 C4 Ta where f 1 0 C1 C2 Cs C4 AVV h YVmax VVmin 2 For the reduced integration element y stabilization see Belytschko et al 14 15 16 has also been investi gated y stabilization refers to the y vectors constructed from the hourglass modes for stabilization While Y stabilization is perhaps more robust than h stabilization this type of hourglass control also requires more operations and storage It is the author s experience that it is relatively more difficult to excite the hourglass modes in an Eulerian computation than in a Lagrangian computation e g a DYNA3D 17 sim ulation However y stabilization still requires fewer operations and less storage than the fully integrated two dimensional bilinear element In three dimensions this is not the case Table 3 1 shows the memory requirements and operations counts for a matrix vector multiply Ku for a variety of element formulations In 2 D y stabilization requires nearly the same storage as the fully integrated element stored in either a compact symmetric element form or in a global row compressed form However this element requires 9 more operations to achieve the matrix vector multiply when compared to the element by element matrix vector multiply with 2x2 quadrature In 3 D y stabilization is about 3 times more expensive to perform than the corresponding global
112. lication of the averaging procedure making it impossible to generate a closed system In order to compute the time averaged velocity and pressure a prescription for the Reynolds stress is required Ultimately the purpose for turbulence models is to provide the prescription of the Reynolds stress and close the system of time averaged equations albeit with varying degrees of complexity and SUCCESS 60 Incompressible Flow Solver 7 3 THE SPATIALLY FILTERED EQUATIONS 7 3 The Spatially Filtered Equations In this section the spatially filtered form of the conservation equations are presented Here an over bar indicates a spatially filtered quantity Applying the spatial filtering procedure Eq 7 7 to Eq 7 11 and 7 12 with the requirement that the filtering operation commutes with differentiation the resulting space filtered momentum equations are Ou Ou u op d T ches F pa oe a hg Pu tU UR M where the strain rate based on the resolved velocity field is 1 Ou S 22l J 7 18 d 2 2 The divergence constraint is gU 7 19 Ox Application of the spatial filter to the nonlinear advection term in the momentum equations produces a new term referred to as the subgrid scale SGS stress tensor t Using Germano s generalized moment ij definitions 84 the SGS stress tensor is Vif p ur uj ne ij and represents the effect of the subgrid scale variations on the resolved scale Like the
113. lifornia 1993 77 Chin Hoh Moeng A large eddy simulation model for the study of planetary boundary layer turbu lence Journal of the Atmospheric Sciences 41 13 2052 2062 July 1984 78 Peter R Voke and Socrates G Potamitis Numerical simulation of a low reynolds number turbulent wake behind a flat plate International Journal for Numerical Methods in Fluids 19 377 393 1995 Incompressible Flow Solver 145 BIBLIOGRAPHY 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 146 D Veynante and T Poinsot Large eddy simulation of combustion instabilities in turbulent pre mixed burners Technical report Center for Turbulence Research Stanford University Palo Alto California 1997 W W Kim and S Menon Large eddy simulations of reacting flow in a dump combustor Technical Report AIAA 98 2432 29th AIAA Fluid Dynamics Conference Albuquerque New Mexico June 1998 C Weber F Ducros and A Corjon Large eddy simulation of complex turbulent flows Technical Report AIAA 98 2651 29th AIAA Fluid Dynamics Conference Albuquerque New Mexico June 1998 F F Grinstein and C Fureby Monotonically integrated large eddy simulation of free shear flows In 36th Aerospace Sciences Meeting and Exhibit Reno Nevada January 1998 American Institute of Aeronautics and Astronautics D Kwak and W C Reynoldsand J H Ferziger Three dimensional time depende
114. loyed in their solution The interested reader may pursue the references included in these chapter for details on the incompressible flow solution algorithms used in LS DYNA and their implementation 2 1 Conservation Equations In the ensuing discussion indicial notation is used with repeated subscripts indicating summation unless otherwise noted The conservation equations are also presented in vector notation in 2 2 2 1 1 Momentum Conservation To begin the conservation of linear momentum is BH Cu cun Pa PP 3c tef 2 1 where u is the velocity 0 is the stress tensor p is the mass density and f is the body force The body force contribution p f typically accounts for buoyancy forces with f representing the acceleration due to gravity The stress may be written in terms of the fluid pressure and the deviatoric stress tensor as Oij pij tij 2 2 where p is the pressure j is the Kronecker delta and 7j is the deviatoric stress tensor A constitutive equation relates the deviatoric stress and the strain rate e g Tij 2uij i 2 3 3 CHAPTER 2 GOVERNING EQUATIONS Here the dynamic viscosity uij is represented as a second rank tensor with no summation on i j in Eq 2 3 Frequently the fluid viscosity is only available as an isotropic viscosity in which case the constitutive relation becomes Tij 2u i 2 4 The strain rate tensor is written in terms of the velocity gradients as 1
115. lso include the average temperature lt T gt and for species transport the mean species concentration lt Z gt Variable Name Quantity LS POST Label X velocity lt X Velocity gt Y velocity lt Y Velocity gt Z velocity lt Z Velocity gt Pressure lt p gt X vorticity lt X Vorticity gt Y vorticity Y Vorticity Z vorticity Z Vorticity Vorticity 2D vorticity Stream Function stream Temperature T Species 1 lt Z1 gt Species 2 lt Z2 gt Species 3 lt Z3 gt Species 4 lt Z4 gt Species 5 lt Z5 gt Species 6 lt Z6 gt Species 7 lt Z7 gt Species 8 lt Z8 gt Species 9 lt Z9 gt Species 10 lt Z10 gt Table 8 1 LS DYNA ISTATS 1 output variables and LS POST labels 8 2 2 Second Moment Statistics The computation of the second moment statistics includes the variables listed in Table 8 2 This data as well as the variables listed in Table 8 1 must be present in the input time average database d3mean in order for the statistics associated with ISTATS 2 to be computed Therefore selection of the level 2 statistics in LS DYNA includes the collection of all the level 1 time averaged data Incompressible Flow Solver 71 CHAPTER 8 FLOW STATISTICS Variable Name LS POST Label Velocity Correlations Velocity Pressure Correlations Velocity Temperature Correlations Velocity Species Correlations lt ujp gt lt u2p gt lt u3p gt lt u iT gt lt uT
116. lt values ITRE Turbulence Model Ci 2 C4 C3 C6 C7 Smagorinsky LESSGS Cs 01 1 1 94 Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 10 INITIAL_CFD Purpose Specify initial conditions for all nodal variables Card Format 1 of 3 emf fe fa fe fos fe fata Variable Default Card Format 2 of 3 PPP pole cae o ed ee Card Format 3 of 3 ms foe fe fa fe fs fe fate Variable Default Incompressible Flow Solver 95 CHAPTER 9 KEYWORD INPUT VARIABLE 96 U Im ud RHO Z2 Z10 EPS DESCRIPTION Initial x velocity Initial y velocity Initial z velocity Initial temperature Initial enthalpy Initial density Initial Species 1 concentration Initial Species 2 concentration Initial Species 10 concentration Intital turbulent kinetic energy Initial turbulent dissipation rate Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS 9 1 11 MAT CFD OPTION The MAT CFD cards allow fluid properties to be defined in a stand alone fluid analysis or in a coupled fluid structure analysis see CONTROL SOLUTION in the LS DYNA Keyword User s Manual The only OPTION currently available is CONSTANT This is material property type 150 It allows constant isotropic fluid properties to be defined for the incompressible low Mach CFD solver Card Format 1 of 3 coa a a a o or a Default ESTE AAA Card F
117. lue of concentration for species 1 The prescribed mass flux rate is Z Em ni f xit 2 19 where f xi t is the known mass flux rate through the boundary with normal n The prescribed flux rate may also be specified in terms of a mass transfer coefficient as Z pDi ij ox nj E hp Zi Zic 2 20 where ho is the mass transfer coefficient and Z is a reference species concentration The boundary conditions for the energy equation Eq 2 13 consist of a prescribed temperature or heat flux rate The prescribed temperature is TOS 2 21 and the prescribed heat flux rate is oT Mir Q xi t 2 22 where q is the known flux rate through the boundary with normal n The heat flux rate may also be prescribed in terms of a heat transfer coefficient oT NIS ni T Te Q 23 6 Incompressible Flow Solver 2 2 CONSERVATION EQUATIONS VECTOR NOTATION where h is the heat transfer coefficient and 7 is a reference temperature Initial conditions take on the form of prescribed velocity species and temperature distributions at t 0 Lgs ui x 0 up xi Z xi 0 Zp xi 2 24 T x 0 T G5 Remark For a well posed incompressible flow problem the prescribed initial velocity field in Eq 2 25 must satisfy Eq 2 25 2 26 see Gresho and Sani 6 If I 0 the null set i e enclosure flows with nju prescribed on all surfaces then global mass conservation enters as an
118. meters are used The sharp gradients in the temperature and velocity fields requires the use of advec tive flux limiters so ifct 1 for the momentum and scalar transport equations The pressure is solved with the SSOR preconditioned conjugate gradient method ipsol 22 and 5 A conjugate projection vectors nvec 5 Representative results for the momentum driven jet are shown in Figures 10 18 10 24 The kinetic energy time history in Figure 10 18 reveals a nearly linear increase in kinetic energy before the starting vortical structure leaves the outflow boundary at the top of the computational domain The nodal temperature time history plots in Figure 10 19 show the effect of the monotonicity preserving advective limiter on the temperature where a step change in the temperature is convected vertically during the starting transient The nodal time history plots in Figure 10 20 show a nearly periodic multi frequency response that is characteristic of the vortex merging process in the shear layer Snapshots of velocity temperature and pressure contours are shown in Figures 10 21 10 24 Example momjet p2 k KEYWORD TITLE Momentum driven jet 9 5 99955955559952252592559555599525259555552552555552552555555555255955555259555 Setup the CFD problem soln 4 CONTROL SOLUTION soln 4 CONTROL CFD AUTO iauto epsdt dtsf dtmax Incompressible Flow Solver 121 CHAPTER 10 EXAMPLE PROBLEMS 2 CONTROL CFD GENERAL insol dtinit cf
119. mpotent i e P P and Q Q Therefore repeated application of the projection operators does not continue to modify the projected results The projection operators are orthogonal and commute i e FO OPE The explicit forms of the continuous projection operators are P 1 V V v gt 4 11 and QC 1 P V V V 4 12 28 Incompressible Flow Solver 4 1 SEMI IMPLICIT PROJECTION METHOD It should be noted that and Q have built in all the appropriate physical boundary conditions Further The eigenvalues of P and Q are either 0 or 1 so that the projections are norm reducing In practice the action of the projection P is to remove the part of the approximate velocity field that is not div free i e u P In effect the projection is achieved by decomposing the velocity field into div free and curl free components using a Helmholtz decomposition The decomposition may be written as u4 V 4 13 where is a non solenoidal velocity field u is its div free counterpart and V is the curl free component ie V x V 2 0 4 1 Semi Implicit Projection Method In LS DYNA the Projection 2 P2 method identified by Gresho 3 forms the starting point for a discus sion of the finite element projection algorithms Before proceeding with a description of the P2 algorithm the semi discrete Navier Stokes equations are presented The spatial discretization of the conservation equations is achieved using
120. n 92 have determined that Cz 0 067 and Cg 0 931 based on the an analysis of the inertial viscous spectrum and SGS closure Unlike the Smagorinsky model which relates the eddy viscosity to the instantaneous filtered strain rate the k 8 model includes some effects of the flow history in the estimation of the eddy viscosity In addi tion there are no explicit assumptions equating the production and dissipation of turbulent kinetic energy making it possible to account for non equilibrium effects However this model still requires the use of flow specific constants In the following section the dynamic procedure for evaluating the constants is outlined 7 4 1 Local Dynamic k model The local dynamic k model was was first proposed by Kim and Menon 93 and later extended to compressible flows by Nelson 94 The dynamic evaluation of Cz and Cg makes these quantities functions of both space and time The dynamic evaluation procedure see for example 95 96 relies on a test scale or high pass filter with a filter scale A gt A where typically amp 2 In the ensuing discussion a circumflex denotes application of the test scale filter Application of the test scale filter further decomposes the flow variables into test scale and subtest scale components e g the subtest stresses are Tij p ui uji 7 31 Application of the test scale filter to the subgrid scale stresses yields a relationship between the
121. ncentration for species 1 5 3 Traction Boundary Conditions The formulation of the weak form of the momentum equations Eq 2 1 yields traction boundary condi tions also known as natural boundary conditions In terms of the shape functions N4 the traction boundary conditions may be computed as F E Na fi xi t dT 5 1 DL where fi are the components of the prescribed traction Here 1 lt i lt Ndim and 1 lt a lt Nnpe where Ndim is the number of space dimensions and Nnpe is the number of nodes per element Alternatively the boundary conditions may be written in terms of the stress e Fe Na in dT 5 2 i where 6 is the prescribed stress and n is the outward normal for the domain boundary In terms of the pressure and strain rate the traction boundary conditions are Fr fro pd j 2u amp jnj 5 3 By default homogenous traction boundary conditions are applied unless other boundary conditions are prescribed i e homogeneous traction natural boundary conditions are the do nothing boundary condi tions The traction and velocity boundary conditions can be mixed In a two dimensional sense mixed boundary conditions can consist of a prescribed normal traction and a tangential velocity For example at the outflow boundary in Figure 2 1 a homogeneous normal traction and vertical velocity on I con stitutes a valid set of mixed boundary conditions A detailed discussion of boundary conditions for the incompre
122. nd at is an approximate discrete velocity field at time n 1 Note that the discrete divergence of is generally not zero i e GL 0 where GT is the discrete divergence operator The functional dependence of F upon the discrete velocity u depends upon whether the algorithm is implicit explicit or semi implicit The dependence on pressure or rather pressure gradient is p F u 4 5 a u w w Lg At At p 4 6 where G is the discrete gradient operator and GT is the discrete divergence operator Applying the discrete divergence operator to Eq 4 6 yields a Poisson equation for the pressure at time level n 1 1 1 E cod you 4 7 By eliminating the velocity at time level n Eq 4 6 yields a relationship for the projected div free velocity field 1 ut grt DES 4 8 Projection Properties The philosophy behind projection algorithms is to provide a legitimate way to decouple the pressure and velocity fields in the hope of providing an efficient computational method for transient incompressible flow simulations Thus given an approximate non solenoidal velocity field i 7 may be projected onto a divergence free subspace such that ou 7 p PU 8 4 9 and Vp Q F 4 10 Here P and Q are the projection operators and they have the following properties P projects a velocity vector onto a div free subspace and Q projects a vector into a curl free subspace Both P and Q are ide
123. ness is a measure of the asymmetry of the fluctuations So gt 0 are predominately positive and So lt 0 are predominately negative The flatness is a relative measure of remotely occurring symmetric spiking fluctuations 68 Incompressible Flow Solver 8 1 DERIVED FLOW STATISTICS The skewness of velocity derivatives is defined as E Sii E 8 17 and the flatness factor of velocity derivatives du W Fi 2 8 18 Ou ey where there is no summation of repeated indices Skewness and flatness of velocity differences have been also measured 98 Au 3 Sai r A x ete Uj and the flatness factor is Au Fair Bu 8 20 A where Au u x r t ui xi t and rj xi r For homogeneous and isotropic turbulence the velocity fluctuations have a Gaussian probability distribu tion Since this distribution is symmetric all odd moments are zero The nonzero moment higher than the variance is the flatness that obtains a value of 3 0 The derivative skewness is 0 3 0 5 and the derivative flatness is 3 4 98 Velocity derivative skewness and flatness are dominated by the smaller scales in the flow so these statistics provide a good measure of how well the numerical scheme is resolving the smaller scales of the flow 8 1 4 Anisotropic Stress Tensor An intrinsic distinction exists between isotropic and anisotropic Reynolds stresses The isotropic stress is 3q 6 and the deviatoric anisot
124. nitialized material by material resulting in built in temperature initial conditions 5 6 Flux Boundary Conditions The heat flux rate at a boundary is prescribed as Qa Na i xi t ni dT 5 8 T where i is the known flux rate through the boundary with normal n Using the constitutive relationship between temperature and heat flux is oT ia ox j where kj is the thermal conductivity tensor The homogeneous form of the heat flux boundary condition is the default and represents a perfectly insulated i e zero heat flux boundary Gi xi t k 5 9 The heat flux rate may also be prescribed in terms of a convective heat transfer coefficient qx t nj h T T 5 10 where h is the heat transfer coefficient and To is a reference temperature For the species transport equations the prescribed mass flux rate for species 1 is me no fh n dT 5 11 T where fi is the known flux rate through the boundary with normal n For simplicity in the presentation only the boundary conditions for species 1 are presented here Using the constitutive relation between species concentration and mass flux rate 2 OZ fi i t epa 5 12 J where 7D is the tensorial mass diffusivity The homogeneous form of the mass flux boundary condition is the default and represents a boundary where the gradient of the species in the boundary normal direction is Zero The prescribed flux rate may also be specified in terms of a convect
125. nly referred to as fractional step pressure correction methods or Chorin s method 28 have grown in popularity over the past 10 years due to the relative ease of implementation and computational performance This is reflected by the volume of work published on the development of second order accurate projection methods see for example van Kan 29 Bell et al 30 Gresho et al 3 4 31 32 Almgren et al 33 34 35 Rider 36 37 38 39 Minion 40 Guermond and Quartapelle 41 Puckett et al 42 Sussman et al 43 and Knio et al 44 The numerical performance of projection methods has been considered by Brown and Minion 45 46 Wetton 47 Guermond 48 49 Guermond and Quartapelle 50 51 and Almgren et al 52 As background a brief review of Chorin s original projection method is presented before proceeding with the finite element form of the projection algorithm The vector form of the momentum equations may be written as ou Pa Vp Flu 4 1 where for a constant viscosity F u f uV u pu Vu 4 2 Now F u may be decomposed into a div free and curl free part where the div free part is ou Uode 4 3 l 3 4 3 and the curl free part is VxVp 0 4 4 Discretizing in space and time the decomposition neglecting the contribution of the pressure gradient 27 CHAPTER 4 THE PROJECTION METHOD yields u At where F u is the spatially discrete analogue of F in Eq 4 2 a
126. nowledge Nielen Stander Ian Do Jason Wang Phil Gresho and Steve Sutton at Lawrence Livermore National Laboratory and Tom Voth at Sandia National Laboratories for their helpful suggestions during the preparation of this manual In addition I would like to thank Suresh Menon at Georgia Tech for his assistance with the k LES models and Sara Lucas for her help with the turbulence statistics Finally I would like to thank John Hallquist for providing the opportunity to contribute to LS DYNA s growing set of simulation capabilities Mark A Christon Albuquerque NM March 2001 vil viii Chapter 1 Introduction The incompressible flow capability in LS DYNA has been developed specifically to attack the class of transient incompressible viscous fluid dynamics problems that are predominant in the world that sur rounds us The goal for for the flow solver has been to achieve high performance across a spectrum of supercomputer architectures without sacrificing any of the aspects of the finite element method that make it so flexible and permit application to a broad class of fluid dynamics problems The incompressible flow solver plays two primary roles in LS DYNA The first is to provide a stand alone incompressible flow solver that complements the existing solid structural boundary element and compressible flow capabilities that comprise LS DYNA s multi physics capabilities The second is to provide the flow solver for fluid soli
127. ns for a variety of flow related physics in LS DYNA The selection of the incompressible low Mach flow physics and related flow solver is determined by the INSOL input on the CONTROL_CFD_GENERAL keyword Currently there are two valid values for IN SOL INSOL 1 selects the explicit time integrator that requires the use of a lumped mass matrix In this case the IMASS THETAK THETAB THETAA and THETAF variables associated with the CONTROL CFD MOMENTUM keyword are ignored INSOL 3 selects the semi implicit DESCRIPTION Set the interval to check and report the grid Reynolds and advective CFL numbers ICKDT lt 0 checks and reports the grid Reynolds and advective CFL numbers but does not modify the time step ICKDT gt 0 modifies the time step according to the prescribed CFL limit and any required stability limits The report of the grid Reynolds and CFL numbers to the screen may be toggled with the grid sense switch EQ 0 ICKDT 10 default Activate the use of full quadrature for certain terms in the momentum and transport equation The accuracy flag improves the accuracy of body force calculations and certain advective convective terms with a modest increase in computational time EQ 0 don t use the increased quadrature rules default EQ 1 use increased quadrature on advective convective body force terms projection algorithm which makes use of these variables 2 This option is only available for the semi implicit implicit solution
128. nt computation of turbulent flow Technical Report TF 5 Thermosciences Division Department of Mechanical Engineering Stanford University Stanford California May 1975 M Germano Turbulence the filtering approach Journal of Fluid Mechanics 238 325 336 1992 J Smagorinsky General circulation experiments with the primitive equations i the basic experi ment Monthly Weather Review 91 3 99 164 1963 Robert S Rogallo and Parviz Moin Annual Review of Fluid Mechanics volume 16 chapter Nu merical simulation of turbulent flows pages 99 137 Annual Reviews Inc Palo Alto California 1984 Won Wook Kim A new dynamic one equation subgrid scale model for large eddy simulations Technical Report AIAA 95 0356 AIAA 33rd Aerospace Sciences Meeting and Exhibit Reno NV January 1995 Won Wook Kim and Suresh Menon An unsteady incompressible navier stokes solver for large eddy simulation of turbulent flows International Journal for Numerical Methods in Fluids 31 983 1017 1999 U Schumann Subgrid scale model for finite difference simulation of turbulent flows in plane channels and annuli Journal of Computational Physics 18 376 404 1975 A Yoshizawa Bridging between eddy viscosity type and second order models using a two scale dia In Ninth Symposium on Turbulent Shear Flows Kyoto Japan August 16 18 pages 23 1 1 6 1993 Helmut Schmidt and Ulrich Schumann Coherent structure of the convective boundary layer
129. of the energy spectrum is resolved is quite high For example for simple channel flow the grid resolution for pure LES scales approximately as Re for simple channel flow see for ex ample 75 76 In contrast the idea of under resolved or VLES very large eddy simulation has shown promising results in a variety of applications see for example 77 78 79 80 81 82 suggesting that the non purist view of LES may provide useful engineering results LS DYNA s CFD capabilities are intended to find broad application and since there is no single universal turbulence model that is suitable for all applications the goal is to provide multiple turbulence models and modeling approaches In the subsequent sections of this chapter the averaging approaches and current turbulence models available in LS DYNA are outlined Refer to Chapter 9 for details on the CON TROL_CFD_TURBULENCE keyword which is used for the specification of turbulence model constants and model options 57 CHAPTER 7 TURBULENCE MODELS 7 1 Averaging Filtering and Scale Separation The broad range of length and time scales associated with the highly disordered three dimensional and time dependent behavior of turbulent flow makes direct simulation for high Reynolds number engineering flows computationally intractable In order to attempt to model the behavior of turbulent flows the fields are typically separated into mean and fluctuating components The specific type of av
130. onal to the sub domain size but cannot smooth error modes that span sub domains Operational experience with both the global and local pressure stabilization formulations has shown that either approach is effective in terms of filtering spurious pressure modes and improving the overall robust ness of the computations Ultimately the combination of the pressure stabilization with preconditioning and the A conjugate projection technique has proven both robust and computationally efficient making it a reasonable alternative to more complicated approaches based on multi grid or multi level algorithms Incompressible Flow Solver 53 CHAPTER 6 PRESSURE SOLUTION METHODS Table 6 3 Effect of the A Conjugate projection on the PPE solution time for the 3 D square cylinder vortex shedding problem The grind times for the explicit algorithm have been normalized with respect to the solution time for 1000 time steps using a direct solver 71 72 N77 is the average number of iterations required per time step to solve the PPE problem Vel 22380 Nnp 25568 N 772 Explicit Algorithm Global Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG AAA EEG No of Grind Grind Grind n Nir Time 7 gt n or 7 cl Nir ame gt Explicit Algorithm Local Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG No Nea Grind Grind Grind vectors TT wr 7 DE Mtr a 7 a Mir ome 7 PE 54 Incompressible Flow Solver 6 4 PPE SOLVER PERFO
131. onditions are pI 2uS n f x r on T 2 44 Incompressible Flow Solver 9 CHAPTER 2 GOVERNING EQUATIONS The traction and velocity boundary conditions can be mixed In a two dimensional sense mixed boundary conditions can consist of a prescribed normal traction and a tangential velocity For example at the outflow boundary in Figure 2 1 a homogeneous normal traction and vertical velocity on I constitutes a valid set of mixed boundary conditions A detailed discussion of boundary conditions for the incompressible Navier Stokes equations may be found in Gresho and Sani 5 Turning attention to the species transport equations boundary conditions for Eq 2 9 may consist of either a prescribed concentration or a mass flux rate In the binary mixture example the prescribed con centration is Z x t Z1 x t 2 45 where Z is the known value of concentration for species 1 The prescribed mass flux rate is pDiVZ n f x t 2 46 where f xi t is the known mass flux rate through the boundary with normal n The prescribed flux rate may also be specified in terms of a mass transfer coefficient as pD VZ n hp Zi ZA 2 47 where hp is the mass transfer coefficient and Z is a reference species concentration The boundary conditions for the energy equation Eq 2 13 consist of a prescribed temperature or heat flux rate The prescribed temperature is T x t f x t 2 48 and the prescribe
132. ormat 2 of 3 CAMERA Incompressible Flow Solver 97 CHAPTER 9 KEYWORD INPUT Card Format 3 of 3 ems x foe a fe fos fis for fs VARIABLE MID RHO MU K CP BETA TREF GX GY GZ DIFF1 DIFF2 DIFF3 DIFF10 98 DESCRIPTION Material identification a unique number has to be chosen Fluid density Fluid viscosity Thermal Conductivity Heat capacity Coefficient of expansion Reference temperature Gravitational acceleration in x direction Gravitational acceleration in y direction Gravitational acceleration in z direction Species 1 diffusivity Species 2 diffusivity Species 3 diffusivity Species 10 diffusivity Incompressible Flow Solver 9 2 SHARED KEYWORDS 9 2 Shared Keywords The incompressible flow solver is an integral part of LS DYNA and because of this it shares keywords with the other physics options available in LS DYNA This section identifies the shared keywords and defines their interaction with the flow solver CONTROL SOLUTION Specify the analysis solution procedure to be used Here SOLN 4 activates the incompressible flow solver in a stand alone mode CONTROL OUTPUT Set miscellaneous output options For the flow solver NPOPT 0 echos flow velocities and pressures to the d3hsp file CONTROL TERMINATION Stop the job The EDENG and ENDMAS parameters are reserved for structural analyses DATABASE BINARY D3PLOT Set the time interval for the output of s
133. ormed according to Eq 3 42 r ave f Keu A u uz 3 41 F a FP 3 42 There are several things to note about the parallel assembly procedure First the parallel assembly is simply a generalization of the sequential assembly procedure that includes inter processor communication Second the algorithm only requires the communication of nodal data at the edges of adjacent sub domains Therefore as the problem size increases the communication overhead scales with number of surface nodes associated with sub domain boundaries Finally this algorithm permits the implementation of vector valued messages in order to avoid start up issues associated with short messages i e for the assembly shown in Eq 3 42 the message length is proportional to the number of sub domain boundary nodes and the number of degrees of freedom per node Finally the use of non overlapping grids implies that nodes in the finite element mesh that lie on sub domain boundaries are stored in multiple processors as shown in Figure 3 4b In contrast over lapping sub domains would require the redundant storage of all the data associated with elements at the sub domain boundaries 3 2 3 The Parallel Explicit Algorithm In this section the issue of solving the PPE in parallel will not be addressed so that attention may be fo cused upon the solution process for the nodal variables The DDMP version of the explicit time integration algorithm proceeds as follow
134. ourglass stabilization EQ 2 y hourglass stabilization viscous form Set the hourglass stabilization multiplier see IHG above EQ 0 EHG 1 0 default Incompressible Flow Solver 9 1 INCOMPRESSIBLE FLOW KEYWORDS Remarks 1 The IMASS variable is only active when INSOL gt 2 on the CONTROL CFD GENERAL keyword 2 The balancing tensor diffusivity should always be used with the explicit treatment of the advection terms This is the default 3 The use of the flux limiting procedures is currently restricted to the explicit advective procedures 4 The time weighting variables only apply to the case when INSOL gt 2 on the CON TROL_CFD_GENERAL keyword 5 The ITSOL keyword for the CONTROL_CFD_TRANSPORT keyword only applies for INSOL gt 2 on the CONTROL CFD GENERAL keyword Incompressible Flow Solver 93 CHAPTER 9 KEYWORD INPUT 9 1 9 CONTROL CFD TURBULENCE Purpose Activate a turbulence model and set the associated model parameters Card Format Ere a a os fe for fe ents pa ao VARIABLE DESCRIPTION TTRB Select the turbulence model EQ 0 Turbulence models are disabled default EQ 1 Smagorinsky LES subgrid scale model C1 C7 Turbulence model constants See the table below for the model specific definition of the constants and their default values Remarks 1 The number and value of the input constants are model specific See the table below for definition of the constants and their defau
135. p 1 0 diff2 diff10 propt nloc beta Tref 710 0 0 5 diff3 diff4 diff5 qr irid icomp Incompressible Flow Solver 10 2 NATURAL CONVECTION IN A SQUARE CAVITY 1 0000000 1 0000000 1 0000000 1 0000000 PART Fluid domain pid secid mid eosid hgid grav adpopt tmid 1 1 1 0 0 0 0 0 8 9 9999999999999999 9 9 9999 9SSSSSSEP UR SS ooo oS SS SSS SS SS USESESUUUSUPSSSSESUUUDETS 5 Handle output of state time history and restart data DATABASE BINARY D3PLOT dt cycl ledt beam nplte 1 0 DATABASE BINARY D3THDT dt cycl ledt beam npltc 1 0e 6 DATABASE HISTORY NODE ndl nd2 nd3 nd4 nd5 nd6 nd7 nd8 49 115 291 162 392 167 227 397 DATABASE BINARY RUNRSF dt cycl ledt beam npltc 1000 99 352553352 55252 25552223255525252525225552255225522525252552255522252525522555255 INCLUDE box dat k END Example box p2 k Sample screen output from the fully implicit computation is shown in Figure 10 10 This output reports the minimum maximum pressure velocities and temperature as well as the current local time truncation error LTE new and old time steps and the limiting minimum maximum nodal differences that the local time truncation error is based on Although time accurate methods were used for this computation a steady state solution is obtained in each case The time history plot of the kinetic energy verifies this as shown in Figure 10 11 The minimum and maximum time steps used by each method over the com
136. pmx 0 320E 01 0 595E 00 1 mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 0931E 02 div u 3 7186E 16 pmn pmx 0 320E 01 0 595E 00 1 mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 0942E 02 div u 3 7329E 16 pmn pmx 0 320E 01 0 595E 00 1 mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 0954E 02 div u 3 7428E 16 pmn pmx 0 320E 01 0 595E 00 a mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 2 0966E 02 div u 3 7498E 16 pmn pmx 0 320E 01 0 595E 00 1 mn mx 0 000E 00 0 150E 01 261E 07 0 294E 07 0966E 02 dt 1 19E 01 write d3plot file 0966E 02 dt 1 19E 01 write runrsf file c 1996 er N H 1997 c N H 1998 c H N 1999 er 3p 2000 E H 2000 2000 co cr cH OS ct Ste tO a CF cat NONI Figure 10 3 Screen output of RMS divergence minimum maximum pressure and velocity values Incompressible Flow Solver 107 CHAPTER 10 EXAMPLE PROBLEMS 0002 62 11 95 15 61 mduoo D X velocity SIMPLE DUCT ENTRANCE REGION Node No e A6 qe 8 85 a 226 0 8 0 6 0 4 0 2 50 100 150 Time 200 Figure 10 4 Nodal time history plot for the x velocity along the centerline at the inlet node 5 and the outlet node 226 108 0002 22 11 ooz Li Aseo O Kinetic Energy SIMPLE DUCT ENTRANCE REGION 3 5 2 5 N A a 0 5
137. putation are summarized with the peak velocities in Table 10 1 The maximum time step for the fully implicit method was limited to 1 76e 3 while the semi implicit method found a maximum time step nearly 4 times larger This is due to the fact that the local time truncation error was limited to 1 0e 3 for the fully implicit method while the time step calculation for the semi implicit method was based solely on the convective stability That is to say the semi implicit method did not have the same temporal accuracy as the fully implicit method The peak velocities computed here compare favorably with the velocities reported in Table I in 100 for a 40 x 40 grid despite the fact that a grid with 20 x 20 elements 21 x 21 nodes was used for these computations Stream function temperature velocity and vorticity contour plots are shown in Figures 10 12 10 16 Incompressible Flow Solver 115 150 151 initialization completed 0 t 0 0000E 00 u_i mn mx 0 tmn mx 0 0 t 0 0000E 00 div u 0 000E 00 0 000E 00 0 1 t 1 0000E 04 u_1 mn mx 434E 02 0 tmn mx 0 000E 00 0 LTE 1 7212E 04 DTSF Min diff 7 6871E 04 Max diff 7 6870E 04 div u 1 2 t 2 5000E 04 div u 3 u_1 mn mx 907E 02 0 tmn mx 0 000E 00 0 LTE 2 0413E 04 DTSF Min diff 8 4254E 04 Max diff 8 4254E 04 3 t 4 7500E 04 div u 7 u i mn mx 166E 01 0 tmn mx 0 000E 00 0 LTE 2 9651E 04 DTSF Min diff 1 3321E 03 Max diff 1 3321E 03 t 1
138. putational Physics 138 734 765 1997 47 Brian B Wetton Error analysis of pressure increment schemes submitted to SINUM April 1998 48 Jean Luc Guermond Some implementations of projection methods for navier stokes equations Mathematical Modelling and Numerical Analysis 30 5 637 667 1996 49 Jean Luc Guermond A convergence result for the approximation of the navier stokes equations by an incremental projection method C R Acad Sci Paris 325 1329 1332 1997 50 Jean Luc Guermond and L Quartapelle On the approximation of the unsteady navier stokes equa tions by finite element projection methods Numerische Mathematik 80 207 238 1998 51 Jean Luc Guermond and L Quartapelle On stability and convergence of projection methods based on pressure poisson equation nternational Journal for Numerical Methods in Fluids 26 1039 1053 1998 52 Ann S Almgren John B Bell and William Y Crutchfield Approximate projection methods Part i inviscid analysis SIAM Journal on Scientific Computing 22 4 1139 1159 2000 53 J C Heinrich S R Idelsohn E Onate and C A Vionnet Boundary conditions for finite element simulations of convective flows with artificial boundaries International Journal for Numerical Methods in Fluids 39 1053 1071 1996 54 Robert L Sani P M Gresho Robert L Lee and D F Griffiths The cause and cure of the spurious pressures generated by certain fem solutions of th
139. quires a unit density with viscosity u 1 Re The DATABASE BINARY OPTION keywords are used to write a graphics state to the d3plot every 10 time units nodal time history every time step and a running restart every 1000 time steps The input files duct dat k duct p2 k and duct fe k may be found in the duct sub directory of the 2 D examples Example duct p2 k KEYWORD TITLE Simple duct entrance region 5 O OO 5 102 Incompressible Flow Solver 10 1 POISEUILLE FLOW Setup the CFD problem soln 4 CONTROL SOLUTION soln 4 CONTROL CFD AUTO iauto epsdt dtsf dtmax 2 CONTROL CFD GENERAL insol dtinit etl ickdt istats tstart iavg 3 0 010 2 0 100 CONTROL CFD MOMENTUM imass iadvec ifct divu thetak thetaa thetaf 1 10 1 0 msol maxit ichkit iwrt ihist eps ihg ehg 20 100 2 0 0 CONTROL CFD PRESSURE ipsol maxit ickit iwrt ihist eps 11 nvec stab beta sid plev lcid CONTROL TERMINATION endtim endcyc dtmin endeng endmas 1000000 00 2000 5 PAPA 5 Define the load curves 5 DEFINE CURVE lcid sidr sfa sfo offa offo dattyp 1 al o1 0 0 1 0 1000000 0 150 DEFINE CURVE lcid sidr sfa sfo offa offo dattyp 2 al o1 0 0 0 0 1000000 0 0 0 5 AAA AAA 5 Setup the boundary conditions 5 Pressure boundary conditions 5 BOUNDARY PRESSURE CFD SET Posee tl psa A toy SRA OS 4 5 Peers 7 8 ssid lcid p 100 T 0 6 Incompressible Flow Solver 103 CHAPTER 10 EXAMPLE PR
140. required per time step to solve the PPE problem Nel 11250 Nnp 11466 N 146 direct solver grind time 4 871 x 107 CPU seconds per element per time step Projection 2 Algorithm No Stabilization JPCG SSOR PCG ESSOR PCG No Nes Grind Grind Grind Vectors AIT nu T 2 TIT T I TT ame ua gt Projection 2 Algorithm Global Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG Ly No of Grind Grind Grind vectors Nir nme a T ore T E Mer Time i PE Projection 2 Algorithm Local Jump Stabilization B 0 05 JPCG SSOR PCG ESSOR PCG AA No of Grind Grind Grind vectors a nu T us IT ome T T Mer ame ua Lo 52 Incompressible Flow Solver 6 4 PPE SOLVER PERFORMANCE Z Vorticity 5 00e 01 gt 4 00e 01 gt 2 00e 01 gt 0 00e 00 2 00e 01 4 00e 01 5 00e 01 a Pressure 2 96e 00 2 00e 00 gt 1 00e 00 0 00e 00 7 52e 01 b Figure 6 6 Snapshot of a cutplane showing the z vorticity field and c isosurfaces of the pressure for the 3 D cylinder in a channel for Rep 100 the A conjugate projection technique is to introduce both long and short wavelength information thus effectively deflating the eigenvalue spectrum of the PPE This effect not only reduces the iteration count and concomitant computational cost but is seen to improve the effectiveness of simple sub domain parallel preconditioners that are effective for wavelengths proporti
141. ressure field and therefore is not necessarily solenoidal Incorporating the pressure gradient contribution the approximation to F u p is ut u wt u i M M MM7 Cp 4 17 At At L EP B Incompressible Flow Solver 29 CHAPTER 4 THE PROJECTION METHOD Eliminating the velocity at time level n yields the discrete statement of the Helmholtz decomposition with the concomitant div free constraint C u 1 viz M C n 1 Mi Bae Val E H 4 18 Here M is the row sum lumped i e diagonalized mass matrix and A At pl p The system of equations in Eq 4 18 is analogous to those introduced by Chorin 28 Although Chorin suggested using a method of successive substitutions to obtain the velocity and pressure fields correspond ing to the div free state the solution of Eq 4 18 via a gradient based iterative method and the solution of the Schur complement of Eq 4 18 has proven more efficient The Schur complement of the projection operator may be formed explicitly from Eq 4 18 yielding IC M Gc anu rt 4 19 where Mz is the diagonalized i e row sum lumped mass matrix and includes the prescription of essential velocity boundary conditions Note that the term projection operator is used loosely here to describe Eq 4 18 since its solution yields a div free velocity field and the corresponding Lagrange multiplier However this operator is not to be confused
142. ropic part is 2 aij lt un gt 340j 8 21 The normalized anisotropy tensor is defined as bj 7 It is readily seen that aj 0 and any anisotropy is seen by non zero values of the off diagonal terms Incompressible Flow Solver 69 CHAPTER 8 FLOW STATISTICS Start Dump Start Dump Start Dump 0 50 100 150 200 250 Time Figure 8 2 Time windows for the collection of the raw time averaged data for LS POST 8 2 LS DYNA Statistics Levels Currently LS DYNA generates three levels of time averaged output data that range from mean quantities level 1 such as mean velocity pressure etc to an intermediate level level 2 that includes correla tions such as velocity and pressure and the highest level level 3 that includes the higher order corre lations lt u gt and lt ut gt The level of time averaging is selected using the STATS parameter on the DATABASE BINARY D3MEAN keyword At run time LS DYNA generates the d3mean data file that contains raw time averaged data that is col lected during the course of a simulation In order to facilitate the computation of a windowed time average for the derived statistics LS DYNA collects the raw time averaged data over a series of time windows as shown in Figure 8 2 Here the time averaging procedure starts at tstart 50 and is restarted after the data for a time window is written to disk This approach permits the exclusion of start up transients in t
143. s 1 Calculate the partial acceleration i e acceleration neglecting the pressure gradient at time level n i All processors calculate the processor local right hand side terms neglecting the pressure gradient according to Eq 3 41 ii Perform the inter processor finite element assembly according to Eq 3 42 In this step pro cessors that share a common sub domain boundary exchange messages containing partially assembled right hand side contributions and accumulate the fully summed terms at the sub domain boundaries iii All processors compute the processor local partial acceleration using a pre calculated lumped mass matrix i e the mass matrix has been fully summed for the nodes at the sub domain boundaries a M 3 43 24 Incompressible Flow Solver 3 2 THE DOMAIN DECOMPOSITION MESSAGE PASSING PARADIGM 2 Solve the global PPE problem for the current pressure field CT M c P ca 3 44 where g 0 has been used for simplicity 3 Update the nodal velocities i The computation of the pressure gradient term in Eq 3 45 consists of the parallel on processor computation of the discrete pressure gradient followed by an inter processor assembly u u Ari Mj 575 CP y 345 4 Repeat steps 1 3 until a maximum simulation time limit or maximum number of time steps is reached Remark For the explicit solution algorithm the global row sum lumped mass is accumulated at all nodes in the finite
144. scribed initial conditions and boundary conditions are tested and if necessary a projection to a divergence free subspace is performed on the initial velocity field u This guarantees that the flow problem is well posed even if the user prescribed initial conditions violate the conditions of Eq 2 25 2 26 presented in Chapter 2 In practice the criterion for performing a projection onto a div free subspace is based upon the RMS divergence error 14 Incompressible Flow Solver CTu CTu 3 8 Nel mS f where Nel is the number of elements and e is a user specified tolerance typically 10719 lt lt 1070 If the RMS divergence error is greater than the specified tolerance for the initial candidate velocity field then the PPE problem in Eq 3 9 is solved for A and a mass consistent projection performed using Eq 3 10 c M cjA c a 3 9 u MF C 3 10 The explicit time integration algorithm must respect both diffusive and convective stability limits Al though the analytical stability limits for the explicit time integration of the Navier Stokes equations in multiple dimensions remain intractable 1 approximate stability computations may be performed using local grid metrics In an unstructured grid with variable element size the calculation of the grid Re Reynolds and CFL Courant Freidrichs Levy numbers uses the element local coordinates and centroid velocities Figure 3 2 shows the canon
145. sense witch code For the incompressible flow solver the sw2 sw5 and swa switches are ignored The active flow solver sense switches are Type sw0 swl sw3 sw4 div grid iter Iprint quit Stop Response LS DYNA terminates immediately A restart file is written and LS DYNA terminates A restart file is written and LS DYNA proceeds A plot state is written and LS DYNA proceeds Toggle the output of the divergence metrics for a projection step Toggle the output of the grid Reynolds and CFL numbers Toggle the output of statistics for the iterative equation solvers Toggle the output of velocity min max LS DYNA terminates immediately LS DYNA terminates immediately 139 APPENDIX A SENSE SWITCH CONTROLS 140 Incompressible Flow Solver Bibliography 1 Philip M Gresho Stevens T Chan Robert L Lee and Craig D Upson A modified finite element 2 3 4 5 6 7 8 ras LL LL method for solving the time dependent incompressible navier stokes equations part 1 Theory International Journal for Numerical Methods in Fluids 4 557 598 1984 Philip M Gresho Stevens T Chan Robert L Lee and Craig D Upson A modified finite element method for solving the time dependent incompressible navier stokes equations part 2 Applica tions International Journal for Numerical Methods in Fluids 4 619 640 1984 Philip M Gresho On the theory of semi implic
146. solution of the PPE problem where the coefficient matrix is fixed and the right hand side changes each time step To address this aspect of solving the PPE an A conjugate projection is integrated with the iterative solution of the PPE in order to use solution information from the previous time steps In the ensuing discussion the PPE problem is cast as Ax b where A CT M ICorA CM IC S x handb CT grt 44 Incompressible Flow Solver 6 3 THE PROJECTION CG METHOD 1 a Un stabilized b Global Jump 1 2 lt 1 ara Q 1 Ne oq De M NEN NEN E MEE NICA Coarse Medium Fine 1 0 100 200 300 400 500 Iterations Figure 6 3 Convergence histories for SSOR preconditioned CG applied to the un stabilized PPE b the PPE with the global jump stabilization and c the PPE with the local jump stabilization for the coarse 2800 elements medium 11200 elements fine resolution 44800 elements vortex shedding mesh see Christon 63 The use of an A conjugate projection as a pre processing step for the solution of the linear system Ax b follows the development presented by Fischer 64 with extensions that permit seeding the A conjugate vectors Related work on solving linear systems with multiple right hand sides may be found in Saad 65 and Chan and Wan 66 To begin the development the idea of a pre processing A conjugate projection step relies on minimizing the distance between the solution at a given time step x
147. solver Incompressible Flow Solver 99 CHAPTER 9 KEYWORD INPUT 100 Incompressible Flow Solver Chapter 10 Example Problems This chapter presents several LS DYNA flow calculations for the first time user who wishes to perform verification computations for comparison before embarking on a full blown analysis For this reason the control files are replicated here with representative results that can be used to check the local LS DYNA installation For the most part sample problems presented here use relatively coarse meshes to minimize run times and provide a starting point for the user who wishes to experiment with code options before attempting any significant calculations The example problems in this chapter may be obtained by contacting Livermore Software Technology Corporation All of the example problems presented in this chapter use the following nomenclature The example input files are prefixed by the problem name e g cylinder File names that contain dat k as a suffix contain the nodal coordinates element connectivity and all sets for the problem The file names with a _p2 k or fe k suffix contain the primary control keywords for the semi fully implicit projection and explicit solution methods respectively For example cylinder p2 k will contain the primary control keywords for the cylinder flow problem and will INCLUDE the cylinder dat k data file Details on the INCLUDE keyword may be found in the LS DYNA Keyword User s Man
148. ssible Navier Stokes equations may be found in Gresho and Sani 5 Note that this boundary condition provides an important component for the coupling between fluid and structural components for coupled problems and is computed internally for this class of problems 5 3 1 Pressure Boundary Conditions In many practical situations the viscous contributions in Eq 5 3 may be neglected This is typically the case in situations where the viscosity is small i e the Reynolds number is relatively large In this case the viscous terms are ignored leaving only the pressure contribution to the traction boundary conditions 5 Na B m 5 4 T Incompressible Flow Solver 37 CHAPTER 5 BOUNDARY CONDITIONS AND SOURCE TERMS where f is the prescribed pressure Due to the common occurrence of this type of limiting traction boundary condition LS DYNA refers to this as a pressure boundary condition which reflects its use in finite element models The BOUND ARY PRESSURE CFD SET keyword is used to prescribe the pressure as a function of time A t on a boundary segment set 5 4 Outflow Boundary Conditions In many problems it is necessary to truncate the physical domain resulting in an artificial outflow bound ary In the presence of strong vortical structures the homogeneous natural do nothing boundary condi tions can result in a global pressure response that is physically unrealistic Heinrich et al 53 present a summary of this problem
149. sure increment highlighting one of the primary differences between the projection and the explicit methods The explicit method segregates acceleration and pressure which can be done legitimately while the projection methods attempt to decouple velocity and pressure Similar to the projection methods the pressure Poisson problem in Eq 3 6 is simply the Schur complement of Eq 6 2 For both algorithms the PPE and the saddle point operators are symmetric and in general singular due in part the fact that the pressure is known only to within a constant representing the hydrostatic pressure 6 2 0100 Stabilization Although the Q1Q0 element has been condemned by theoreticians for its weakly singular modes this element has long been the workhorse for incompressible flow and continues to be widely used The unstable modes of the 0100 element for incompressible flow have been investigated by Sani et al 54 55 and more recently by Griffiths and Silvester 56 Griffiths and Silvester demonstrate that for problems of physical interest the Q1Q0 element will converge to the true solution in the limit as h gt 0 A new convergence proof for the Q1Q0 element may be found in Gresho and Sani 5 The ensuing discussion summarizes the pressure stabilization techniques in LS DYNA that are compat ible with the Q1Q0 element and that circumvent the usual BabuSka Brezzi div stability condition In LS DYNA there are two pressure stabilization techniques one t
150. t universal the eddy viscosity does not vanish in laminar regions the limiting wall behavior is incorrect and the model is excessively dissipative for transitional flows This section outlines the LES models in LS DYNA that are based on the subgrid scale kinetic energy The SGS kinetic energy is xz s ub EE 3 uu UU 7 26 Following Menon and Kim 87 88 the transport equation for the SGS kinetic energy is 9k385 B OKS85 o a Ou P Ly A EE LLL _ 9 585 E pe Ala or 5 p ox Ox EA Ox Ox eel where e is the subgrid scale dissipation This transport equation is based on the difference between the equation of mechanical energy and its spatially filtered form and it is similar to the k transport equation of Schumann 89 and follows the form used by Yoshizawa 90 62 Incompressible Flow Solver 7 4 SUBGRID KINETIC ENERGY MODELS In Eq 7 27 the SGS stresses are modeled using the Boussinesq hypothesis with the filter scale A providing the length scale and the SGS kinetic energy providing the velocity scale Thus the SGS stresses are VP 2CAV KS Si E FUE 7 28 where t 2K 55 and the eddy viscosity is Ly CAV RS 7 29 In order to close Eq 7 27 the dissipation rate is modeled as ses 3 2 e ES UE 7 30 In order to solve Eq 7 27 the constants C and Cz must be prescribed Schmidt and Schumann 91 suggest Cz 0 0856 and Cz 0 845 Alternatively Chakravarthy and Meno
151. tate databases for plotting purposes For incompressible flow calculations LS DYNA generates an augmented state database that contains flow specific variables e g vorticity stream function etc DATABASE BINARY D3THDT Set the time interval for the output of time history data For incom pressible flow calculations LS DYNA generates a separate binary time history database that contains flow specific variables e g vorticity stream function etc This new database is the d3thins database DATABASE BINARY D3MEAN Define the starting time the time window in time steps and the level of statistics to be used in creating a time averaged database of flow statistics DATABASE BINARY D3DUMP Define the output frequency in cycles for the binary restart files DATABASE BINARY RUNRSE Define the output frequency in cycles for the running binary restart files DATABASE HISTORY NODE Control which nodes are output to the binary database D3THDT and the ASCII file VODOUT DEFINE CURVE Define a load curve to be used for the prescription of time dependent boundary conditions SECTION SHELL Define the section properties for 2 D shell elements This keyword is used for 2 D Cartesian flow problems to select the element formulation used with the Navier Stokes solver SECTION SOLID Define the section properties for the 3 D elements This keyword is used for 3 D flow problems to select the element formulation used with the Navier Stokes
152. the Boussinesq approximation The initial conditions consist of a div free velocity field with a free field temperature of 300K an inlet temperature of 400K The working fluid is air with a unit Prandtl number Pr v a resulting in Pe 4000 where Pe RePr Here the relatively large Froude number indicates that the influence of buoyancy forces is small compared to the inertial forces i e a momentum driven jet Figure 6 5 shows snapshots of the temperature vorticity and pressure fields from the momentum driven jet simulation In this computation natural boundary conditions were applied at the far field boundary removing the hydrostatic mode from the PPE Table 6 1 and 6 2 show the results of the jet computations in terms of the average iteration count and grind time for 1000 time steps Here the grind times are normalized with respect to the grind time for the same computation using the PVS direct solver 71 72 where a single factorization is performed during initialization and a single re solve is carried out at each time step All computations in this comparison were performed using a 500 MHz single processor DEC Alpha The solution time due to the PPE re solve was 15 2 for the semi implicit projection method and 60 0 for the explicit algorithm using the direct solver This relatively low percentage is due to the fact that a single off line factorization is performed with one re solve per time step For the PPE the ideal situation would be
153. the Q1Q0 element with bilinear support for velocity and piecewise constant support for the pressure in two dimensions In three dimensions the velocity support is trilinear with piecewise constant support for pressure The methods for obtaining the weak form of the conservation equations are well known and will not be repeated here see for example Gresho et al 7 Hughes 8 and Zienkiewicz and Taylor 9 The spatially discrete form of momentum conservation Eq 2 1 and the divergence constraint Eq 2 7 are Mu A uju Ku MM Cp F 4 14 Clu g 4 15 where M is the unit mass matrix A u and K are the advection and the viscous diffusion operators respec tively F is the body force and g accounts for the presence of prescribed velocity boundary conditions C is the gradient operator and C is the divergence operator i e C and CT are nearly the discrete finite element analogues of G and G7 discussed above In order to simplify the nomenclature u and p are understood to be discrete approximations to the continuous velocity and pressure fields Following the development of the projection method introduced above qt u M F u p 4 16 where 7 u p may involve an explicit or implicit dependence upon u and is understood to be a discrete vector quantity From this it is clear that the intermediate velocity corresponds to an approximate velocity field that has not yet felt the influence of the current p
154. tional Fluid Dynamics March 1996 LLNL UCRL JC 123727 33 Ann S Almgren John B Bell Phillip Colella and Louis H Howell An adaptive projection method for the incompressible euler equations In Eleventh AIAA Computational Fluid Dynamics Confer ence pages 530 539 AIAA 1993 34 Ann S Almgren John B Bell and William G Szymcyzk A numerical method for the incom pressible navier stokes equations based on an approximate projection SIAM Journal for Scientific Computing 17 2 358 369 March 1996 35 Ann S Almgren John B Bell Phillip Colella Louis H Howell and Michael L Welcome A con servative adaptive projection method for the variable density incompressible navier stokes equa tions Journal of Computational Physics 142 1 46 1998 36 William J Rider The robust formulation of approximate projection methods for incompressible flows Technical Report LA UR 3015 Los Alamos National Laboratory 1994 37 William J Rider Filtering nonsolenoidal modes in numerical solutions of incompressible flows Technical Report LA UR 3014 Los Alamos National Laboratory Los Alamos New Mexico September 1994 38 W J Rider D B Kothe S J Mosso J H Cerutti and J I Hochstein Accurate solution al gorithms for incompressible multiphase flows Technical Report AIAA 95 0699 AIAA Reno Nevada January 1995 39 William J Rider Approximate projection methods for incompressible flow implementation vari an
155. to the standard finite element formulation include the use of single point integration a row sum lumped mass matrix hourglass stabilization and balancing tensor diffusivity The benefits promised by one point integration are tremendous in computational fluid dynamics problems because of the requisite mesh sizes for interesting problems and the concomitant memory requirements The reduction from 8 quadrature points to 1 in three dimensions reduces the computational load by a factor of about 6 to 7 and reduces memory requirements by over a factor of 2 for the basic gradient operators Neglecting the storage costs associated with the PPE the total storage requirements for the explicit algorithm is 60 words per 3 D element Further it has been demonstrated that the convergence rate of the one point elements is comparable to the fully integrated elements at a fraction of the computational cost see Liu 16 With the element formulation defined attention is now turned to the parallel aspects of the explicit algorithm when a domain decomposition message passing paradigm is used 3 2 The Domain Decomposition Message Passing Paradigm This section describes the domain decomposition message passing DDMP implementation of the ex plicit time integration algorithm for the the incompressible Navier Stokes A brief overview of domain decomposition and the associated sub domain to processor mapping is discussed first Next the parallel right hand side assembly
156. tory is off default EQ 1 Convergence history is on EPS Set the convergence criteria for the iterative equation solver EQ 0 EPS 1 0e 5 default IHG Set the type of hourglass stabilization to be used with the momentum equations This only applies to the explicit treatment of the momentum equations INSOL 1 EQ 0 IHG 1 default EQ 1 LS DYNA CFD viscous hourglass stabilization EQ 2 y hourglass stabilization viscous form EHG Set the hourglass stabilization multiplier see IHG above EQ 0 EHG 1 0 default Remarks p Ze The IMASS variable is only active when INSOL gt 2 on the CONTROL_CFD_GENERAL keyword The balancing tensor diffusivity should always be used with the explicit forward Euler treatment of the advection terms This is the default The use of the flux limiting procedures is currently restricted to the explicit advective procedures DIVU sets the ceiling on the discrete divergence that is permitted during a simulation when IN SOL 1 If the divergence at a given time step exceeds the value set by DIVU then an intermediate projection is performed to return the velocity to a div free state The time weighting variables only apply to the case when INSOL gt 2 on the CON TROL_CFD_GENERAL keyword The MSOL keyword for the CONTROL_CFD MOMENTUM keyword only applies for INSOL gt 2 on the CONTROL_CFD_GENERAL keyword Incompressible Flow Solver 87 CHAPTER 9 KEYWORD INPUT
157. trategies the combination of a stabilized pressure Poisson operator with an A conjugate projection and SSOR preconditioned conjugate gradient method has been found to yield the overall best performance relative to the re solve cost of a high performance direct solver 6 1 The Saddle Point Problems The pressure Poisson equation PPE in the projection algorithms presented in Chapter 4 arises due to the use of a Helmholtz velocity decomposition that decomposes the velocity into div free and curl free components In a finite element context the Helmholtz velocity decomposition with the concomitant div free constraint CT u g is written as M C f ut d Ma 6 1 cr 0 A pe where in general C7G 4 g l The Schur complement of the projection operator in Eq 6 1 yields the PPE in Eq 4 19 In the ensuing discussion the term projection operator is used loosely to describe Eq 6 1 However the projection operator is not to be confused with the real discrete projection operators P I M C C MPEG 3 and Q I P see Gresho and Chan 4 In the explicit time integration algorithm presented in Chapter 3 a saddle point problem analogous to Eq 6 1 may also be formulated in terms of the pressure p partial acceleration and acceleration as MLC f _ f ma sy C e 41 CHAPTER 6 PRESSURE SOLUTION METHODS In this case the acceleration and pressure are used rather than velocity and pres
158. ts and robustness Technical Report LA UR 2000 Los Alamos National Laboratory Los Alamos New Mexico July 1995 40 Michael L Minion A projection method for locally refined grids Journal of Computational Physics 127 158 178 1996 41 J L Guermond and L Quartapelle Calculation of incompressible viscous flow by an uncondition ally stable projection fem Journal of Computational Physics 132 12 23 1997 42 Elbridge G Puckett Ann S Almgren John B Bell Daniel L Marcus and William J Rider A high order projection method for tracking fluid interfaces in variable density incompressible flows Journal of Computational Physics 130 269 282 1997 43 Mark Sussman Ann S Almgren John B Bell Phillip Colella Louis H Howell and Michael L Welcome An adaptive level set approach for two phase flows Journal of Computational Physics 148 81 124 1999 44 Omar M Knio Habib N Najm and Peter S Wyckoff A semi implicit numerical scheme for reacting flow ii stiff operator split formulation pre print submitted to Journal of Computational Physics 1999 Incompressible Flow Solver 143 BIBLIOGRAPHY 45 David L Brown and Michael L Minion Performance of under resolved two dimensional incom pressible flow simulations Journal of Computational Physics 122 165 183 1995 46 Michael L Minion and David L Brown Performance of under resoloved two dimensional incom pressible flow simulations ii Journal of Com
159. u MF CA 4 27 The semi implicit projection algorithm must respect a relaxed convective stability limit where CFL lt O 5 10 as described in Gresho and Chan 4 Because of the semi implicit treatment of viscous terms there is no diffusive stability constraint Computational experiments have demonstrated that CFL lt 5 is a reasonable tradeoff between accuracy and computational cost Incompressible Flow Solver 31 CHAPTER 4 THE PROJECTION METHOD In an unstructured grid with variable element size the calculation of the grid Re Reynolds and CFL Courant Freidrichs Levy numbers uses the element local coordinates and centroid velocities as de scribed in Chapter 3 The second order semi implicit projection method uses one additional modification to the finite element formulation that derives from the explicit treatment of the advective terms For advection dominated flows it is well known that the use of a backward Euler treatment of the advective terms introduces excessive diffusion Similarly Gresho et al 1 have shown that forward Euler treatment of the advective terms results in negative diffusivity or an under diffusive scheme In order to remedy this problem balancing tensor diffusivity BTD derived from a Taylor series analysis to exactly balance the diffusivity deficit is adopted In the one point quadrature element the BTD term is simply added to the kinematic viscosity in Eq 4 28 to form the tensorial dif
160. ual 27 10 1 Poiseuille Flow Poiseuille flow in a channel is characterized by a balance between a constant pressure difference and viscous shear forces 1 e 1 Ou Op Reoy dx where Re is the Reynolds number based on channel height H and u p and x y are the non dimensional x velocity pressure and coordinates The 2 D channel flow problem considered here consists of a 2 D duct with a 5 1 aspect ratio Re 100 and dp dx 0 12 This choice of Reynolds number and pressure gradient results in a steady parabolic velocity profile with Uma 1 5 and ugyg 1 0 Here the non dimensional equations were obtained using a length scale H velocity scale U and time scale U H Figure 10 1 shows the mesh with node and segment set id s for the entrance region 10 1 In order to represent the constant pressure gradient an inflow pressure boundary condition of p 0 6 is prescribed with v 0 At the inflow boundary the pressure boundary condition is prescribed using the 101 CHAPTER 10 EXAMPLE PROBLEMS Node Set ID 1 Em Segment Set ID 100 Node Set ID 1 P4 Node Set ID 3 Figure 10 1 Mesh for 2 D Entrance Region BOUNDARY PRESSURE CFD SET keyword with segment set id 100 Homogeneous outflow bound ary conditions are used at the outflow resulting in p 0 At the top and bottom wall no slip and no penetration boundary conditions u v 0 are applied using node set id 1 There are multiple flow solution algor
161. uced IntegrationOperators ee 3 2 The Domain Decomposition Message Passing Paradigm 3 2 Domain Decomposition lees iii 3 2 2 Parallel Assembly via Message Passing 3 2 3 The Parallel Explicit Algorithm 3 2 4 Communication Costs 004 4 The Projection Method 4 Semi Implicit Projection Method 4 2 Fully Implicit Projection Method ADA SAA oo Ehre fex io eR NO 422 First Time Step s AE EX 4 2 3 General Time Step 5 4 xe RICE C Ha SH AE Boundary Conditions and Source Terms 5 1 Node and Segment Sets 4 40a ww AREA ARA 5 2 Nodal Boundary Conditions 5 3 Traction Boundary Conditions 2 4 5 3 1 Pressure Boundary Conditions 5 4 Outflow Boundary Conditions ID Body FOrCES uu a uo he auc AE BE OR 5 6 Flux Boundary Conditions 2 20 4 Did PIBSSUFC LEVIS o mute y yt ctr wie ee ee wee d ele Pressure Solution Methods 6 1 The Saddle Point Problems 6 2 QlQOStabilization ll 6 3 The Projection CG Method p42 xp re ek Saree 6 4 PPE Solver Performance 4 299229 ee RS Sx RUE iv 27 29 32 32 33 33 35 35 35 37 37 38 38 39 40 41 7 Turbulence Models 7 1 Averaging Filtering and Scale Separation llle LLI Spatial PIEDRA S ES ee eee RES ee SRS SE
162. ural dynamics thermal analysis boundary element and smooth particle hydrodynamics keywords This chapter is intended to be used in conjunction with the LS DYNA Keyword User s Manual 27 which fully documents the keywords and their associated parameters Note that LS POST support for d3mean and d3thins files will only be available in LS POST 2 0 and later versions d3thins is the new time history database for the incompressible flow solver Selecting what data goes into the d3mean database is accomplished via the DATABASE_BINARY_D3MEAN keyword documented in the LS DYNA Keyword User s Manual 9 1 Incompressible Flow Keywords In LS DYNA node sets are used to prescribe nodal boundary conditions while segment sets are used to apply known flux and traction boundary conditions The use of node and segment sets follows the conventions defined in the LS DYNA Keyword User s Manual 27 75 CHAPTER 9 KEYWORD INPUT 9 1 1 BOUNDARY OUTFLOW CFD OPTION Available options include SET or SEGMENT Purpose Define passive outflow boundary conditions for the incompressible flow solvers These boundary conditions are active only when SOLN 4 or SOLN 5 on the CONTROL SOLUTION keyword In the incompressible flow solver the role of the outflow boundary conditions is to provide a computa tional boundary that is passive particularly in the presence of strong vortical flow structures Typically this boundary condition is applied at boundaries t
163. variables into grid resolved and subgrid variations 0 7 2 The Time Averaged Equations In this section the time averaged form of the conservation equations are presented As a starting point the momentum equations and the divergence constraint are Qui Ou op oO ace jn 4 2uS 7 11 P3 Tpuj Ox a ari y ij pfi and J d cs 7 12 Ox Here the viscosity is treated as a scalar for simplicity but internally LS DYNA always treats the viscosity as a second rank tensor for generality The time averaged equations are derived by applying the additive decomposition u u and applying the time averaging procedure from 7 1 where an over bar indicates application of the time averaging procedure in Eq 7 3 The time averaged momentum equations are Ou Qu op d pu 2u8 Tij i 7 13 Pog a oe ao eU ee ay where T ae uj uj Si z 3 7 14 d 2 and the Reynolds stress tensor is Tjj puju 7 15 The divergence constraint is du SEP dk 7 16 Ox The Reynolds stress tensor is a symmetric second rank tensor T 1j with six independent components Thus six independent Reynolds stress components have been added to v and p as unknowns However there are only four equations 3 momentum equations and the divergence constraint indicating that the system is not closed In fact the nonlinear terms in the momentum equation will yield new unknowns with each app
164. ver 10 2 NATURAL CONVECTION IN A SQUARE CAVITY THERMAL CAVITY RA 1000 Time 19 152 Fringe Levels Contours of X Velocity min 3 66637 at nodes 115 3 6668 00 max 3 66637 at node 227 2 933e 00 2 200e 00 1 467e 00 7 333e 01 1 192e 07 7 333e 01 1 467e 00 2 200e 00 gt 2 933e400 MA of 3 666e 00 Figure 10 14 Thermal cavity x velocity contours THERMAL CAVITY RA 1000 Time 19 152 Fringe Levels Contours of Y Velocity min 3 67413 at nodes 297 3 674e 00 max 3 67413 at node 55 2 939e 00 2 204e 00 1 470e 00 7 348e 01 uM a 1 192e 07 3 7 348e 01 v i ee 1 470e 00 RR 2 204e 00 2 939e 00 P ades 3 674e 00 E mr d Pi 2 F pn d Y lt x Figure 10 15 Thermal cavity y velocity contours Incompressible Flow Solver 119 CHAPTER 10 EXAMPLE PROBLEMS THERMAL CAVITY RA 1000 52 Time 1941 Fringe Levels pris asa re 121 4 026e 01 max 40 2799 at node 432 3 308e 01 E E A tj 2 588e 01 1 868e 01 1 147e 01 4 272e 00 z 2 929e 00 N 1 013e 01 li 1 733e 01 Y 2 453e 01 3 74e401 Pd A i gt dl Y Le A IEEE a Figure 10 16 Thermal cavity z vorticity contours 120 Incompressible Flow Solver 10 3 THE MOMENTUM DRIVEN JET 10 3 The Momentum Driven Jet This computation shows the starting vortical structure associated with a heated slot jet ent
165. wrt ihist 0 0 iwrt ihist 1 0 endeng endmas thetaa thetaf eps ihg ehg 1 13 0 eps 1 0e 05 5 PAPA offo dattyp offo dattyp Define the load curves DEFINE CURVE lcid sidr sfa sfo offa 1 0 0 1 0 9999 0 1 0 DEFINE CURVE lcid sidr sfa sfo offa 2 0 0 0 0 9999 0 0 0 SSS SSS SSS SS SSS SSS SS SS 9999999 9 999 9 SSS SSS SSS SSS SSSSSSSSSSSSSSSSSSSSSSSSSSsss Setup the boundary conditions nsid 1 no slip no penetration surface nsid 2 on axis inlet nsid 3 top bottom inlet ssid 4 no slip no penetration surface ssid 5 on axis inlet surface ssid 6 top bottom surfaces Sheath chamber walls no slip no penetration BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 101 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 102 2 1 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 1 103 2 1 0000000 x inlet surface BOUNDARY PRESCRIBED CFD SET Incompressible Flow Solver 133 CHAPTER 10 EXAMPLE PROBLEMS nsid dof lcid sf 2 102 0 2 0 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 2 103 0 2 0 0000000 top bottom x z velocity BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 3 101 0 2 0 0000000 BOUNDARY PRESCRIBED CFD SET nsid dof lcid sf 3 103 0 2 0 0000000 Prescribed inlet pressure BOUNDARY PRESSURE CFD SET ssid lcid p 2 1 11000 0 BOUNDARY PRESSURE CFD SET ssid lcid p 3 1 10000 0 59 9
166. y for node 289 at x y 1 34e 2 3 71e 2 Time Incompressible Flow Solver 127 CHAPTER 10 EXAMPLE PROBLEMS MOMENTUM DRIVEN JET f Time 0 7758 0 O Fringe Levels nin 1 20864 at mocos 547 1 313800 max 1 3131 at node 303 1 061e 00 8 088e 01 0 9 5 566e 01 3 3 3 044e 01 I d 5 223e 02 1 999e 01 8 4 521e 01 A 7 043e 01 E 4 9 565e 01 EN 1 209e 00 QO QO OW QO 89 Y Ke o Figure 10 21 X velocity snapshot at t 0 7758 s MOMENTUM DRIVEN JET V 7 Time 0 7758 i ringe Levels pary Tre t nodes 436 5 133800 max 5 13854 at node 3722 4 553e 00 3 967e 00 3 382e 00 IN 2 796e 00 2 210e 00 1 625e 00 V 1 039e 00 4 534e 01 1 322e 01 7 179e 01 Y i Figure 10 22 Y velocity snapshot at t 0 7758 s 128 Incompressible Flow Solver MOMENTUM DRIVEN JET Time 0 7758 Contours of Temperature min 300 at node 344 max 400 019 at node 201 5 x z a ANS pore m Sis 10 3 THE MOMENTUM DRIVEN JET Fringe Levels 4 000e 02 3 900e 02 3 800e 02 3 700e 02 3 600e 02 3 500e 02 3 400e 02 3 300e 02 3 200e 02 3 100e 02 3 000e 02 Figure 10 23 Temperature snapshot at t 0 7758 s radi MOMENTUM DRIVEN JET Time 0 7758 Contours Pressure min 5 023 max 0 947254 at node 10652 Fringe Pa 9 473e 01 3 502e 01 2 468e 01 8 439

Download Pdf Manuals

Related Search

Related Contents

T910 User Manual pdf  exhibitor service kit - nj - Florida Council of Independent Schools  La finance ser encore l`écono La finance sert  StarTech.com SNPMHDHD  DESCRIPCION Ambientador elaborado a base de  WARNING Kebo Ottoman    Samsung SGH-P300 manual de utilizador  NANO 2  Manual do Usuário  

Copyright © All rights reserved. Failed to retrieve file