VIEWS: 20 PAGES: 18 CATEGORY: Childrens Literature POSTED ON: 11/28/2009
Three-dimensional-free-surface-modelling-in-an-unstructured
Ninth International PHOENICS User Conference September 23 – 27 2002, Moscow, Russia Two-dimensional free surface modelling for a non-dimensional Dam-Break problem M. Ben Haj(1) , Z. Hafsia(2) , H. Chaker(3) and K. Maalel(4) (1) and (3) : Ecole Nationale d’Ingénieurs de Tunis. Laboratoire de Modélisation Mathématique et Numérique dans les Sciences de l’Ingénieur. B.P. 37 – Le Belvédère, 1002, Tunis, Tunisie. (2) and (4) : Ecole Nationale d’Ingénieurs de Tunis. Laboratoire de Modélisation en Hydraulique et Environnement. B.P. 37 – Le Belvédère, 1002, Tunis, Tunisie. E-Mails: (1) (3) : Mehdi1@webmails.com : Hedia.Chaker@enit.rnu.tn (2) (4) : Zouhaier.Hafsia@enit.rnu.tn : Khlifa.Maalel@enit.rnu.tn Date: 12/29/2001 Hardware: Sun Sparc Station 5 PHOENICS version: 2.10 Abstract This work presents numerical computations for the analysis of Dam-Break Flow using two-dimensional equations in a vertical plane. The numerical results are obtained using the SEM and HOL methods, implemented in the PHOENICS code, to predict the time variation of the free surface. In this code the finite volume method is used to discretize Navier-Stokes and transport equations. The time evolution of the flow depth and the pressure distribution at the dam site are investigated for dry-bed condition. The effect of the initially non-hydrostatic state on the long term surface profile and wave velocity are studied. The simulated results showed that this effect is significant. The simulated dry-bed tip velocity after the Dam-Break, compared with the analytical results. Contents list 1. Introduction 2. Mathematical description 2.1. The fluid flow equations 2.2. The free surface model 2.2.1 The Scalar Equation Method (SEM) 2.2.2 The Height of Liquid Method (HOL) 3. Results and discussion 3.1 Free surface profile 3.2 Propagation of the tip 3.3 Evolution of Flow depth at dam site 3.4 Evolution of pressure distribution 4. Conclusions 1 1. Introduction The analysis of Dam-Break Flow (DBF) is a part of dam design and safety criteria and it is also crucial for flood plain management and before locating infrastructures in river valleys (Stoker 1952 and Sentürk 1994). The basic formulation of DBF mechanics date back to the earliest attempt by Ritter (1892). He derived an analytical solution for the hydrodynamic problem of instantaneous Dam-Break in a frictionless and horizontal channel of rectangular shape. It has been well established that the pressure distribution is non-hydrostatic immediately after the dam failure (Mohapatra and al. 1999). Hence, the DBF mechanics must be analyzed regarding the hydrostatic pressure assumption. This work concerns a detailed numerical study of DBF for dry-bed condition using two methods implemented in PHOENICS code (SEM: Scalar Equation Method and HOL: Height Of Liquid Method) for tracking the free surface. The turbulence effects are not considered to describe DBF (the turbulence would cause diffusion effects, which are minor when compared with the strong convective nature of the flow). Various studies concerning numerical simulations of DBF using two-dimensional equations have been reported: Harlow and Welch 1965; Amsden and Harlow 1970; Hirt and Nichols 1981; Jha, Akiyama and Ura 1995. In all these studies DBF was used only as a test case for free surface tracking model. In this paper, we first present a brief description of the SEM and HOL methods used to predict the free surface evolution. The numerical method used for this purpose is a composite one based using the finite volume method to discretize Navier-Stokes and transport equations and the Simplest algorithm (Spalding 1981) for the treatment of pressure-velocity coupling. Mechanics of DBF for dry-bed condition are analyzed based on numerical simulations. The time evolution of the flow depth at the dam site and the evolution of the pressure distribution are compared relatively to a hydrostatic state. The effect of the initially non-hydrostatic state on the long-term surface profile and the wave velocity are also studied. The numerical results for the wave speed and the time interval of the initial nonhydrostatic state are also compared with the analytical results. 2. Mathematical description The two fluids separated by the interface are treated, as both belonging to a single-phase, i.e. there is only one value of each velocity component, density, dynamic viscosity, etc., for each computational cell. The relevant governing equations need to be solved, accounting for the discontinuities in the physical properties of fluids at the interfaces. The mathematical model will need to: Locate the unknown inter-fluid boundaries; Satisfy the field equations governing conservation of mass, momentum; Be consistent with the boundary conditions. 2.1. The fluid flow equations The general equations that govern two-dimensional transient fluid flow are given by the momentum and mass continuity equations which are written in their vector from (Pericleous and al. 1998): ρu . (ρu u ) . ( μu ) S t (1) 2 ρ . (ρu ) 0 t (2) where u is the velocity vector. The other variables at time t are ρ and μ , the single-phase density and dynamic viscosity. The source term S in Eq. (1) represents the gravitational force and pressure. These equations could be written as one equation: ρ . (ρ u ) . ( μ) S (3) t For 1 , we find the continuity equation (2) and for u , we find the momentum equation (1). Time integration of Eq. (3) over a control volume in the y-z plane gives (figure 1): t ρ . (ρ u ) . ( μ) S dVdt 0 t VP (4) A cell containing node P has four neighboring nodes identified as South, North, Low and High nodes (S, N, L, H). The notation, s, n, l and h designed the south, north, low and high cell faces respectively. High (H) h P = Center of the cell n North (N) South (S) s P z l y Low (L) Figure 1: Control volume. Eq (4) gives in discrete and implicit formulation (Rosten and Spalding 1987), P aN N aSS aHH aLL aTT S aN aS aH aL aT aP (5) 3 The values of in the neighboring cells N, S, H and L are, by default, those which prevail at the end of the current time step; T denotes the value of in the P cell at the end of the previous time step. For each scalar variable (ie. other than velocities), there is one such equation for each of the NY*NZ cells in the domain of integration and the values of these variables are given in the center of each cell. The velocity components are given on the faces of each cell. For scalar variables, the coefficients expressions which appeared in the equation (5) are: aN max( 0 , dn α mn ) max ( 0 , mn ) aS max( 0 , ds α ms ) max ( 0 , m ) s aH max( 0 , dh α mh ) max ( 0 , mh ) aL max( 0 , dl α ml ) max ( 0 , m ) l In the above expressions, d denotes the diffusive coefficient, and m denotes the convective mass flux across cell faces. The parameter α is equivalent to the parameter DIFCUT in Q1 file, which governs the relative contribution of convection and diffusion. Its value is set by the PIL variable DIFCUT which has a default value of ½, and which ensures that the contribution of diffusion is omitted when m/ d 2 . This is commonly known as the “hybrid” scheme. The expressions used for the coefficients of the v-velocity equation are as follows: v v v aN max( 0 , dN α mN ) max ( 0 , mN ) v v aS max( 0 , dP α mP ) max ( 0 , mv ) P v v v aH max( 0 , dhn α mhn ) max ( 0 , mhn ) v v aL max( 0 , dln α mln ) max ( 0 , mv ) ln The expressions used for the coefficients of the w-velocity equation are as follows: w w w aN max( 0 , dhn α mhn ) max ( 0 , mhn ) w w aS max( 0 , dhs α mhs ) max ( 0 , mw ) hs w w aS max( 0 , dhs α mhs ) max ( 0 , mw ) hs w w w aH max( 0 , dH α mH ) max ( 0 , mH ) w w aL max( 0 , dP α mP ) max ( 0 , mw ) P The solution procedure is based upon the iterative SIMPLEST algorithm (Spalding 1981), which is a modified version of SIMPLE algorithm (Patankar 1980 ; Versteeg and Malalasekera 1995) and which ensures more rapid convergence. 2.2. The free surface model The Scalar Equation (SEM) and Height of Liquid (HOL) are two methods incorporated in PHOENICS code both essentially single-phase treatments (The PHOENICS Encyclopedia). These methods employ 4 a scalar to represent the liquid volume fraction in a control volume cell. In a gas cell, 0 , in a liquid cell, 1, and for a partially filled cell takes on an intermediate value between 0 and 1. The value of the scalar is used to update the fluid properties. The density and the viscosity fields can be calculated by the following expressions: ρ ρG(1 ) ρL () (6) (7) G(1 ) L () Where ρ G , ρ L are respectively the gas and the liquid densities and μ G , μ L are respectively the gas and the liquid dynamic viscosities. At the interface, the properties ρ and μ are assumed to vary linearly with from those of one fluid to those of the other. 2.2.1 The Scalar Equation Method (SEM) A- The Governing Equation The dynamics of the free surface in The Scalar Equation Method (SEM) are governed by the transport equation of the liquid volume fraction denoted by : On the assumption that the transport equation of is governed only by convection. B- Van Leer discretisation of the scalar-convection terms Using Green formula, time integration of Eq. (8) over the control volume gives: Φ . (Φ.u ) 0 t (8) Φ dVdt Φu ds 0 t i tV s (9) which is written in discrete form: n n ΦP 1 ΦP Faces n 1 i i 1 u . ni Ai Δt n Φi(U) VP (10) where the subscript i sums over the faces of the control volume P . The second term on the right hand n side of Eq. (10) is the flux summation term. The value Φi(U) is the value of Φ at the face, which, in the PHOENICS code, is taken to be the upwind cell value denoted by the subscript U. The velocity n 1 ui is the latest value of the face velocity and n is the surface normal of the face i . The "fully implicit upwind" scheme considers the transport of a variable, Φ , across a cell face (e.g. the north face) according to the following formula (Norris 1996): Φn= ΦP Φn= ΦN for Φ n > 0 for Φ n < 0 (11) where Φ P and Φ N are the values of Φ at the grid nodes at the end of the time step. 5 The term ui . ni Ai ΔtVP1 in the flux summation term of Eq. (10), is the Courant-Friedrichs-Lewy (CFL) stability number, n 1 i ui n 1 . ni Ai Δt VP (12) therefore, Eq. (10) can be rewritten in terms of i , n n ΦP 1 ΦP Faces i 1 i n Φi(U) (13) When first order differencing schemes (hybrid, upwind, etc.) are used to solve the convection Eq. (8), the interface smears and becomes less distinct (under- or overshoots). Various numerical methods can be used to overcome this numerical smearing known as numerical diffusion. In this study, the Van Leer TVD scheme is adopted, as being the simplest to implement and most economical of the higher order schemes available (Pericleous and al. 1995). The Van Leer approach is an explicit finitedifference scheme which modifies the first order upwind formulation using the theory of characteristics: Φ n = Φ P + (d Φ /dy)P . (dy - undt )/2 Φ n = Φ N - (d Φ /dy)N . (dy + undt )/2 for for Φn>0 Φn<0 (14) where Φ P and Φ N are the values of Φ at the end of the time step. The d Φ /dy is based on values at the start of the time step, which implies that the scheme is explicit in nature. This gradient is specified such that the scheme reduces to the purely upwind form in the presence of local extrema. This Total Variation Diminishing (TVD) approach offers the advantages of a higher-order scheme while avoiding under- or overshoots. Due to the explicit formulation, the Courant criterion places a maximum limit on the time increment for the stability of the solution: dt = min (dy/v, dz/w ) where this minimum is with respect to every cell in the grid, not just at the interface. 2.2.2 The Height of Liquid Method (HOL) In this method the liquid volume fraction denoted by is defined like the average between the liquid volume in the cell and the total volume of the cell: Φ VL VT (15) Where VL and VT are respectively the liquid volume in a cell and the total volume of this cell. 6 For this purpose, the domain is decomposed to many columns in the vertical direction and for each column the variation of the liquid quantity with time is determined. This variation is calculated from a mass balance equation: Φ MT M L L VT (16) Where, ρ L : liquid density; ρ G : gas density; : volume fraction of liquid; VT : total cell volume; MT : total mass of liquid in a vertical column of cells. For all subsequent times, MT is determined from: MT = MT0 + ΔF MT0 : total mass of liquid in a vertical column at the old time step. ΔF : sum of inflows of liquid to the cell column minus the sum of outflows of liquid from the cell column. The inflows and outflows are computed from the volumetric flows obtained as convective fluxes ( ui . Ai . t ), multiplied by the upwind values of ( . ρ L ). The sign of ΔF can be positive or negative. The contribution of the flux in the column direction is not considered in ΔF (only the fluxes in the other directions are calculated). F Ai t i L u i n i Ai t i L ui n i (17) in out Where u i : is the value of an inflow or outflow velocity of a cell face from the column; A i : is the surface of a face; n i : is the surface normal of the face i and t : is the time interval in which the flow is effectuated. ML : is the mass of liquid in a column of full cells below the free surface (figure 2): M L (L VT) E MT L VT E : is the floor function. (18) Finally, the position of the interface is defined by calculating the volume fraction of liquid, in every cell, from equation (16). The interface is defined to lie in the cells with a in the range 0.0 < < 1.0 . 7 7 z 6 5 MT - M L 4 3 2 1 ML Figure 2: Definition of variables. When ML is greater than MT, equation (16) produces a negative , which is of course unrealistic, and should be increased to zero. Therefore we take = max ( , 0.0). Similarly, when MT is greater than ML, could be greater than one, which should be reduced to one. Therefore we take = min ( , 1.0). 3. Results and discussion In our study, an instantaneously dam break flow is analyzed. The location of the dam site is at y = 0 m. The location of the inflow section is –300 m and the location of the outflow section is +300 m (Figure 3). The initial upstream flow depth in the reservoir before the dam failure is h 1 = +10 m and the initial downstream flow depth is equal to zero (dry downstream bed). At t = 0, the velocity in inflow zone is U1 = 0. The channel bed slope is equal to zero in all the computations. The initial conditions of zero velocity and known depths at the inflow and the outflow boundaries are kept unchanged as the boundary conditions during the transient computations. The horizontal cells number upstream the dam site is NY1 = 60 for SEM and NY1 = 300 for HOL and downstream the dam site is NY2 = 60 for SEM and NY2 = 300 for HOL. The vertical cells number are NZ1 = 20 for both SEM and HOL methods. The computations are performed for a time of 15 s and with a time step ∆t = 0,2 s for SEM and ∆t = 0,04 s for HOL. The Q1 files for the Dam-Break Flow with the two methods SEM and HOL are reported respectively in the appendices 1 and 2. 8 10 m U1 = 0 y=0 300 m 300 m Figure 3: The problem position. All computations in this study are presented using non-dimensional variables in which the references scales adapted are: L0 = 1 m : Length scale, g0 = 9,8 m/s2 : Gravity scale, V0 = 3,13 m/s : Velocity 0 scale, P0 = 9,8 pa : Pressure scale, T0 = 0,319 s : Time scale. L0 h1 , g0 g , V0 L0.g0 , P l.V02 et T0 L0 . L0 is the upwind water depth, g 0 : the gravity, V0 : Velocity scale, P0 : Pressure scale and T0 : V0 time scale. Note that the wave does not arrive at the outflow boundary for the simulation time considered. 3.1 Free surface profile In this section free surface profiles are presented. We compare the numerical results for long-term free surface profiles with the analytical solutions obtained using the hydrostatic pressure assumption. Figures 4 a) and b) shows the numerical results for the long-term free surface profiles at t* = 15 using respectively the SEM and HOL superposed with the non-dimensional analytical solution (Wu and al. 1999): h 1 h1 if if y* 1 1 y* 2 y* 2 (19) h 1 2 y* 2 ) h1 9 ( h 0 h1 Where y* if y and h1 is the initial upstream flow depth in the reservoir. t gh1 Figure 4, showed the non-dimensional flow depth Z as a function of the non-dimensional distance from the dam Y. The location of the tip in the cases of SEM [Fig. 4 (a)] and HOL [Fig. 4 (b)] is under predicted by the analytical model, as compared with the numerical result. The difference between the two is about 6% and 13% in case of SEM and HOL respectively. The height and the location of the positive wave front is more approached using SEM. 9 a) b) Figure 4: Non-dimensional Free Surface Profiles a) For SEM method b) For HOL method. 3.2 Propagation of the tip Figures 5 a) and b) present the non-dimensional front velocity for Dam-break Flow given by SEM and HOL methods. In these figures, t represents the non-dimensional time after the collapse of the dam, and Yw the non-dimensional distance traveled by the tip. The slopes of these curves give the velocities of the tip in both cases. It can be clearly seen from Figures 5 a) and b) that the two dimensional effects reduce the rate at which the tip advances on a dry bed compared to the analytical solution. It is found from Fig. 5 a) that, on average, the tip advances with a non-dimensional speed of 1,75 for t<t* = 3 and 1,55 for t>t* = 3 rather than a value of 2 as suggested by Ritter (1892). From Fig. 5 b) it can be seen that, the tip advances with a non-dimensional speed of 1,06 for all times which is smaller than the value 2. These results indicate a significant long-term effect of non-hydrostatic pressure distribution, in the case of dry-bed conditions. a) b) Figure 5: Non-dimensional Front Location a) For SEM method b) For HOL method. 3.3 Evolution of Flow depth at dam site Ritter’s (1982) solution, which use the hydrostatic assumption, predict that the flow depth at the dam site attains a constant value of 4/9 instantaneously upon the dam break (Graf 1996). From equation (19) and for y* 0 , we obtain that h/h1 4/9 . However, with the SEM and HOL methods, the flow depth at the dam site takes some times to attain this constant value which are close to the numerically value of eight non-dimensional time units. The 10 variation with time, t, of the flow depth at the dam site, h, obtained using the SEM and HOL methods is shown in Fig. 6. a) b) Figure 6: Time Variation of Flow Depth at Dam Site a) For SEM method b) For HOL method. Figures. 6 (a) and 6 (b) show, that the flow depth at the dam site, h, attain the constant value of 4/9 for six and seven non-dimensional time units with SEM and HOL respectively. The simulated fixed point localized at the dam site was reproduced by using the two methods SEM and HOL. 3.4 Evolution of pressure distribution Figures 7 (a) and (b) show the pressure history on the floor at the dam site for both SEM and HOL. In these figures, the non-dimensional pressure is plotted as a function of the non-dimensional time. In both cases the pressure is not equal to the hydrostatic pressure at the beginning due to the streamline curvature. It eventually approaches the hydrostatic value as time progress. a) b) Figure 7: Pressure History at Dam Site a) For SEM method b) For HOL method. The pressure distribution at non-dimensional times t = 1, 5, and 10 at the dam site for both methods SEM and HOL is respectively presented in figures 8 and 9. It is clear from these figures that the pressure distribution at the dam site is not equal to the hydrostatic value obtained immediately after the dam failure but it’s greater than the hydrostatic value at the bottom of the dam site and the pressure distribution approaches the hydrostatic value as time progress. Figures 8 and 9 (for SEM and HOL) showed that the flow depth at the dam site attains a constant value before five non-dimensional time units. 11 Figure 8: Evolution of Pressure Distribution at Dam Site for SEM method. Figure 9: Evolution of Pressure Distribution at Dam Site for HOL method. 4. Conclusions Analysis of the Dam-Break Flow mechanics for dry-bed condition is presented is based on the numerical prediction of free surface evolution obtained using two methods (SEM and HOL) in the PHOENICS code. A comparative study of simulated results in term of non-dimensional variables with the analytical solution allowed the following conclusions: 12 1. The location of the tip in the cases of SEM and HOL is under predicted by the analytical model, as compared with the numerical result. 2. The two dimensional effects reduce the rate at which the tip advances on a dry bed for SEM and HOL which is smaller than a value of 2 as suggested by Ritter 1892. These results indicate a significant long-term effect of non-hydrostatic pressure distribution, in the case of dry-bed condition. 3. Ritter’s (1982) solution, which use the hydrostatic assumption, predict that the flow depth at the dam site attains a constant value of 4/9 instantaneously upon the dam break. However, with the SEM and HOL methods, the flow depth at the dam site takes some times to attain this constant value. 4. In both cases of SEM and HOL, the pressure is not equal but greater than the hydrostatic pressure at the beginning due to the streamline curvature. It eventually approaches the hydrostatic value as time progress. References [1] A.A. AMSDEN and F.H. HARLOW. (1970). "A Simplified MAC Technique for Incompressible Fluid Flow Calculations." Journal of Computational Physics 6, pp 322-325. [2] CHAM-Concentration Heat and Momentum. (2001). "The PHOENICS Encyclopedia." [3] W.H. GRAF. (1996). "Hydraulique fluviale: Ecoulement non permanent et phénomènes de transport." Tome 2, presses polytechniques et universitaires romandes. [4] F. H. HARLOW and J. E. WELCH. (1965). "Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface" The Physics of Fluids, vol. 8(12), pp: 2182-2189. [5] C. W. HIRT, B. D. NICHOLS. (1981). "Volume-Of-Fluid (VOF) method for the dynamics of free Boundaries." Journal of Computational Physics, 39, 201-225. [6] A. K. JHA, J. AKIYAMA & M. URA. (1995). " First- and Second-Order Flux Difference Splitting Schemes for Dam-Break Problem." Journal of Hydraulic Engineering, Vol. 121, No. 12, pp. 877-884. [7] P.K. MOHAPATRA, V. ESWARAN and S. MURTY BHALLAMUDI. (1999). " Two-dimensional analysis of dam-break flow in vertical plane." Journal of Hydraulic Engineering, Vol. 125, N. 2. [8] S. E. NORRIS. (1996). "A comparison of differencing schemes for time marching finite volume methods." School of Mechanical Engineering, University of Sydney, Sydney, Australia. [9] S.V.PATANKAR. (1980). "Numerical Heat Transfer and Fluid Flow." Series in Computational Methods in Mechanics and Thermal Sciences, Minkowycz and Sparrow, Editors. [10] K.A. PERICLEOUS, K.S. CHAN and M. CROSS. (1995). "Free surface flow and heat transfer in cavities: the sea algorithm." Numerical Heat Transfer, Part B, 27, pp: 487-507. [11] K.A. PERICLEOUS, G.J. MORAN, S.M. BOUNDS, P. CHOW AND M. CROSS. (1998). "Three-dimensional free surface modelling in an unstructured mesh environment for metal processing applications." Applied Mathematical Modelling 22, pp: 895-906. 13 [12] A. RITTER. (1892). "The propagation of the water waves" V.D.I Zeitschrift des Vereines Deutscher Ingenieure, v.36, pp: 947-954. http://www-il.usgs.gov/proj/feq/fequtl/fequtl.chap4_7.html [13] H. I. ROSTEN and D. B. SPALDING. (1987). "The PHOENICS EQUATIONS." CHAM TR/99, Document Revision 02, Software Version 1.4 [14] F. Sentürk (1994). "Hydraulics of dams and reservoirs." Water Resources Publications, Highlands Ranch, Colorado, USA. [15] D. B. SPALDING. (1981). "A general purpose computer program for multi-dimensional one- and two-phase flow." Mathematics and Computers in simulation XXIII (1981) 267-276, North-Holland Publishing Company. [16] J. J. STOKER. (1952). "Water waves." Wiley Classics Library Edition Published 1992. [17] H. K. VERSTEEG and W. MALALASEKERA. (1995). "An introduction to computational fluid dynamics: The Finite Volume Method." Longman Group Ltd, Malaysia. [18] C. WU, G. HUANG and Y.ZHENG (1999). "Theoretical solution of Dam-Break shock wave." Journal of Hydraulic Engineering, Vol. 125, No. 11, pp. 1210-1215. Appendices 1 SEM (Q1 file) TALK=T;RUN( 1, 1);VDU=X11-TERM GROUP 1. Run title and other preliminaries TEXT(Dam-Break by SEM) INTEGER(NY1,NY2,NZ1,NZ2) NY1=60;NY2=60;NZ1=20;NZ2=1 NY=NY1+NY2 NZ=NZ1+NZ2 GROUP 2. Transience; time-step specification STEADY=F; LSTEP=75 ;TLAST=15 GRDPWR(T,LSTEP,TLAST,1.0) GROUP 4. Y-direction grid specification NREGY=2 IREGY=1;GRDPWR(Y,NY1,300,1.0); IREGY=2;GRDPWR(Y,NY2,300,1.0); GROUP 5. Z-direction grid specification NREGZ=2; IREGZ=1;GRDPWR(Z,NZ1,10,1.0); IREGZ=2;GRDPWR(Z,NZ2,0.1,1.0); GROUP 7. Variables stored, solved & named STORE(DEN1,PRPS); SOLVE(VFOL,SURN) SOLUTN(P1,Y,Y,Y,N,N,N) SOLUTN(V1,Y,Y,N,N,N,N) SOLUTN(W1,Y,Y,N,N,N,N) GROUP 8. Terms (in differential equations) & devices 14 ** activate the "gas-and-liquid algorithm", ie volumetric continuity equation, and allow convection fluxes to be modified in GROUND GALA=T; TERMS(VFOL,N,N,N,N,P,P) TERMS(SURN,N,N,N,N,P,P) GROUP 9. Properties of the medium (or media) ** signal that density is to be computed by the HOL method and set the densities of the liquid and gas respectively RHO1=GRND10 ENUL=GRND10 GROUP 11. Initialization of variable or porosity fields FIINIT(P1)=0.0;FIINIT(V1)=0.0;FIINIT(W1)=0.0; FIINIT(SURN)=0.0; FIINIT(PRPS)=0. FIINIT(DEN1)=1.189; INIADD=F PATCH(LIQUID,INIVAL,1,NX,1,NY1-1/3,1,NZ1,1,1) init(LIQUID,SURN,ZERO,1) init(LIQUID,VFOL,ZERO,1) init(LIQUID,PRPS,ZERO,67) init(LIQUID,DEN1,ZERO,1000.5) GROUP 13. Boundary conditions and special sources ** the pressure is held to zero along the open top boundary PATCH(REFP,CELL,1,NX,1,NY,NZ,NZ,1,LSTEP) COVAL(REFP,P1,FIXVAL,ZERO) ** provide for the gravity-force source of w1 PATCH(GRAV,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(GRAV,W1,FIXFLU,-9.81) GROUP 15. Termination of sweeps LSWEEP=50; LITER(SURN)=1 GROUP 16. Termination of iterations SELREF=T;RESFAC=0.01 GROUP 17. Under-relaxation devices RELAX(P1 ,LINRLX,0.7) RELAX(V1,FALSDT,1) RELAX(W1,FALSDT,1) GROUP 19. Data communicated by satellite to GROUND ** provide for the dumping of field data at each time step, for us by PHOTON IDISPA=5;IDISPB=1;IDISPC=LSTEP;CSG1=M IPRPSA=67;IPRPSB=0; SURF=T; RLOLIM=0.4;RUPLIM=0.6 VARMIN(SURN)=0.0; VARMAX(SURN)=1.0 15 GROUP 22. Spot-value print-out TSTSWP=-1 ECHO=T IYMON=2*NY/3;IZMON=NZ/2 GROUP 23. Field print-out and plot control OUTPUT(P1,N,N,N,N,Y,Y);OUTPUT(V1,N,Y,N,Y,Y,Y) OUTPUT(W1,N,N,N,N,Y,Y);OUTPUT(SURN,N,N,N,N,N,N) OUTPUT(DEN1,N,N,N,N,N,N) OUTPUT(VFOL,Y,N,N,N,Y,Y) OUTPUT(PRPS,N,N,N,N,N,N) NTPRIN=75 NYPRIN=1 NZPRIN=1 IPROF=2 XZPR=T STOP Appendices 2 HOL (Q1 file) TALK=T;RUN( 1, 1);VDU=X11-TERM GROUP 1. Run title and other preliminaries TEXT(Rupture de barrage par HOL) INTEGER(NY1,NY2,NZ1,NZ2) NY1=300;NY2=300;NZ1=20;NZ2=1 NY=NY1+NY2 NZ=NZ1+NZ2 GROUP 2. Transience; time-step specification STEADY=F; LSTEP=375 ;TLAST=15 GRDPWR(T,LSTEP,TLAST,1.0) GROUP 4. Y-direction grid specification NREGY=2 IREGY=1;GRDPWR(Y,NY1,300,1.0); IREGY=2;GRDPWR(Y,NY2,300,1.0); GROUP 5. Z-direction grid specification NREGZ=2; IREGZ=1;GRDPWR(Z,NZ1,10,1.0); IREGZ=2;GRDPWR(Z,NZ2,0.01,1.0); GROUP 7. Variables stored, solved & named STORE(DEN1,PRPS,RHO1,ENUL,IMB1) SOLVE(VFOL) SOLUTN(P1,Y,Y,Y,N,N,N) SOLUTN(V1,Y,Y,N,N,N,N) SOLUTN(W1,Y,Y,N,N,N,N) GROUP 8. Terms (in differential equations) & devices ** activate the "gas-and-liquid algorithm", ie volumetric continuity equation, and allow convection fluxes to be modified in GROUND * Y in TERMS argument list denotes: * 1-built-in source 2-convection 3-diffusion 4-transient 16 * 5-first phase variable 6-interphase transport TERMS (P1 ,Y,Y,Y,N,Y,N) TERMS (V1 ,Y,Y,N,Y,Y,N) TERMS (W1 ,Y,Y,N,Y,Y,N) DIFCUT = 5.000E-01 ; GALA=T; TERMS(VFOL,N,N,N,N,P,P) GROUP 9. Properties of the medium (or media) ** signal that density is to be computed by the HOL method and set the densities of the liquid and gas respectively RHO1=GRND10;ENUL=GRND10 GROUP 11. Initialization of variable or porosity fields FIINIT(P1)=0.0;FIINIT(V1)=0.0;FIINIT(W1)=0.0; FIINIT(VFOL)=0.0 FIINIT(DEN1)=1.189 INIADD=F ** place the initial "pile" of liquid of which the slumping is to be simulated ** redefine patch for reduced grid PATCH(LIQUID,INIVAL,1,NX,1,NY1,1,NZ1,1,1) init(LIQUID,VFOL,ZERO,1.0) init(LIQUID,DEN1,ZERO,1000.5) GROUP 13. Boundary conditions and special sources ** the pressure is held to zero along the open top boundary PATCH(REFP,CELL,1,NX,1,NY,NZ,NZ,1,LSTEP) COVAL(REFP,P1,FIXVAL,ZERO) ** provide for the gravity-force source of w1 PATCH(GRAV,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP) COVAL(GRAV,W1,FIXFLU,-9.81) GROUP 15. Termination of sweeps LSWEEP=50 GROUP 16. Termination of iterations SELREF=T;RESFAC=0.01 GROUP 17. Under-relaxation devices RELAX(P1 ,LINRLX,0.7) RELAX(V1,FALSDT,1.0); RELAX(W1,FALSDT,1.0) GROUP 19. Data communicated by satellite to GROUND ** provide for the dumping of field data at each time step, for us by PHOTON HOL=T IHOLA=3 IDISPA=LSTEP;IDISPB=1;IDISPC=LSTEP;CSG1=F IPRPSA=67;IPRPSB=0 GROUP 22. Spot-value print-out TSTSWP=-1;ECHO=T; IYMON=2*NY/3;IZMON=NZ/2 GROUP 23. Field print-out and plot control 17 OUTPUT(P1,N,N,N,N,Y,Y);OUTPUT(V1,N,Y,N,Y,Y,Y) OUTPUT(DEN1,N,N,N,N,N,N) OUTPUT(VFOL,Y,N,N,N,Y,Y) OUTPUT(PRPS,N,N,N,N,N,N) NTPRIN=75 NYPRIN=1 NZPRIN=1 IPROF=2 XZPR=T STOP 18