Search **web-archive-uk.com**:

Find domain in archive system:

Find domain in archive system:

web-archive-uk.com » UK » C » CHAM.CO.UK Total: 682 Choose link from "Titles, links and description words view": Or switch to
"Titles and links view". |

- CONVERGENCE MONITORING AND CONTROL

equations generates a field of f values f new Linear relaxation replaces this with f new f old a f new f old 5 where f old is the value from the previous sweep a between 0 and 1 is the relaxation coefficient A relaxation coefficient of 0 prevents any change from the previous sweep a value of 1 applies no relaxation False timestep relaxation False timestep relaxation modifies the finite volume equations by adding an additional pseudo transient term mass in cell f old f new dt f where dt f is the false timestep A large value of dt f makes the additional term small light loose relaxation a small value makes the additional term large heavy tight relaxation Choice of relaxation Only linear relaxation can be applied to pressure and relaxation is usually only needed in BFC cases False timestep is usually applied to velocities The value of dt f is usually based on a characteristic time scale for the problem this is often a cell residence time cell dimension velocity but other choices are sometimes more appropriate buoyancy based viscosity based etc depending on which effects are dominant in the flow Experience and trial and error are the only real guides Slowly varying spot values may indicate excessively tight relaxation over oscillatory spot values may indicate insufficiently heavy relaxation Relaxation can often be lightened after the early stages of a calculation Linear relaxation can be applied to stored variables such as density this can sometimes help to control divergence Automatic relaxation CONWIZ CONWIZ is the current default relaxation mechanism for cases set through the VR Editor The main features of CONWIZ are It starts by making guesses about reference values of length velocity density and temperature From these it deduces and sets some initial values of variables including the DVELDPs i e the rates of change of velocity with pressure difference at every point It sets linear under relaxation factors for all variables including the DVELDPs It sets maximum values to the increments per sweep for some variables It activates whole field solution for all velocities The ideas underlying CONWIZ No theory of convergence promotion has been invented but the following ideas have provided guidance SIMPLE type CFD algorithms differentiate the momentum equations to deduce DVELDP AREA sum of CONVCOs MASS DT DTF In steady flow problems the true time step DT equals 0 at the start convection coefficients CONVCOs may be 0 so a finite false time step DTF is needed to keep DVELDP finite and so to prevent small pressure differences from producing enormous velocities but these are not needed once CONVCOs are finite because motion has started the much used constant DTF relaxation is thus now recognised as intrinsically bad it should no longer be used An initial DTF should be based on an initial estimate of what the CONVCO sum will be then linear under relaxation allows the momentum equation based values to take over gradually The fact that PHOENICS users have

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/numerics/converge.htm (2016-02-15)

Open archived version from archive - NUMERICAL CONVECTION SCHEMES IN PHOENICS

f 0 5 F c F d for Pe 2 12 F f F c for Pe 2 The cell Peclet number is defined as Pe r abs U f A f D f 13 in which A f cell face area and D f physical diffusion coefficient When Pe 2 CDS calculations tend to become unstable so that the HDS reverts to the UDS Physical diffusion is ignored when Pe 2 The HDS scheme is marginally more accurate than the UDS because the 2nd order CDS will be used in regions of low Peclet number 4 HIGHER ORDER DISCRETISATION SCHEMES 4 1 Classification of Higher Order Schemes LINEAR SCHEMES are those whose coefficients are not direct functions of the convected variable when applied to a linear convection equation Linear convection schemes of 2nd order accuracy or higher may suffer from unboundedness and are not unconditionally stable NON LINEAR SCHEMES analyse the solution within the stencil and adapt the discretisation to avoid any unwanted behaviour such as unboundedness These two types of scheme may be presented in a unified way by use of the FLUX LIMITER formulation 4 2 Flux Limiter Formulation The FLUX LIMITER formulation calculates the cell face value of the convected variable from F f F c 0 5 B r F c F u 14 where B r is termed a limiter function and the gradient ratio r is defined as r F d F c F c F u 15 The generalisation of this approach to handle non uniform meshes has been given by Waterson 1994 From equation 14 it can be seen that B r 0 gives the UDS and B r r gives the CDS 4 3 Linear Higher Order Schemes PHOENICS provides the following linear higher order schemes CDS Linear upwind scheme LUS Quadratic upwind scheme QUICK Fromm s upwind scheme Cubic upwind scheme CUS These are unified for implementation purposes as members of the Kappa class of schemes B r 0 5 1 K r 1 K 16 where K 1 gives CDS K 0 5 gives QUICK K 1 gives LUS K 0 gives Fromm and K 1 3 gives CUS These schemes are plotted in the Flux Limiter Diagram FLD in the next panel which takes the form of a plot of B r against r The two main regions of the flux limiter diagram are given by r 0 indicating an extremum and r 0 indicating monotonic variation Linear flux limiter diagram 4 4 Non Linear Higher Order Schemes The following QUICK based non linear schemes have been provided SMART bounded QUICK piecewise linear B r max 0 min 2 r 0 75 r 0 25 4 17 H QUICK harmonic based on QUICK smooth B r 2 r r r 3 18 UMIST bounded QUICK piecewise linear B r max 0 min 2 r 0 25 0 75 r 0 75 0 25 r 2 19 CHARM bounded QUICK smooth B r r 3r 1 r 1

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/numerics/scheme.htm (2016-02-15)

Open archived version from archive - TURBULENCE MODELS

where l m mixing length and has to be prescribed Near a wall in free flows in a duct Nikuradse s l m formula is usually appropriate in a boundary layer a l m ramp distribution is often suitable The model is suitable for simple shear layer dominated or pressure driven flows simple and computationally economical unsuitable for complex flows because it is very difficult to estimate the distribution of the mixing length cannot account for the transport effects of turbulence In PHOENICS set ENUT GRND2 EL1 GRND1 2 3 as appropriate and GENK T to activate calculation of the velocity derivatives LVEL TURBULENCE MODEL The LVEL turbulence model is applicable only when walls are present and it calculates ENUT via Spalding s law of the wall which covers the entire laminar and turbulent regimes thus where k 0 41 E 8 6 and The local parameter X is determined iteratively from a non linear expression involving k E and the local Reynolds number Re which is based on a typical local velocity and the normal distance to the nearest wall WDIS WDIS is calculated from an algebraic function of a scalar variable L and its gradient while L itself is obtained by solving with L fixed to zero in solids TURMOD LVEL activates the LVEL turbulence model The model is useful for situations in which fluid flows through spaces cluttered with many solid objects electronics cooling In such cases the grid density between nearby solids is often too coarse for the k e model to be meaningfully employed One Equation Models PRANDTL S K L MODEL Prandtl s one equation model may be summarized as follows where g i is the gravity vector l m has to be prescribed and the empirical constants are C m 0 5478 C d 0 1643 and PRT KE 1 0 1 equation models account for V s transport which is important in developing and relaxing flows and in flows controlled by free stream turbulence The main shortcomings concern L s its transport is not considered which is important in separated flows and it must be prescribed which is difficult for complex flows The trend has been to use 2 equation models TURMOD KLMODL activates the k L model and buoyancy production is introduced into the model by the following additional PIL commands PATCH KEBUOY PHASEM 1 NX 1 NY 1 NZ 1 LSTEP COVAL KEBUOY KE GRND4 GRND4 The formula for L s must be selected by setting EL1 to one of the built in options Two Equation Models The velocity scale V s is calculated from solution of a transport equation for KE The dependent variable of the 2nd transport equation is not usually L s itself but rather the variable There are many possible choices for the second turbulence variable resulting in different values of m and n These include Harlow Nakamaya KE EP model variants m 1 5 n 1 0 Saffman Spalding KE VOSQ model m 1 0 n 2 0 Kolmogorov KE OMEG model m 0 5 n 1 0 Rotta KE KE L s model m 1 0 n 1 0 The KE EP model has proved the most popular mainly because it does not require a near wall correction term THE KE EP MODEL The standard high Re form of the KE EP model employs the following turbulence transport equations The kinematic turbulent or eddy viscosity and the length scale L s are given by The model constants are C m 0 5478 C d 0 1643 PRT KE 1 0 PRT EP 1 314 C 1e 1 44 C 2e 1 92 and C 3e 1 0 The buoyancy production is 0 for stably stratified layers so that KE is reduced and turbulence damped For unstably stratified layers is positive and turbulence is augmented The constant C 3e has been found to depend on the flow situation It should be close to zero for stably stratified flow and close to 1 0 for unstably stratified flows The default value of C 3e is 1 0 and it can be changed by using SPEDAT SET KEBUOY GCEB R 0 3 TURMOD KEMODL activates the KE EP model and may be introduced by the following additional PIL commands PATCH KEBUOY PHASEM 1 NX 1 NY 1 NZ 1 LSTEP COVAL KEBUOY KE GRND4 GRND4 COVAL KEBUOY EP GRND4 GRND4 C 3e 0 may be effected by simply not using the COVAL statement Two equation models account for transport effects of V s and L s and the distribution of L s is determined by the model The model forms a good compromise between generality and economy of use for many CFD problems The models have limited universality mainly due to the rather simple modelling of the creation and destruction of the L s determining variable Hence the large number of k e variants involving a modified EP equation The use of the eddy viscosity concept imposes isotropy i e the same values are taken for the various Reynolds stresses This assumption may be too simple for example when body forces are present gravity rotation magnetic fields and in flows with streamline curvature The two equation models presented thus far are valid only at high turbulent Reynolds numbers and therefore they are not valid in the near wall viscosity affected region The Yap KE EP Model For separated and reattaching flows the standard KE EP model encounters difficulties e g as boundary layers develop towards separation the predicted near wall length scales become too large and the length of separation regions is underpredicted in flows past steps and obstacles Several variants of the KE EP model have been proposed for dealing with this deficiency including the Yap 1987 modification to the EP equation where and with YN is the distance from the wall and The correction can produce substantial improvements particularly in heat transfer predictions in separated flows The Yap correction is activated by using TURMOD KEMODL YAP The RNG KE EP Model The KE EP model is known to be too dissipative the turbulent viscosity in recirculations tends to be too high thus damping out vortices TURMOD KERNG selects the RNG k e model This attempts to correct this deficiency by using slightly different constants and by adding the following volumetric source term to the EP equation where and the dimensionless parameter The changes to the model constants are PRT KE 0 7194 PRT EP 0 7194 Although very good for separation and reattachment its predictions for jets and plumes are inferior to the standard model The Chen Kim KE EP Model TURMOD KECHEN selects the modified KE EP model of Chen and Kim which is another attempt to cure the over dissipative nature of the KE EP model This variant of the KE EP model uses slightly different constants and introduces an additional source term into the EP equation i e where C 4e 0 25 The changes to the model constants are PRT KE 0 75 PRT EP 1 15 C 1e 1 15 and C 2e 1 9 It predicts reattachment points and vortices almost as well as the RNG version but preserves the good behavior for jets and plumes of the standard version A low Re extension can be selected by TURMOD KECHEN LOWRE Two Scale KE EP Model TURMOD TSKEMO selects the two scale split spectrum KE EP model This model splits the turbulence energy spectrum into 2 regions namely the production P region and the transfer T region Spectral equilibrium is assumed between the T region and the region in which turbulence is dissipated i e the dissipation D region Two transport equations are solved for each of these regions namely KP and EP for the P region and KT and ET for the T region EP is the energy transfer rate from the P to T range and ET is the dissipation rate The total turbulent kinetic energy KE KP KT The two scale model embodied in PHOENICS employs variable partitioning of the spectrum so that the location of the partition the value of KP KT is calculated automatically by the model The transport equations for KP and KT are The transport equations for EP and ET are given by The turbulent eddy viscosity is computed from where and The functional relationship for determines the location of the partition between the P and T regions Note that for turbulent flows in local equilibrium P k ET and ET EP so that The model constants are PRT KP 0 75 PRT EP 1 15 PRT KT 0 75 PRT ET 1 15 C P1 0 21 C P2 1 24 C P3 1 84 C T1 0 29 C T2 1 28 and C T3 1 66 The model may be used in combination with equilibrium GRND2 or non equilibrium GRND3 wall functions The advantage of the 2 scale KE EP model lies in its capability to model the cascade process of turbulent kinetic energy and to resolve the details of complex turbulent flows better than the standard k e model The disadvantage is that it requires 4 turbulence transport equations as opposed to the 2 equations required for the standard k e model The recommendation is that the standard k e model or one of its variants be used in the first instance However in cases where these models are clearly giving poor predictions the 2 scale model should be used to see whether better predictions can be obtained Documentation see the PHENC entry TWO SCALE KE EP TURBULENCE MODEL Library Examples see cases T400 T406 inclusive Reynolds Stress Transport Model In general the Reynolds stress transport model RSTM involves the solution of the following turbulence transport equations 3 normal stresses 3 shearing stresses The dissipation rate EP 3 enthalpy fluxes and 3 mass fluxes per scalar variable The 3 normal stresses are always solved but the shearing stresses are solved only when the appropriate velocity components are solved The Reynolds stress transport equations may be written in the following symbolic tensor form T ij C ij D ij P ij R ij E ij where T ij transience C ij convection D ij diffusion P ij production R ij redistribution and E ij dissipation The R ij term is the most important term requiring closure in the RSTM It is commonly known as the pressure strain term and PHOENICS provides for 3 different closure models which are selected by IRSMHM as follows IRSMHM 0 selects the IPM model the default 1 QIM 2 IPY The IPM and QIM models are the most commonly used in the literature while the IPY model performs similarly but offers advantages for swirl flows These pressure strain models require wall correction terms which are provided by PHOENICS as default through the TURMOD command The corrections account for wall effects on the redistribution process The RSTM itself is activated from the Q1 file by the PIL command TURMOD REYSTRS DTFS WALL1 WALL2 The wall patches WALL1 WALL2 etc must be previously defined through laminar wall PATCH and COVAL settings for the velocity variables The TURMOD command will then create COVAL settings for the velocities and EP using GRND2 wall functions create PATCH and COVAL settings for the stress equation variables so as to represent the wall damping sources associated with the chosen pressure strain model create RELAX commands for all SOLVEd variables with a false time step of DTFS The TURMOD command also activates solution of the dissipation rate EP the normal stresses U2RS V2RS W2RS and as appropriate the shearing stresses UVRS UWRS and VWRS When heat and or mass transfer is activated the default is that the turbulent fluxes of energy and scalars are represented by a simple gradient diffusion model through the default setting IRSMSM 0 i e where and 0 255 The generalized gradient diffusion model can be activated by setting IRSMSM 1 in which case the Reynolds fluxes are computed from where the coefficient C t 0 3 The command IRSMSM must appear before the TURMOD command in the Q1 file The full Reynolds flux transport model can be activated by setting IRSMSM 2 before the TURMOD command in the Q1 file The flux transport equations may be written in the following symbolic tensor form where T i transience C i convection D i diffusion P i production and R i pressure scrambling The R i term requires wall correction terms which are provided by PHOENICS as default through the TURMOD command If H1 or TEM1 are stored TURMOD automatically sets STORE UTRS if NX 1 STORE VTRS if NY 1 STORE WTRS if NZ 1 and NOT PARAB Likewise if SC1 SC2 SCn etc are STOREd storage is provided as appropriate for USC1 USC2 USCn VSC1 VSC2 VSCn WSC1 WSC2 WSCn etc If IRSMSM 2 the procedure is the same excepting that SOLVE replaces STORE because transport equations will now be solved for each of the individual flux components In addition wall damping sources are created for these components and if U1 is SOLVEd and NX 1 then UTRS USC1 USC2 USCn are also SOLVEd The advantage of the RSTM is that it provides a more rigorous and realistic approach and it captures anisotropic effects automatically However RSTMs are much more complex than eddy viscosity models and they are computationally more expensive and less stable Documentation see the PHENC entry REYNOLDS STRESS TURBULENCE MODEL Library Examples see cases T600 T609 inclusive Near Wall Approaches The wall no slip condition ensures that over some region of the wall layer viscous effects on the transport processes must be large The representation of these processes within a CFD model raises 2 problems How to account for viscous effects at the wall How to resolve the rapid variation of flow variables which occurs within this region without using a very fine computational mesh There are two methods the WALL FUNCTION method and the LOW REYNOLDS NUMBER method For extensive documentation see the PHENC entries WALL FUNCTIONS and LOW REYNOLDS NUMBER TURBULENCE MODELS Wall Function Method The viscous sublayer is bridged by employing empirical formulae to provide near wall boundary conditions for the mean flow and turbulence transport equations These formulae connect the wall conditions e g the wall shear stress to the dependent variables at the near wall grid node which is presumed to lie in fully turbulent fluid Strictly wall functions should be applied to a point whose Y value is in the range 30 Y 130 where and U S is the friction velocity The user may elicit printout of Y in the RESULT file by setting YPLS T or WALPRN T in the Q1 file If in addition STORE YPLS appears UP is written to both the PHI and RESULT files The advantages of this approach are that it escapes the need to extend the computations right to the wall and it avoids the need to account for viscous effects in the turbulence model Equilibrium Wall Functions Equilibrium wall functions use log law formulae at near wall grid point and are invoked by setting WALLCO GRND2 in the Q1 file before the WALL or CONPOR command where E 8 6 k 0 41 and U r is the resultant velocity Wall friction is introduced via momentum source terms involving the friction factor where s max s turb s lam and and the local Reynolds number Roughness effects can be taken into account by setting WALLA to the equivalent sand grain roughness height in the Q1 input file PHOENICS then takes E as a function of the roughness Reynolds number For heat and mass transfer the wall flux QW of the variable f is calculated via the Stanton number St is determined from St max St turb St lam with with Pr T and Pr L denote the turbulent and laminar Prandtl Schmidt numbers In PHOENICS the equilibrium wall functions are activated by WALL INAME AREA IF IL JF JL KF KL 1 LSTEP COVAL NAME H1 GRND2 H1WALL or CONPOR NAME 0 type IF IL JF JL KF KL The flags that wall functions are to be generated on that surface Non Equilibrium Wall Functions A generalization of the equilibrium treatment to non equilibrium conditions is also in PHOENICS which may be invoked by setting WALLCO GRND3 in the Q1 file before the WALL or CONPOR command The log law is extended to non equilibrium conditions as follows where and The turbulent friction factor and Stanton number are now given by The value of KE at the near wall point is calculated from its own transport equation with the diffusion of energy to the wall being set equal to zero The mean values of P k and EP over the near wall cell are represented as in the transport equation for KE However in the formula for the near wall eddy viscosity EP is calculated from The non equilibrium wall functions will give better predictions of heat transfer coefficients at a reattachment point Fully rough wall functions PHOENICS also provides wall functions that are suitable for a fully rough near wall layer in local equilibrium defined in terms of the effective roughness height as for example in the atmospheric boundary layer These wall functions may be activated by setting WALLCO GRND5 in the Q1 file before the WALL or CONPOR command These wall functions are the same as those described above for GRND2 equilibrium turbulent wall functions except for the form of the logarithmic wall law which is now given by Ur UTAU LN Y Y0 k where Y0 is the effective roughness height This height is related to the size of the roughness elements on the surface e g a typical value for grass is 0 01m whereas for forest it is about 1m For heat and mass transfer the Stanton number is computed from a modified Reynolds analogy i e St s Prt Library example see case

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/general/turb.htm (2016-02-15)

Open archived version from archive - Two-phase flow

R1in RHO1in W1in VALUE IN P2 R2in RHO2in W2in VALUE IN W1 W1in VALUE IN W2 W2in wherein R1in R2in are the volume fractions in the incoming stream Obviously R1in R2in 1 0 Outlets If both phase densities are constant and both phases can pass in or out outflows can be simply specified using OUTLET OUT HIGH 1 NX 1 NY NZ NZ 1 LSTEP VALUE OUT P1 Pext VALUE OUT P2 Pext if Pext is not zero If only one phase is allowed to pass then a PATCH must be used with a COVAL for P1 or P2 PATCH OUT HIGH 1 NX 1 NY NZ NZ 1 LSTEP COVAL OUT P1 1E3 Pext if only Phase 1 is allowed to pass or COVAL OUT P2 1E3 Pext if only Phase 2 is allowed to pass Walls The WALL command assumes that phase 1 is the continuous phase and will only generate wall functions for this phase If wall effects are required for phase 2 these have to be explicitly added in Note that in the examples above P2 is only a flag for phase 2 mass sources it is NOT a real pressure and it does NOT have to be STOREd or SOLVEd Further Details Of Mass Flow Conditions For single phase flows the mass flow associated with a PATCH and COVAL for P1 is m Type C p1 V p1 P1 where Type is the geometrical multiplier specified through the PATCH type In two phase flows the relationship for P1 or P2 is m Type C p1 V p1 P1 for P1 V Pi inflows m Type R 1 C p1 V p1 P1 for P1 V Pi outflows i e it is the same for inflows but has an additional volume fraction multiplier for outflows This means that for equal phase COefficients and VALues outflows are proportional to in cell volume fraction Unfortunately in many practical cases outflow is more realistically represented as proportional to in cell mass fraction This can be achieved by setting the COefficient to the required COefficient multiplied by the in cell phase density The OUTLET command thus arranges that the COVALs generated are COVAL OUT P1 1E3 Pext or COVAL OUT P1 RHO1 1E3 Pext COVAL OUT P2 RHO2 RHO1 1E3 Pext or COVAL OUT P2 RHO2 1E3 Pext If either RHO is set to a GRND number both COefficients are set to 1E3 and it is left to the user to arrange for appropriate multiplication either by average outlet densities in Q1 or by the actual in cell densities in GROUND Within Phase Diffusion The normal within phase diffusion treatment for two phase flow is very similar to the single phase The diffusion coefficient for each phase is written as G f i RHO i ENUL PRNDTL f ENUT PRT f The term includes both laminar and turbulent mixing Note that there is only one laminar and one turbulent viscosity These are the phase 1 viscosities If different viscosities are required in the second phase the Prandtl numbers can be used to introduce the appropriate ratios If one of the phases normally phase 2 is a dispersed one say droplets or particles then within phase turbulent effects are absent PRT f should be set to a large number say 1E10 Phase Diffusion The phase diffusion term accounts for turbulent dispersion of particles or droplets due to turbulence in the continuous phase and is present in both the phase continuity and phase phi equations The diffusion coefficient is written as G ri RHO i ENUL PRNDTL R i ENUT PRT R i Phase diffusion effects can be removed by setting both Prandtl numbers for both phases to large numbers say 1E10 Turbulence Modelling The turbulence variables KE and EP are always phase 1 variables and the generation rate is calculated from phase 1 velocity gradients If required it is possible to modify the GREX coding to switch the turbulence model to phase 2 though this should not normally be necessary The presence of particles will normally have a damping effect on the turbulence field Strictly speaking this effect is not completely modeled in the built in k E model Only the generation term and dissipation terms are reduced by the volume fraction but there is no explicit damping This can be introduced via additional source sink terms Two forms of these terms are provided in the GX routine GXDISP which is activated by a PATCH whose name starts with KEDI COef and VALue for KE and EP are set to GRND1 0 0 or GRND2 0 0 to activate the models of Mostafa Mongia or Chen Wood respectively Interphase Source Term The interphase source contains the diffusive e g friction heat transfer and convective mass transfer related links between the phases It is formulated in the general form as follows General Form S IP f f i m ji f i int f i where f f i interphase diffusive transfer coefficient kg s m ji net mass transfer rate between the phases kg s maximum of 0 0 and the quantity enclosed f i int value of f at the interface between the phases and f i is a bulk phase value of conserved variable f Notes The unit of S IP is kg s unit of f i e if f is velocity m s the units of S IP are Newtons and if f is enthalpy J kg the units of S IP are Watts S IP obviously appears in the equations for each of paired variables f 1 and f 2 i e S IP 1 f f 1 m 21 f 1 int f 1 S IP 2 f f 2 m 12 f 2 int f 2 wherein m 21 m 12 The last argument of the TERMS command determines whether the interphase source Sip is active between a pair of variables Pairing is between odd and even numbered variables The Interface Value PHINT f The concept of the interface value f i int is the value of f i at the i j interface It can be a function of space time or local bulk phase values It is a property and not a variable obtained from a conservation equation In terms of a steam water system PHINT H1 and PHINT H2 would represent the saturation enthalpies of steam and water at the local temperature and pressure The near interface situation is shown on the picture above In the case where CO2 diffuses between water and air and C1 represents the concentration in the air and C2 represents the concentration in the water then PHINT C1 and PHINT C2 would represent the concentrations at the air water interface In other cases interphase momentum transfer for example the concept of an interface value is not helpful and it is more meaningful to think of direct transfer between bulk phase values Settings for PHINT f The PIL variable controlling f i is PHINT f There are many options for setting PHINT values Bulk to bulk transfer interface values In this case the interface value of current phase is made equal to the bulk value of another phase i e f i int f j Then interphase source becomes S IP f f i m ji f j f i where f i is current phase bulk value of f and f j is another phase bulk value of f Usage Recommended practice for momentum transfer For energy transfer via overall interphase heat transfer coefficients Activation No actions required i e PHINT f 1 default and PHINT f 2 default Known interface values If the values of f at the interface between the phases are known say f i int Value 1 and f j int Value 2 interphase source becomes S IP f f i m ji Value i f i Usage Whenever the interface values are thought to be the given constants as for example saturation interface enthalpies of steam and water at the prevailing pressure Activation PHINT f 1 Value 1 and PHINT f 2 Value 2 Notes A few built in option are provided for H1 and H2 to deal with pressure dependent saturation enthalpies at the interface They are described in PHINT entry of Encyclopaedia User specified formulae for interface value calculations are accepted via GROUND provided that PHINT f i GRNDn as appropriate Known difference in interfacial values If the difference in values of f is known at the interface between the phases say diff f i int f j int the actual interface values can then be expressed via phase bulk values as f i int f j diff f j int f i diff and resulting interphase source becomes S IP f f i m ji f j f i diff Usage Whenever there is a jump conditions at the interface as for example a constant heat of evaporation linking the steam and water saturation enthalpies Activation PHINT f 1 default and PHINT f 2 diff Known ratio of interfacial values If the ratio of f values is known at the interface between the phases say Value 1 f i int f j int the actual interface values can then be expressed via phase bulk values as f i int Value 1 f j f j int f i Value 1 and resulting interphase sources become where C 1 and C 2 are bulk to interface transfer coefficients to be defined in the following considerations Usage Whenever there is a floating interface condition as for example the ratio of specific heats for introduction of interphase heat transfer via enthalpy difference Activation PHINT f 1 Value 1 and PHINT f 2 default Interphase Transfer Coefficient FIP f f i kg s featured in S IP is actually the interphase flux per cell per unit conserved variable difference For built in treatment the specific interphase fluxes are made proportional to f IP which is termed the reference interphase transfer coefficient f IP or FIP kg s Ns m is the force per cell per unit velocity difference It appears as multiplier in the built in formulations of all specific f f i and m ji This multiplier is introduced because of the analogy between friction heat and mass transfer Reference Interphase Transport Coefficient CFIPS The coefficient f IP is controlled by the PIL variable CFIPS If CFIPS is any constant other than a GRND flag f IP is set as follows f IP CFIPS RHO1 R1 R2 Vol for CFIPS 0 0 f IP CFIPS RHO2 R2 R1 Vol for CFIPS 0 0 In the above Ri is the maximum of Ri and RLOLIM This ensures that as long as RLOLIM 0 0 the reference interphase transport coefficient remains finite even though the volume fraction of either phase falls to zero Vol is the cell volume This has reasonable consequences when friction is concerned as it ensures that very small bubbles travel at the same velocity as the carrying liquid Interphase Friction Built in Options Ignoring the interphase mass transfer for bulk to bulk transport of momentum the interphase source term becomes S IP f IP vel j vel i where S IP has units of Newton and f IP is now the interphase drag coefficient between phases in units of Ns m or kg s CFIPS of the above options in units 1 s is now proportional of the product of phase slip velocity and projected area per unit volume If STORE CFIP is set in the Q1 file then f IP is stored for output purposes and also for under relaxation purposes via the RELAX command The built in options provided are described fully under the Encyclopaedia entries CFIPS and INTERPHASE DRAG MODELS and so only a brief description of the major ones will be provided here CFIPS GRND7 selects the DISPERSED FLOW drag model f IP 0 75 C d RHO1 R2 R1 Vol V slip CFIPA D p wherein phase 1 is the carrier and phase 2 is dispersed Here CFIPA is a parameter used to define the minimum value allowed for slip velocity V slip between the phases Cd is a dimensionless particle drag coefficient and Dp is the particle diameter defined by CFIPB The PIL variable CFIPD selects the Cd correlation 0 Standard drag curve default 1 Stokes regime drag correlation 2 Turbulent regime drag correlation 3 Subcritical regime drag correlation 4 Distorted bubble dirty water drag correlation 5 Spherical bubble dirty water drag correlation 6 Ellipsoidal bubble clean water drag correlation Notes For CFIPD 4 or 6 the surface tension must also be defined via CFIPC If CFIPB Dp the minus sign activates the removal of R1 from the f IP formula a form which is often used for a dilute suspension of particles CFIPD 7 selects the particle fluidization drag model This model uses a f IP formula of the above form only when R1 0 8 otherwise it uses a formula based on the well known Ergun correlation for packed beds CFIPS GRND8 selects the same drag models as for CFIPS GRND7 but with phase 2 as carrier and phase 1 dispersed i e f IP 0 75 RHO2 R1 R2 Vol V slip CFIPA D p Notes The commands STORE SIZE REYN VREL CD APRJ EOTV WEB will allocate 3D storage for the particulate phase diameter its Reynolds number relative slip velocity drag coefficient projected area Eotvos and Weber numbers when appropriate The options are restricted to rigid particles they do not allow for particles droplets or bubbles to change shape and hence influence the drag The options do not account entirely for multi particle systems e g the continuous phase viscosity is used in evaluating Re rather than the apparent mixture viscosity The options are restricted to spherical shapes but can be extended to irregular shaped particulates by interpreting Dp as an effective sphere diameter or by using a different drag correlation The options do not account for compressibility effects in which case Cd becomes a function of both Mach number and Reynolds number Interphase Heat Transfer CINT In many cases the driving force is the temperature difference and available heat transfer correlations are based upon it With solving for enthalpies H1 and H2 and ignoring mass transfer the interphase sources become S H1 h 12 A S T 1 int T 1 S H2 h 21 A S T 2 int T 2 where h 12 bulk1 to interface heat transfer coefficient W m 2 K h 21 bulk2 to interface heat transfer coefficient W m 2 K A S total interface area m 2 The interface temperatures T 1 int and T 2 int can be eliminated via overall heat transfer coefficient and bulk to bulk temperature difference as follows S Hi 1 A S T j T i 1 C 1 1 C 2 where C 1 h ij and C 2 h ji This formulation is generalised so that the interphase transfer coefficient f IP is taken to be harmonic average of the bulk to interface coefficients C 1 and C 2 Thus f IP 2 1 C 1 1 C 2 Note if C 1 C 2 C then f IP C Built in provisions have been made to calculate C 1 and C 2 via their ratios to the reference transfer coefficients as follows C 1 CINT f 1 FIP C 2 CINT f 2 FIP Notes The default values of CINT are 1 If either CINT is set to a GRND flag there is no multiplication by FIP for that phase and the GROUND set value is used for C i directly Two options are provided internally for CINT H1 and CINT H2 see Encyclopaedia entry CINT for details Interphase Heat Transfer Built in Options In PHOENICS two phase mode it is common to solve for enthalpy It is required a source in Watts in the form S Hi h ij A S T j T i where h ij the heat transfer coefficient W m 2 K A S total surface area of particulates m 2 and the subscripts i and j refer o the current and other phase For the particulates of spherical shapes of diameter D p their surface area is given by A S 6 R2 Cell Volume D p The heat transfer coefficient is then computed from the local Nusselt Number Nu h ij D p k thus h ij k Nu D p where k heat conductivity of the carrier phase W mK Two options for computing Nu have been provided in PHOENICS CINT H i GRND7 This computes Nu from fit to the experimental data on heat transfer from spheres and it is valid over the complete Re and Prandtl number range CINT H i GRND8 This computes Nu from the so called Ranz Marshall correlation which is valid for 0 Re 200 Notes In both cases the particle diameter is set through CINH2C However if CFIPS GRND7 or GRND8 and SIZE is STOREd then this will be used instead The command STORE NUSS will store the calculated Nusselt Number in a 3 D storage for plotting Two variants of CINT H i GRND7 or GRND8 deal with phase specific heats as follows For constant and equal phase specific heats the following settings should be made If phase 1 is taken as continuous phase CINT H1 GRND7 or GRND8 CINT H2 1 E20 CFIPS GRND7 CFIPA minimum slip velocity CFIPB particle size If phase 2 is taken as continuous phase CINT H2 GRND7 or GRND8 CINT H1 1 E20 CFIPS GRND8 CFIPA minimum slip velocity CFIPB particle size For constant and different equal phase specific heats If phase 1 is taken as continuous phase only CINT H1 GRND7 or GRND8 CINT H2 1 E20 CFIPS GRND7 CFIPA minimum slip velocity CFIPB particle size PHINT H1 C P1 C P2 Interphase Mass Transfer CMDOT The interphase mass

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/ipsa/ipsa.htm (2016-02-15)

Open archived version from archive - combust.htm

can be shown graphically as Linear Relationships For Fuel Oxidant and Product Mass Fraction Linear Relationships For Fuel Oxidant and Product Mass Fractions 2 5 Reaction rate for fixed f and h 2 6 Features relevant to PHOENICS The f equation has no source term The h equation can has a source term resulting from absorption and emission of radiation but this is often small enough to be neglected If so temperature can be deduced from h via where the mixture specific heat is calculated from the individual species specific heats The mass fraction of oxidant mox can be deduced from The mfu equation has a reaction rate source of the form source of fuel const mfu x mox x function of T When the reaction rate constant is large there is no need to solve for mfu it depends only on f and can be deduced from This is the mixed is burned assumption The mass fraction of products can be deduced from mprod 1 mox mfu Density can be calculated from the Ideal Gas Law as where the mixture molecular weight W is calculated from the individual species molecular weights 2 7 Reaction rates in turbulent flows In turbulent flows with denoting time averages rate rctd does NOT equal rate rctd because of the non linearity of rate and mfu mox does NOT equal mfu mox because mfu and mox may be finite only at different times Since CFD codes compute only rctd mfu and mox but need rate new formulae are needed for turbulent flow No certain or universal method has yet been discovered but some guesses have proved lucky Experiment shows that reaction rates are often more greatly affected by local turbulence than by chemical kinetic factors The eddy break up model Reference 2 rests on the hypothesis that ONLY turbulence and fuel concentration affect the rate The first version was rate const x rho x mfu x abs vel grad The second version was rate const x rho x mfu x epsilon k Both work fairly well but they are only approximate 3 Implementation in PHOENICS The SCRS provided as a built in option in PHOENICS is a very simple version Only a single fuel and single oxidant have been provided for and there is no explicit dilutant stream It is of course possible to program in much more complex models in GROUND using the built in options as templates The f and mfu equations are created from two of the spare concentration equations The mox and mprod mass fractions are deduced as above and are then stored for printing and plotting For the rate controlled model various expressions for the reaction rate are coded into the subroutine GXCHSO They are activated by appropriate PATCH and COVAL commands for the mfu equation The calculation of density from the Ideal Gas Law but with a mixture specific heat is coded into GXRHO whilst the deduction of temperature from the SCRS definition of enthalpy can be found

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/general/combust.htm (2016-02-15)

Open archived version from archive - INTRODUCTION TO BODY FITTED COORDINATES

often forgotten by beginners Velocity Components With BFCs as with regular grids scalar variables are stored at the cell center i e the arithmetic mean of the cell corners Velocities are stored at the centers of cell faces The direction of the stored velocity is along the line joining adjacent cell centers The magnitude of the stored velocity is obtained by resolving the total velocity vector by projection into the direction of the line joining the cell centers The SOLVEd velocity U1 V1 W1 is therefore known as a RESOLUTE Cartesian Velocity Components Cartesian velocity components UCRT VCRT WCRT are also calculated and stored automatically These are the components of the total velocity in the directions of the reference Cartesian axes They should be regarded as being stored at cell centers Velocity Definitions Q is the total velocity vector U1 V1 and W1 are the SOLVEd velocity resolutes UCRT VCRT and WCRT are the deduced vector components in the Cartesian coordinate system UCMP VCMP and WCMP are the deduced vector components in the grid line directions Fundamentals Of Grid Generation Define Nodes The Cartesian Coordinates of all grid nodes cell corners must be specified Grid Orthogonality Lines joining cell centers should intersect the dividing face at angles as close to right angles as possible Angles less than about 30 degrees should be avoided if possible This can be checked by turning the Grid check ON on the Match Grid dialog Grid Spacing Grid nodes should be spaced finely where flow properties are expected to change rapidly with distance coarsely where they don t Smooth changes of aspect ratio are advised Grid Orthogonality The Steps In BFC Grid Generation The process of specifying a BFC grid involves Specify POINTS provide the Cartesian coordinates of key features of the geometry to be meshed Specify LINES join the points by lines which are divided into segments corresponding to the number of cells which lie along the line Lines may be straight arc segments or spline curves through defined points Create FRAMES link the lines to make 2 D FRAMES MATCH Grid to Frames match the grid plane to appropriate frames Form the 3 D Grid link 2 D grid planes to form a 3 D grid Points In the grid generator POINTS are created by clicking on the Points New Points buttons Mouse point allows the point to be placed with the mouse Type x y z brings up a data entry panel Once created points can be moved with the Modify Points button The PIL command for creating points is GSET P point name x y z where x y z are the Cartesian coordinates of the point Lines In the grid generator LINES are defined by clicking on Lines New Lines then New Straights New Arcs or New Curves then clicking on the points defining the line Once a line has been created the number of cells along it can be set or modified with the Modify Lines button The PIL

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/general/bfcintro.htm (2016-02-15)

Open archived version from archive - USE OF BFC'S

for a north boundary PATCH BFCIN NORTH 1 NX NY NY 1 NZ 1 1 COVAL BFCIN P1 FIXFLU GRND1 COVAL BFCIN U1 ONLYMS GRND1 COVAL BFCIN V1 ONLYMS GRND1 COVAL BFCIN W1 ONLYMS GRND1 COVAL BFCIN UCRT ZERO UIN COVAL BFCIN VCRT ZERO VIN COVAL BFCIN WCRT ZERO WIN Here UIN VIN and WIN define the Cartesian resolutes of the external velocity vector The density at the inlet is set as BFCA RHOIN if RHO1 is not STOREd or COVAL BFCIN RHO1 ZERO RHOIN if it is Potential Flow Solutions PHOENICS provides two alternative methods for the performance of potential flow simulations the DIRECT method and the DARCY method The DIRECT method involves using the TERMS command for the variable PHI i e the velocity potential to deactivate all terms of the transport equation except the diffusion term TERMS PHI N N Y N N N The resulting equation reduces to Laplace s equation for the potential in incompressible flow if If the whole field linear equation solver is invoked for PHI in the SOLUTN command the solution for the potential can be obtained in one sweep provided LITER is set sufficiently large On the second sweep the velocity field may be calculated from the gradient of the potential PHOENICS provides for this in subroutine GXPOTV called from GREX3 when POTVEL T Library case 116 exemplifies the use of this method for BFC F The method has been little explored for BFC T The Darcy method is the alternative way of solving potential flows in BFCs The DARCY method exploits the analogy that exists between the potential flow equations and the equations that govern flow in a highly resistive medium i e Darcy s law where k is a resistance coefficient having units of s 1 The similarity is exploited by solving the usual momentum and continuity equations but with DARCY T in the Q1 file The resistance coefficient k is set via the real PIL variable DARCON for which the default value is 1E4 Setting DARCY T causes the velocities to be solved cuts out the convection and diffusion terms for economy and introduces high resistances for the velocities However the user must also introduce the command SOLVE P1 It is important to remember that the velocity potential is The variable P1 computed and printed should thus be interpreted as Convergence of the Darcy solution is usually rapid 10 to 20 sweeps as the equations are linear It is possible to restart just the velocity fields from a Darcy solution using RESTRT U1 V1 W1 The pressure field will take the FIINIT P1 value but pressure correction will quickly get the right levels Rotating Coordinate Systems For a coordinate system rotating with constant angular velocity omega the apparent body force can be resolved into the local grid directions These forces are incorporated into the mathematical model as additional volumetric momentum sources Cases 744 radial equilibrium test and 747 geostropic balance provide examples of the use of rotating coordinates

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_lecs/general/bfcuse.htm (2016-02-15)

Open archived version from archive - GENERAL COLLOCATED-VELOCITY METHOD (GCV)

up a multi block case the user should specify grids for each block and use the READCO command to assemble the blocks Extra cells in computational space must be reserved to provide links between the blocks The geometry for these extra cells does not have to be defined Only those block faces which are linked to another block or blocks need an extra layer of cells Blocks may be connected in any arbitrary way to other blocks and or to themselves There is no limit to the number of times a block can be linked to another block or to itself All blocks must be right handed in physical and computational space It is not possible to recalculate the normal links during a transient run All links are defined as pairs of commands using the following format MPATCH n MBLnABCD typ1 if1 il1 jf1 jl1 kf1 kl1 1 1 MPATCH m MBLmEFGH typ2 if2 il2 jf2 jl2 kf2 kl2 1 1 where n m block numbers ABCD EFGH arbitrary letters to distinguish one patch from another typ1 typ2 patch types EAST WEST NORTH SOUTH LOW HIGH i1f i1l and i2f i2l patch boundaries in local block IJK co ordinates If the blocks share the same orientation in IJK space no further settings are needed If the blocks are rotated relative to each other in IJK space the block alignment must be specified This is done with the following command SPEDAT SET GCV MBLmEFGH C ABC where MBLmEFGH name of the second patch ABC character string which defines the rotation SET GCV C required settings The string defining the block rotation consists of three letters which may be any of E W N S H L The individual characters define the relative orientation between the N E and H faces of the first block of the pair of blocks and the current block respectively The default setting is SWL denoting natural connections The string ENL would define the connection NORTH EAST EAST NORTH and HIGH LOW Examples E N Block 2 S W N H W Block 1 E N Block 3 S S L Block 1 2 MPATCH 1 MBL1A NORTH MPATCH 2 MBL2A WEST SPEDAT SET GCV MBL2A C WNL Block 1 3 MPATCH 1 MBL1B EAST MPATCH 3 MBL3B NORTH SPEDAT SET GCV MBL3B C LNE E N Block 2 S W N H W Block 1 E N Block 3 S S L Block 2 3 MPATCH 2 MBL2C WEST MPATCH 3 MBL3C HIGH SPEDAT SET GCV MBL3C C SLE A full description of the multi block settings is given in Multi block Grid and Block Linkage Definition for GCV 4 Sliding Grid Settings One group of blocks is allowed to rotate relative to another group of blocks The sliding interface must always be between one stationary and one rotating block A patch with the name MBS should be specified for each rotating block These patches should cover the complete rotating blocks including any dummy edge cells

Original URL path: http://www.cham.co.uk/phoenics/d_polis/d_enc/enc_gcv.htm (2016-02-15)

Open archived version from archive

web-archive-uk.com, 2016-10-25