Home

SwelFlo User Manual Version 2.01

1. Line 98 Coll Output Curve Rock Properties Simulation Direction f f A i Produce an Yes Allow Heat Flow Yes Topdown Bottomup Output Curve No To From Country 9 No No of points in Thermal Cond 2 0000 Output Curve 0 W K m j Range of flowrates 0 TETU ETT for aip ov a0 kg cubic m a Heat Capacity 1000 0 G kg K Range of 0 000000E 00 Hor bottomhole P s Time Se 0 60480E 06 touse Bara 0 000000E 00 Spacing for Total well length m 1000 00 data output fm 50 hera ma sv 5 Error Reporting Error reports from SwelFlo simulations Current Directory C Users Markm Documents SwelFlo GenericWell project Figure 1 SwelFlo windows visible on a first startup 3 1 A Topdown Run The built in template is for a topdown simulation of the downhole profile of a producing vertical well of length 1000m with wellhead pressure P 8 bara flowrate Q 20 kg s temperature T 169 75 C and mass fraction of carbon dioxide XCO 0 001 Press the Run button hear the bottom of the main window or the Run command in the menu of that window to perform a topdown simulation starting at the wellhead values specified It should run very quickly and show the resulting pressure versus depth in the well The downhole pressure profile should look like that portrayed in Fig 2 Clicking on the T and V tabs reveals temperature and mass weighted velocity spinner v
2. 2 _ y2 G lt v gt Spu 1 S pju G Be E Zo Spy 1 S p 52 This is only useful for constant cross sectional area since otherwise the term G should not be inside the derivative 25 The frictional pressure drop is modelled using a two phase correction factor dj which is described below The term is computed as in 2 66 and 7 31 in Chisholm ane dP _ foG G dz lio 4r py where fp is the Darcy friction factor given implicitly by the Colebrook equation see next subsection and G is computed in the usual way for two phase flow G Q A p l Su py Sty C 3 1 Friction Factor Solutions are needed to the implicit Colebrook White equation 9 to find for turbulent flow the Darcy friction factor fp that solves 1 210 q 2 512 SS g Tp 10 3 6D Revyfp where D is well diameter and Re puD p where u is fluid velocity Newton iteration can be used to iterate to find a solution as is done in GWELL An alternative adopted in SwelFlo is to use Serghides Solution 31 accurate to within 0 0023 and fast since it is explicit It is given by A MO ote 9 B glom Ga 10 C2 21g 555 11 fo a a 12 where e is roughness in the same units as wellbore diameter D As suggested in 86 the extended Colebrook equation is used in cases where roughness over diameter exceeds the value 0 05 as can happen with the Duns and Ros correlation for mist f
3. 3 4 Output Curve Since a productivity index is known for the feedpoint an output curve may be computed Turn off the Downhole Profile option this is just for better speed here it is generally possible to compute both a profile and an output curve in the same run and click Yes for the Produce an Output Curve option in the SETUP tab Now RUN causes the simulator to compute only an output curve with just 12 successful points if bottomhole P range is the default of zero as illustrated in Fig 5 I oc Simulation O Data from header Q kg s Figure 5 SwelFlo first steps showing the output curve computed from the template setup with no bottomhole pressure range specified The number of points in the output curve may be refined by looking at RESULTS Show Detailed Output and noting the range of bottomhole pressures near the end of the output 13 Q kg s 4 5 5 5 6 P Bara Figure 6 SwelFlo first steps showing the improved output curve computed from the template setup with refined bottomhole pressure range file that gave successful output points Then for example setting a range of 40 to 52 bara here for bottomhole pressures and running again gives a much faster computation and many more points on the output curve as illustrated in Fig 6 Similarly setting a range like this also affects bottomup computation of a downhole profile and gives a much faster bottomup match to wellhead values if
4. If solving from the top downwards values P and TT for pressure and temperature are known at node n Note that when carbon dioxide is present these variables may be used in both single phase and two phase regions since saturation varies with temperature for a gassy two phase region Then it is desired to solve for the next values P and T 1 A nonlinear Newton Raphson search method modified by Armijo s Rule is used for this search which varies P 1 and T 1 until momentum and energy values at node n 1 match those at node n to within absolute errors of 1 0E 08P and 1 0E 09 respectively where E is the total energy at node n Because of the Newton search the method is implicit Because simple arithmetic averages of densities and friction values at the two nodes are used for the region between them the numerical method is equivalent to an implicit method called the trapezoidal rule This means the method is second order accurate and is the most accurate of the second order A stable linear multistep methods It also means the direction of iteration topdown or bottom up should not affect stability 55 Figure 16 A sketch of two sections of the wellbore set at different angles to the vertical with subsections In particular integrating equation from node n to node n 1 it follows that for constant area A that is constant G mu dP ff d oa aP J r dz zie lt u gt lt p gt gsin 0 fro Fa dz whi
5. floating point number the time in seconds since started discharging for transient heat flow calculations between well and adjacent country rock reservoir suggested is one week 6 048E05 seconds line 13 NDATAF NOTOPS NUSEC integer number of wellbore casing sections that follow NUSEC must not be more than 500 and must be at least one 44 lines 144 NDATAF NOTOPS to 13 NDATAF NOTOPS NUSEC SECL RAD EPS DELZ DEV TRUEV CASINGCHOICE section length MD radius of this section of casing roughness subdivision length desired for computations MD deviation angle from horizontal in degrees vertical is 90 0 the true vertical extent of this section casing choice All numbers are floating point separated by space s except Casing Choice which is an integer There should be NUSEC of these lines The units for lengths and pipe roughness in the input file are m whereas in the displayed grid windows these depend on the user s choices for units If a positive value is entered for true vertical extent this is used together with the measured length arclength along the wellstring to calculate and override the deviation angle If a value for CASINGCHOICE integer greater than one is entered the values on this line for radius and roughness are over ridden by values read in from the database file CasingLinerDatabase csv in the main SwelFlo directory The number for Casing Choice is taken to be the line number in the csv file The number on
6. P T V Since a run always generated an output file but does not always generate a downhole profile there may be gaps in the numbering of the plot files so that they always have the same number as the corresponding output file The text file SWFoutcrv001 csv is created in the current or working folder if an output curve is requested and obtained containing an output curve ready for Grapher input with flowrate in column 1 wellhead pressure in column 2 flowing enthapy in column 3 then any flowrates and WHPs and flowing enthalpies in columns 4 5 and 6 that may have been provided in data files Units used for user provided flowrates pressures and enthalpies depend on what the user has chosen for this project Numbering is chosen to match that for the corresponding SWFoutput file This file is also in the correct format 36 to be added to the Data files that SwelFlo reads and displays in graphs along with simulated results In that case SwelFlo will only read the first two columns to extract simulated flowrates and WHPs This allows a user to compare output curves from different projects on the same graph Or from different runs in the same working folder differentiated by number 4 Three files are created summarising simulation results if a downhole profile is requested suitable for immediate input as data files say for comparing two different projects under the same well directory There is a pressure file a temperature file and
7. SIMRANOK the most recent simulation ran OK is used to decide whether to display graphs and results display if true GOTOUTC computed an output curve OK is used to decide whether to display an output curve READYTOPLOT is set to true if there are data points and or simulation points ready to be plotted GETPROFILE is true if the user wants a downhole profile computed 42 ALLOWHELO is true if the user want to allow heat to flow between the casing and the country QMAXOC and QMINOC set the upper and lower limits for wellhead flowrate in the plot of the output curve in the current units chosen for flowrates kg s or t hr NSIMOUTCPTS records the number of output curve points that have been successfully simulated in the latest simulation SPACNG is the number of metres depth spacing MD desired between points output to files after simulation for future simulations to easily read in and plot as data alongside simulated results For P T and V versus depth Files are called PDataSim txt TDataSim txt and VDataSim txt where the asterisks stand for numbers 001 to 999 corresponding to the simulation number line 5 NDATAF integer The number of data file names to be read from this input file up to a maximum of 20 files Each data file may contain up to 10 000 data points line 6 data file namel character type of file integer The name of the first data file to be read Currently up to 256 characters This must be
8. are Orkiszewski Armand Hadgu with just the phase velocities affected full Hadgu with momentum conservation also affected Duns and Ros with just the phase velocities affected Duns and Ros with momentum conservation also affected and the Hagedorn and Brown correlation These correlations were originally developed for the petroleum industry rather than the geothermal industry More details on these correlations and how they have been modified for this simulator appear in the appendices to this manual 18 It is unclear which correlation is preferable despite a number of studies comparing them The Hadgu correlation is likely to be the basis for the WELLSIM simulator and is based on Hadgu s thesis results 19 and also draws on the GWELL code 1 5 that was available in the early 1990 s The Orkiszewski correlation is quite popular in geothermal wellbore simulation 4 5 Wellhead Conditions a a eee O jaca PE WELLHEAD CONDITIONS jf e E a File Edit Search Save SaveAs Clear DSS Q11Da P Bara Q kg s T oC H kJ kg flowing H kJ kg x flowing X xco2 8 00000 1 0000 166 153 0 7143851E 03 0 88379030E 03 0 67787631E 02 0 10062390E 00 0 0010 0 0010 Row 1 Column 2 i Figure 8 The Wellhead Conditions grid window in SwelFlo Fig shows the Wellhead Conditions grid window that appears when the Wellhead button is clicked It is important to set the wellhead pressure flowrate XCO2 and one
9. regime He chooses the correlation dP Ge fo pat 79 80 where the Darcy friction factor fp is a function of Reynolds number D Re pM 7 pu and is given by the Serghides approximation to the Colebrook equation accurate to 0 0023 A 2logio sto a 80 B 2logi ote 81 Ca 2005 555 Re 82 fo 4 2 83 where e is roughness in the same units as wellbore diameter D F 4 4 Gravitational Pressure Drop For bubble flow Hadgu uses the same form as SwelFlo lt p gt gsin 6 For slug flow this is modified slightly to lt p gt 1 8 gsin This is done in SwelFlo by using 6 S028 For annular flow the form is lt p gt gsin 6 and the saturation value used in lt p gt is computed using the Baroczy formulation in eqn 63 F 4 5 Momentum Flux Term The momentum flux term is also called the accelerational pressure drop in Hadgu s thesis In bubble flow it is the same as SwelFlo usually uses for Armand and Orkiszewski correlations G lt v gt in the third form presented in eqn 8 Hadgu makes the same sign error on the momentum flux term as is made in the GWELL manual In his PhD thesis the error originates in the transition from his correct eqn 8 5 to his incorrect eqn 8 6 In slug flow Hadgu uses the results of Fernandes et al the pressure gradient required to accelerate downflowing fluid between the slug and the wall to the speed Uzrpg gives the gradient of momentum
10. reservoir pressure is provided and is positive The productivity index is important for reliable output curves since it provides a measure of feed productivity that takes proper account of changes in the viscosity of the fluid at the feed as the enthalpy pressure and saturation of the fluid change at different flowrates Simply using linear or quadratic drawdown parameters is only accurate near the flow conditions that produced them they should vary with viscosity but don t while the productivity index is independent of viscosity Using a constant productivity index gives more accurate and reliable output curves since it allows for the effects on flowrate of any changes in viscosity at feedpoints as operating conditions change In bottomup mode if an output curve is requested if a productivity index is provided then that is used at a feed If flowrate is specified or linear drawdown at a feed then a productivity index is calculated from the downhole profile and used to produce the output curve If quadratic drawdown is specified the Forcheimer term A is computed as well as the 27 productivity index and both are used to compute the output curve Otherwise A is set to zero giving linear drawdown at that feed If a dry steam power law is specified at a feed no attempt is made to compute a productivity index there an output curve is still generated Dry steam accounts explicitly for variations in steam density so the productivity ind
11. the value of S output is Shu s otherwise for 8 large enough this inequality ceases to hold and the regime is slug flow with S output as Sslug Note that experience shows that it is best if the slope SL is not too small as rapid changes in vapour phase velocity can then result at saturations near bubble slug transition These cause apparently random convergence issues with the Newton step for solving the conservation equations depending on whether this transition region is encountered or not In practice a value of 0 3 is found to give reasonably smooth curves without the sharp kink seen in Fig with a value of 0 01 5 Otherwise if Ls lt vyp lt Ly the regime is transition to mist or annular flow with relative velocity of vapour changing linearly in vgp from slug value to zero The 68 0 8 0 6 y 0 4 0 2 0 0 2 0 4 0 6 0 8 l B Figure 21 An illustration of finding 6 where y P and y Lp intersect plotted against 3 for particular values of densities and flowrates calculations made are a ratio ra 35 Uy vr ratio Uy 36 GX how S p 37 Uv Pv 6 Otherwise the regime is mist with u vr since both phases move at the same speed and GX flow UvPv S An illustration of how the fluid looks in the four flow regimes bubble slug transition mist in a snapshot view of the flow taken from 29 is presented in Fig 22 The result of applying this procedure is illustrated
12. 1 0 15 0 0018 CASING 4 1 2 M 9 5 4 1 4 5 0 0018 CASING 6 5 8 B 28 0 5 8 6 6 0 0018 Changes or additions to the Casing_Liner_Database csv file will not show up in the simulator unless it is re started after the changes are made since the simulator reads this file only on startup The Edit Row menu item in this grid window allows the user to insert or delete entire rows or clear their contents A maximum of 500 rows is currently allowed The Clear menu item allows the user to clear the contents of the current cell the selected cells or the current row There is no need to Save this window before closing There is no prompt when the user closes since setup data is not lost However there is a Save and a Save As menu item on this window which allows the user to either save to the current Project directory the default wellhead conditions csv file SWFCasing csv or to save the setup as a SWFCasing csv where could be digits or to save the csv file with any other name if desired This is to give some flexibility in saving casing geometry setups for example if the user wishes to import read open them from another project from within the corresponding Casing grid window Depths and radii are saved as scaled variables using the scaling currently chosen by the user These saved csv files are what may be opened by the CASING GEOMETRY window menu File gt Open Trying to Open an incorrectly shaped csv file may result in the exe
13. 107 becomes Apo n AVin Use which looks the same as Torrens equation with the exception that p replaces p Both this comparison and the one made just after equation indicate that Torrens computation and the one used in Hagedorn and Brown s original paper become very close to each other when vapour flow dominates since then R 0 and pn gt py This is when the acceleration terms are significant The gravitation pressure gradient in the Hagedorn and Brown correlation uses pseudo saturation S in place of saturation in the usual form e Sps 1 S pig sin 96 F 7 1 Pseudo saturation Correlation Hagedorn and Brown originally presented graphical correlations for pseudo liquid holdup 1 S We use our own slightly adjusted versions of the DeGance and Atherton fits S 1 Y exp 3 6372 0 8813 In Nnoia 0 1335 In Mpota 0 018534 In Nhota 0 001066 1n Mnoia provided that if Nnoiq is outside certain bounds Nhold g _ 1 pa 211 5480627974693 Nhold gt 1000 1 0 00131950127 gt Mia lt 0 1 I have slightly altered these conditions from those in 11 to ensure that a smooth transition is seen at the bounds to more than eight decimal places as Npoig is varied The resulting graph of the ratio of liquid holdup to Y R w 1 S y is plotted against Nho in Fig 33 and looks the same as the equivalent plot in 11 Also it was found ne
14. P increases The bubble slug transition region slope used was SL 0 01 70 F 3 Hadgu Correlation The Hadgu Correlation is described in detail in 19 It is likely to be the main correlation in WELLSIM This section summarises the equations comprising the correlation as coded into SwelFlo It affects the shapes of the regions describing the nature of the flow in the wellbore bubble slug churn annular and the relative velocities of vapour and liquid phases in those regions The option Hadgu v only means the forward Hadgu correlation is used to calculate the velocities of vapor and liquid phases depending on the flow regime This is in contrast to the option Full Hadgu where in addition the Hadgu correlations affect the computation of momentum balance or pressure gradient by using frictional pressure drop formulae from Hadgu s thesis that differ from Chisholm s modified Martinelli formulae that are used for the Armand and Orkiszewski correlations and when v only options are chosen F 3 1 Phase Velocities The correlations that determine vapor and liquid phase velocities in line with Hadgu s thesis 19 are described here These correlations are invoked by choosing Hadgu v only in the correlation text box in the Input tab This choice leaves the momentum formulae unchanged from those described in subsection First we describe the forward problem with saturation a static property given and with flowi
15. Pa Q CE cnet Usually n 0 75 and C ranges from 0 4 for typical producers to 0 75 for good producers to 1 5 for the very best producers Note that in this formula however pressure in the reservoir and in the well are first converted to bars before squaring etc and the calculated flowrate is taken to be in t hr and is then converted to kg s internally afterwards So C and n must be chosen accordingly for pressure in bars and flowrate in t hr Dryness may be specified but this does not really make sense for a dry steam feed Drawdown Feeds The extra parameters to be provided are c the linear drawdown in bars t hr which is in fact one over the usual productivity in t hr bar and d the nonlinear drawdown term for flowrate in t hr and pressures in bara in the drawdown equation Pa Pwel cQ dQ Q A 25 To get just linear drawdown you can set d 0 in the input file in column nine Note that here as for dry steam pressures are first converted to bara and the calculated Q is assumed to be t hr so is then converted to kg s afterwards As noted above there is no need to Save this window before closing There is no prompt when the user closes since setup data is not lost However there is a Save and a Save As menu item on this window which allows the user to either save to the current Project directory the default wellhead conditions csv file SWFFeeds csv or to save the setup as a SWF Feeds csv where could b
16. Y is the productivity index which is needed if drawdown information is to be independent of viscosity changes at a feedpoint for output curve simulations v m s is the effective kinematic viscosity for flow between the feed and the reservoir and A is a nonlinear Forcheimer term If A is set to zero this equation is just Darcy s law for linear drawdown at a feed Forcheimer s equation is an extension of Darcy s law to situations where turbulent flow near the wellbore leads to lower flowrates and nonlinear dependence on pressure drop 28 Please see Appendix for details about effective viscosity and productivity index and for a discussion on their relationship to the more usual productivity in t hr bar The productivity index and perhaps also A for quadratic drawdown is needed for correct output curve calculations since the effective viscosity changes with pressure and especially also if flashing occurs at a feedpoint at higher flowrates and this changes the flowrate If FICODE is 2 flowrate fixed the whole line should contain F_DEPTH D FICODE I RESV 1 1 RESV 3 D RESV 12 FDRY D RESV 6 D FLOW I where the new variable is FLOW I the fixed flowrate desired at the feedpoint in kg s floating point Positive flow is for production If FICODE is 3 dry steam power law the whole line should contain F_DEPTH I FICODE I RESV 1 1 RESV 3 D RESV 12 FDRY 1 RESV 6 RESV 10 1 RESV 11 1 The new variabl
17. a downhole profile is requested If the user wants to examine the full downhole profile at any point on the output curve they can specify the corresponding wellhead values for that point and request a Downhole profile to see graphs of P T and V versus depth giving those wellhead values 4 The SETUP Tab This section contains more detailed information about the working of the SETUP tab 4 1 Simulating Downhole Profiles This section gives more detail about the ways to simulate and obtain a downhole profile that is values of pressure temperature saturation enthalpy density and flowrates versus depth in the well In the SETUP tab the user chooses Yes for Produce a profile chooses a Flow correlation and a Simulation Direction Then the four buttons labelled Wellhead Reservoir Temperature Casing and Feed Points may be used to set up the simulation This requires setting wellhead values a reservoir temperature profile if heat transfer between wellbore and country is to be included casing geometry and feed point properties These four buttons are discussed in more detail in later subsections The ranges that are allowed for the simulation variables are listed in Table 1 and the maximum dimensions of variables are listed in Table 2 14 Variable min value max value units pressure 1 0 1000 Bara temperature 1 374 C well length 1 5000 m casing section length 1 5000 m wellbore radius 0 10 m roughness 1 10 m subdivision
18. and fluid velocity is m s Any output test data provided may have up to 3000 lines each of which must take the form Mass flowrate at wellhead wellhead pressure flowing enthalpy Inclusion of enthalpy can be avoided by ending the line with a or in fact it can simply be omitted and there can be just two values per line flowrate and whp in the units currently chosen by the user Note that the csv output curve files written by simulations are compatible with being read in as data files for other simulations without alteration The first line gives a read error but such lines are ignored by the simulator which just tries to read the next line The wellhead information flowrate wh pressure given in the first line from each of the downhole profile data files is also added to output test data for display in the GRAPH of the output curve if an output curve is requested by the user The Re Plot commands found on the DATA INPUT tab and on the INPUT tab will read in any data files added before re plotting the most recent simulation 30 6 Results Tab A text view of the results of the last simulation run in a folder is presented in the RESULTS tab showing wellhead values flash depth if any and feedpoint flow properties either imposed as boundary conditions or obtained as simulation results depending on whether the simulation was topdown or bottomup A button in this tab also allows the user to view the latest output file generated
19. and temperature P and Tn in node n the aim is to find values in the next node n 1 such that Ful Pn 1 n 1 Pas Tn 0 17 Fg Pn da Ln tel ned a 0 18 where steady state momentum conservation says Gn 1 Gn 2 Fm Pny1 Tapi Pah Tn Php Phd lt U gt n41 lt U gt n lt p gt nti lt p gt n gsin Az pef J ome 2 a o 97 Nr Dire and steady state energy conservation says Se Sal Bos 50 0d Qul Az gsindAz 20 Fr Pris Lets Pas Tn lt h gt In two phase regions the variables actually used in the Newton step are T and X although T and P could still be used since there is always carbon dioxide present in the simulated fluid The method used is Newton s method in two dimensions modified by Armijo s Rule One step of this iterative method is described as follows an initial guess P and Tp at values for Pa and T 1 is updated to the values Prt Pay AJT Fir Fnet Tapi Pas Tn 21 That Toyi Fg P 1 gt Taris Eha Tn n where the Jacobian matrix J of first partial derivatives is Po Ta AP Ty a IF E P oTe 33 Pea Ca Fe dl 1 and where the size of A initially one is halved if the resulting updated values of pressure and temperature do not reduce the values of Fy and Fg to less than their starting values This reduction by the factor A is Armijo s Rule which
20. average enthalpies of the different components weighted by their individual mass fractions H 1 Geo Hy GucorlHcos Heoin 5 21 and H 1 Gyoo9 He yco2 Hco2 5 22 where a mass fraction CO H enthalpy kJ kg The subscripts 1 v w s 1CO2 vCO2 and soln stand for liquid phase gas phase water steam CO in liquid phase CO in gas phase and solution respectively 27 The total mass fraction of the gas phase x is calculated using a mass balance on CO Q Te co2 1c02 5 23 vco2 T 1co2 5 1 5 VISCOSITY 5 1 5 1 Carbon Dioxide CO Fits to the viscosity of CO as a function of pressure and temperature are based on the data tabulated by Vargaftik 1975 For the viscosity Uco in poise Moon A BT CT DT ET x 1078 5 24 where T temperature in C The pressure dependent coefficients are found from the linear interpolation between the tabulated values in Table 5 3 5 1 5 2 Mixture The viscosity of the liquid phase is assumed to be approximately equal to the viscosity of pure water since the amount of dissolved CO is small For the gas phase a weighted average is used Hi Hgo 5 25 By coo Mco2 1 Uco Hey 5 26 5 1 6 SURFACE TENSION The effect of CO on the surface tension of water has been studied by Heuer 1957 as part of his Ph D dissertation He measured the interfacial tension of H O CO at different temperatures and partial pressur
21. been corrected in SwelFlo Feedpoint viscosity is computed as an integral average over inverse values viscosity from wellbore to reservoir pressures as described in 26 126
22. depth mMU True vert depth m KEI fr p 4 Column 1 Figure 9 An example of the Reservoir Temperatures grid window in SwelFlo 4 7 Casing Geometry The Casing button opens a CASING GEOMETRY window which is essential to the simulation and the user enters the details of one or more up to 500 sections of casing length measured depth MD of the section an optional Casing Identifier selection radius roughness and angle to the horizontal in degrees so that 90 0 is vertical or alternatively the vertical extent of the section If vertical extent and length measured depth MD along the casing are provided and are nonzero the simulator will calculate the angle and override any angle values given The newly calculated angles will be displayed if the Preprocess button is pressed If vertical extent is set to zero or less it is ignored f amp CASING GEOMETRY File Edit Search Open Save SaveAs EditRow Clear wa aa Sec length mMU Casing Identifier radius m roughness m subdiv length m angle to horiz 0 vertical extent m no selection 0 1000000 0 00000 20 00000 90 000 0 00000 1 Column Figure 10 Screen snap of the Casing Geometry grid window The user also specifies the size of subdivisions of a section which is the mesh used for imposing conservation laws in the simulator If subdivisions are set to too large a size the conservation equations may fail to hold within the desired accuracy
23. ds o ae e ee ed 3 First Steps 4 5 31 A Topdow a RUN Je ec e a da Ee da da 3 1 1 Consistent Wellhead Values 0 000000 004 32 ne WE gem ae a eee Sw oh le Beg ee ID OOS he oe Ae eee Ge GE A EI ANA PEA Se oe Ee Sees The SETUP Tab 4 1 Simulating Downhole Profiles 0 o 4 2 Topdown and Bottomup modes cosmos een die ae a Ek Ad eS 4 2 1 A single bottom up simulation o e e ER A IE III 4A Correl tions IA 4 5 Wellhead Conditions lt 5 a2 4 oo oe doe haw des 4 5 1 Wellhead Thermodynamic Setup Priorities a a a ao a th ad GM Soe A He eine as ee ak OS re 4 7 Casing Geometry 4 6 26 04 4 9 ob ed Sed Oe BAD Se ee SE we me Od O E ee a na ea 4 9 Pre Process Run Stop Re Plot o 0 o LEAR A AO 4 11 Rock Properties i 45 642 25 44542 4 428 AA Er aR 4 12 Depth Spacing for data Output lt vie a ieee D4 eG ee Eo Data Input Tab 6 Results Tab Graphs Tab 8 Menu Items IA lt aeeie oe oe ee ee ee es Oe Oe Oe Oe ee yee 8 1 1 New Project Folder ceci ew aora eR aa ER 8 1 2 Import Existing SwelFlo setup o 0002 eee oe EE E A ee e a ae ad ad ee 8 1 4 Switch to Existing Project o eee ee ee AA AE a hs ye loa a oe Bee A 10 11 13 14 14 17 17 18 18 19 19 20 21 22 26 26 28 28 29 31 31 AE A hee amp De eee we eee we ee ee oe wre 34 8 3 Pre Proces
24. everything in the simulation setup has already been saved to SWFinput txt 8 1 6 Change Units This allows the user to choose different units for displaying results in the RESULTS and GRAPHS tabs and in the main output file These units are saved in the text file SWFunitsChoices txt in the project directory so applying only to the current project The chosen units are also the units that are assumed to apply to any data files used by the project and are also used in all saved csv files and output files in the project directory Internal units are SI units these are also the units saved to the project file SWFinput txt when writing pressures and flowrates etc Changing units after setting up data files may lead to incorrect display of data files since they are assumed to use the current unit system If plots are re displayed by using the Re Plot button immediately after changing units they will also show incorrect units since the output file being used to generate the plot is still in the previously chosen units A rerun will fix this It is up to the user to ensure that data files use the same units as the simulator 8 2 Help This menu option offers the option to display a simple pdf help file that contains this Manual It also gives access to an About menu item which will display the version number 34 and a brief message about the simulator 8 3 Pre Process Run Stop These have the same effect as the corresponding buttons at th
25. feedpoint which the code ensures occurs at a node a virtual node is added at the same depth in which to mix any incoming fluid with the wellbore fluid Mixing is assumed to occur at the wellbore pressure such that flowing enthalpy is conserved and any inertial effect of influx or efflux of fluid at a feed is not known Then with regard to the sketch in Fig 17 conservation of mass of water at a feed at node n takes the form Qn 1 E Qn Qf i n 1 Figure 17 A sketch of a section of the wellbore with a feed at the depth of node n for a topdown simulation An extra node of zero volume is added between node n and n 1 in which to mix any incoming reservoir fluid at the wellbore pressure P Details on how the flowrate Qs is computed for single and two phase flow are to be found in subsection 4 8 Four different types of flow response computations are offered at a feed productivity index fixed flowrate dry steam and linear or quadratic drawdown While the above equation holds no matter what the signs of the flows the values of XCO2 carbon dioxide mass fraction and flowing enthalpy lt h gt of the fluid in the unknown node n 1 do depend in the same way on flow directions There are six different possibilities illustrated in Fig 18 60 Case One Case Four Figure 18 A sketch showing the six different flow situations that are possible at a feed point The arrows in this figure show the ac
26. flux rather than the actual momentum flux as dP U U En Urma arg LIB M LS y 84 AL ic lsu Note that this is already in divergence form it is not an estimate of momentum flux but rather of the divergence of momentum flux 81 The length of a slug unit ls is calculated by using fs 18 and fs 20D lsu The RHS of eqn is always positive Being in divergence form means that it should change sign as the direction of iteration changes between topdown and bottomup In his thesis Hadgu appears to ignore set to zero the momentum flux term in the annular flow region Since our experience with simulations is that momentum flux can be very significant in the upper parts of a geothermal well we use the Baroczy correlation he favours for annular flow inverted in equations 64 to 66 to evaluate the usual momentum flux G lt v gt F 5 Duns and Ros Correlation The Duns and Ros correlation 14 is based on the Ros correlation developed from tests on a vertical pipe rigged in the laboratory and using air water and air liquid hydrocarbon mixtures As for the Orkiszewski correlation the Duns and Ros 1963 correlation 14 starts with known values for the superficial velocities Uga Sty QXfow PyA and Us 1 Su Q 1 Xgow p1 A and with wellbore diameter and liquid viscosity values for dimensionless liquid velocity dimensionless gas velocity dimensionless pipe diameter and dimensionless liquid visco
27. in Fig note that the transition regions are small but necessary to ensure a smooth and monotonic dependence of saturation on 6 Note that there is no attempt in the Orkiszewski correlations to account for the deviation of the well from vertical The correlations are vertical ones Subsequent studies of inclined flows suggest this is OK for angles up to 20 or 30 from vertical Further inclination towards horizontal tends in practice to remove the bubble regime as bubbles tend to hit the upper casing surface and coalesce and this is not reflected in the Orkiszewski correlation Counterflow where vapour flows upwards and liquid downwards as in a heatpipe setup is not accounted for in this formulation either except in the slug regime since it would give negative 3 values and flowing enthalpy dryness would depend on which end of the well was being considered 69 DIRECTION OF FLOW DIRECTION OF FLOW gt ANNULAR SLUG ANN LAR MIST TRANSITION BUBBLE SLUG Figure 22 An illustration of the four flow regimes in the Orkiszewski correlation taken from his paper 29 Figure 23 An illustration of the Orkiszewski correlation showing saturation S as a function of volumetric gas flow 3 for a flow of 20 kg s and fluid properties corresponding to pressure 23 bara and temperature 215 C The four red disks show the places where bubble changes to bubble slug then to slug then to transition and finally to mist as
28. in GWELL and is based on the lecture notes by Kjaran and Eliasson 22 The choke velocity is estimated as km UCH os Pm where km is a two phase fluid bulk modulus defined by km ko ky 1 S 1 8 where each individual phase bulk modulus for subscripts v for vapor and I for liquid is given by If the two phase mean velocity Su 1 S u exceeds the choke velocity ucy the user is warned in the last column of the main output file with the word CHOKED at the depths that choking occurs If the two phase mean velocity exceeds twice the choke velocity the simulation is stopped and the user is warned that choking occurred Choking is more commonly observed in bottom up simulations The subroutine that checks for choking is a separate calculation included as an extra check on the solution obtained I expect that solving the conservation equations becomes difficult when choking is encountered giving non convergence of the Newton method separately to the choking check being discussed here 103 H Non condensable Gases The effect of non condensable gas in geothermal fluid is to alter the saturation curve for pure water due to the effects of partial pressure of the non condensable gases The most prevalent such gas is carbon dioxide which is used in SwelFlo as a proxy for all such gases In particular if the temperature of liquid phase liquid is raised it begins to boil at a temperature the bubble point that is below th
29. in a number of ways We note here some alternative forms for the momentum flux based on rewriting CuvSpyu Call S pu Gum Gv where Um is a flowing average velocity Um Cov X Uy Call uy and X is the flowing steam quality the steam mass flow rate divided by the total mass flowrate Ms _ AS pyUy Q Q The second form uses the effective specific volume CwX KCal X 1 X A ee Pu Pl X K where K is the slip Ul It is useful to note that Uy X 1 5 Pl when deriving the above alternative forms for part of the momentum flux Then the momentum flux term can also be written in the forms Lig ae 118 since AG Q is a constant or using the second form above see also eqn 14 E 79 AG Ue or La TA Jz QGv gt or bringing the constant Q outside the derivative o G5 Gue Then conservation of momentum is often written as an expression for the pressure gradient down a geothermal well as the sum of the divergence of the momentum flux the frictional losses and the gravitational or static pressure gradient Op 1 2f Qum lt pulvl gt pmgsin 112 az 1g Q q lt Pl gt Pmg 112 and since AG Q is the steady mass flowrate Q can be brought outside the derivative in the momentum flux term so that momentum conservation can be written as Op o 2f n Go Um 5 pvlv gt pmgsin 113 Chisholm s version of the momentum flu
30. in the current folder The depth if any at which liquid flashes to two phase is displayed just below the wellhead values In both topdown and bottom up modes the simulator detects changes in phase between liquid flow two phase flow and vapour flow In particular the depth at which a change from up flowing liquid to up flowing two phase fluid occurs the flashpoint is located to within the minimum mesh size accuracy of 2 cm by subdividing nodes when such a change is detected until the minimum node size is reached This helps to reduce the discretization effects of coarse meshing on output curve computations SETUP DATA INPUT RESULTS GRAPHS Wellhead P Bara Q kg s F flowing H kJ kg in place H kJ kg flowing X in place X XCo2 0 800000E 01 20 000 169 75 920 00 837 56 0 10060 0 59600E 01 0 00100 Flash Depth 704 320 Feedpoints Note if no feeds are input will at least show bottomhole here Z mMD Pwell Bara Q kg s Tres Hres kJ X XCo2 Type Pres Bara P1 kg s M 1000 0 48 505 20 000 1692 e2s8se 0 000 0 100E 02 prod indx 85 000 30 792 Figure 14 A view of the RESULTS tab in the main SwelFlo window In the detailed output file are any downhole values computed depth pressure temperature static dryness enthalpies of liquid and vapour phases velocities of liquid and vapour phases densities of liquid and vapour phases total static fluid enthalpy wellbore r
31. just the FULL filename with full path information The file type is an integer 1 for pressures 2 for temperatures 3 for velocities interpreted from spinner data 4 for output curve data There can be only one type of data in a single data file but there can be many pressure data files for example The units for depths pressures and flowrates in all data files are taken to be as set by the user in the simulator menu and data files are assumed to contain data in the currently chosen units line 7 data file name2 character type of file integer The name of the next data file to be read and type line 5 NDATAF data file name character type of file integer The name of the next data file to be read and type line 6 NDATAF NOTOPS integer The number of lines of wellhead values provided next Exactly 2 lines are assumed present and the second one is ignored if important terms are out of range line 7 NDATAF P Q T H Hflo X Xflo XCO2 floating point P is wellhead pressure Pa Q is flow at the wellhead positive for production negative for injection kg s 43 T is temperature H is static enthalpy and Hflo is the flowing enthalpy of the fluid at the wellhead X is the static steam mass fraction dryness Xflo is the flowing steam quality and XCO2 is the mass fraction in the total fluid of CO2 at the wellhead These values determine the starting conditions for topdown simulations line 8 NDA
32. os is gogn unu gtp report UNU GTP 1983 02 pdf C T Kelley Iterative Methods for Linear and Nonlinear Equations SIAM 1995 Malanin S D 1959 The system water carbon dioxide at high temperatures and pressures Geochemistry No 3 pp 292 306 M J McGuinness The impact of varying wellbore area on flow simulation in the Proceedings of the Thirty Eighth Workshop on Geothermal Reservoir Engineering Stanford University Stanford California February 11 13 2013 SGP TR 198 Six pages M J McGuinness Feedpoint viscosity in geothermal wellbore simulation Geothermics 50 April 2014 24 29 DOI 10 1016 j geothermics 2013 07 007 McGuinness M J 2013 Correct Energy Conservation in Geothermal Wellbore Simulation Submitted to Geothermics 2013 minor revisions requested Nield D A amp Bejan A Convection in Porous Media Springer Verlag NY 1999 Orkiszewski J Predicting Two phase Pressure Drops in Vertical Pipe J Petroleum Tech June 1967 pp 829 838 Pritchett J W Rice M H amp Riney T D Equation of state for water carbon dioxide mixtures implications for Baca Reservoir DOE ET 27163 8 1981 Serghides T K Estimate Friction Factor Accurately Chem Eng March 5 1984 91 pp 63 64 32 Sutton F M Pressure temperature curves for a two phase mixture of water and carbon dioxide N Z Journal of Science 1976 19 pp 297 301 33 Takenouchi S amp Kennedy G C The binary system H2O C O a
33. same and the momentum flux term is related to total pressure gradient as dp Usa Us Uscom dp dp TP MF 114 2 P dz dz did HA acc where P is pressure Pa and 2 is the total pressure gradient Since the total pressure gradient term in the momentum conservation equation does not depend explicitly on A z this method used to relate momentum flux gradient only to overall gradient cannot be used when area depends on z since then momentum flux depends on z through both pressure and area That is the full Duns and Ros correlation cannot be upgraded in any obvious way to allow for area varying with depth for mist flow The Hadgu correlation takes the usual form o for bubble and mist regimes in the original constant area formulation so I have modified it in the same way as above that is I have altered it to read o However in the slug and churn regimes a special correlation is used to correlate directly the gradient of the momentum flux which depends on the velocity of the down flowing fluid around the Taylor bubble following the Fernandes et al development The correlation depends on the wellbore diameter and does not adapt to conform to the variable A z derivation here Hence no attempt is made to do so for these two regimes under the Hadgu correlation It is uncommon to see area A included as a variable function of depth in the geothermal wellbore simulation literature A is usually assumed to b
34. setup between several different projects Copying and pasting between cells of a grid window and or text files is supported 4 6 Reservoir Temperatures Reservoir temperatures versus depth may optionally be entered into this grid window Either the distance along the wellbore or the true vertical depth of a reservoir temperature may be entered The simulator expects one or the other of these columns to be used but not bits of both If both are used then either the column with more values is used or if both columns have the same number of values the true vertical depth column is used The reservoir temperature profile entered here if any is used to calculate transient heat loss or gain between casing and country The mathematics of that calculation may be found in the Appendices The Edit Row menu item in this grid window allows the user to insert or delete entire rows or clear their contents A maximum of twenty rows is currently allowed The Clear menu item allows the user to clear the contents of the current cell the selected cells or the current row It is possible to Save the Reservoir Temperature profile to a file that may be Opened by another project from its Reservoir Properties grid window over riding any pre existing reservoir temperature profile just as for the Wellhead grid window The default name for this file is SWFTres csv 20 Pe RESERVOIR TEMPERATURES o E x File Edit Search Save SaveAs EditRow Clear DE aa
35. the simulator inserts a new node at the feed depth to make this so Feeds must be in depth order with the shallowest at the top The following checks are made 1 if DELZ does not divide evenly within 1E 07 into SECL it is altered so that with the nearest integer number of sections it does 23 if LENGTH does not match the sum of section lengths it is reset to equal the sum of section lengths the user is warned in the error and progress windows and the new value of LENGTH is displayed in the main input window If bottom feed depth does not match LENGTH of well within 0 01 length is set to bottom feed depth However for a topdown simulation the bottom feed is allowed to be absent from the input information If so it is assumed to be at the well bottom a distance LENGTH from the wellhead and information about this bottom feed flowrate enthalpy wellbore pressure etc is available in the detailed output file Show Detailed Output but only those feeds whose information is provided as input to the simulator have their information summarised in the main RESULTS tab If bottom feed depth does not match sum of section lengths within 0 01 simulation stops with detailed warnings If any feed depth does not match up with a subsection DELZ within 0 1 simulator inserts a new node to match the feed depth This mismatch might be caused by previous adjustments and is best fixed by altering the input file Feedpoints must
36. the previous section and coded into Fortran Note that there are ranges of pressure gradient especially for lower pressure values and higher flowrates where there are two values of Uga for one value of pressure gradient That is the correlation is sometimes noninvertible given a pressure gradient there may be two values of Usa matching it This corresponds to the observation made by Hagedorn and Brown 21 p 480 that sometimes there is a reversal in the curvature of the pressure depth traverses near the top of the well and that the correlation for 4 matches this reversal with a high degree of accuracy Figure shows more detail near zero Usg and for a slightly different value of S and illustrates that S does go to zero as Usg goes to zero and that the pressure gradient goes closer to liquid static which dominates the friction gradient in contrast to what happens near pure vapour where friction gradient dominates Figure shows computed total pressure gradient versus Usg together with the values of friction pressure gradient gravitational pressure gradient pseudo saturation S and homogeneous saturation It is clear from the figure that the pseudo saturation is in a one to one relationship with Us It also clear from this figure that the interaction between an increasing friction pressure gradient with increasing Usg and a decreasing gravitational pressure gradient is the cause of the noninvertibility of the overall pressure gra
37. trying to conserve energy and momentum near here This is smoothed by inserting a bubble slug transition in the range S 0 28 0 32 Usz is set to a linear function of S between these values moving smoothly from bubble to slug regions The jump from negative Us values in slug flow to positive values in mist flow is more problematic and occurs at physical S values for very low flow rates The condition that the mist value of Ug gt 1 0 is sometimes not met until saturation reaches values where slug values of Us become negative or sometimes is not met for any real values of saturation This and smaller jumps seen at higher flowrates are removed by introducing another linear transition region called transition between slug churn and mist flow The right hand end of the transition region is given by the mist flow values for which U 1 0 say S and U The left hand end of the transition region we choose to have a slug flow value Usz that is a little greater than Ug that is that Ug Ug where 0 01 Then the forward problem is solved to find the corresponding value S for S The transition region is the straight line connecting S Ug 6 to S U For very small flowrates it is necessary to check that Ug is not greater than the maximum possibly value G p and if so to reset Us to G p 76 2 57 1 54 Usl 0 54 0 0 2 0 4 0 6 0 8 1 S Figure 26 An example of Ug versus saturation fo
38. usually much too large and many of the simulations fail to complete since there is no steady state solution Refining this range speeds up the computations and helps provide a good selection of output curve points to be plotted If a range of flowrates is specified only those results lying in this range will be plotted in the output curve In topdown mode if feed information is provided at bottomhole optionally the reservoir pressure may be used to determine a productivity index given the calculated flow calculated wellbore pressure and effective viscosity at bottomhole Then if an output curve is requested together with a topdown run a series of bottomup simulations are run to produce an output curve using this productivity index over riding any nonzero value that the user has previously set at bottomhole The value of bottomhole productivity index and enthalpy shown in the input window is as set by the user the values of feedpoint enthalpy and productivity index displayed under the RESULTS tab are the values computed in the topdown run and used to get the output curve values detailed under the output button in the RESULTS tab that is in the main output file corresponding to the output curve results are also the actual values used as computed from the topdown run A productivity index is also calculated for any other types of feedpoints for which 1 one is not already provided and 2 it is not a dry steam power law case and 3
39. viscosity m s and flowrate and pressures are now in kg s and Pa This productivity index is related to the more usual productivity P by after converting P to SI units Pyv 3 6 x 105 The best way to calculate an output curve is to use equation to calculate the flow Q at a feedpoint given the pressure in the wellbore since this properly accounts for the effects of changes in viscosity on feed flowrate the productivity index is independent of these while the usual productivity does depend on the current value of viscosity The simulator will calculate a value for productivity index at bottomhole for all feedpoint cases except that of dry steam drawdown If an output curve is requested this computed value will be used unless a value has been provided explicitly in the input file and the type of feedpoint is type one see the later section on input file setup This formulation is straightforward to apply when v is roughly the same in the reservoir and in the wellbore When v depends strongly on pressure as it does when flow becomes two phase some care is needed in calculating an effective value of v as it may vary considerably between reservoir and wellbore due particularly to changes in fluid density Carelessly calculated viscosity values can show up as anomalous drops in the flowrate as pressure at bottomhole is decreased through the flashpoint giving unusual looking output curves 121 Equation 115 still hold
40. 2 Otherwise 78 If S lt S transition flow is used to find Ug Note that this will not arise if S lt 0 32 Otherwise Mist flow is used to find Us Inverse Problem An iterative method is used to find phase velocities and hence saturation given flowing steam quality Xgow The Newton method modified by Armijo s Rule was attempted but the transition regions between regimes can confuse the Newton search so Brent s method is used This iteration repeatedly calls the forward problem described above Given Xgow and P T Q the iteration actually returns the values of u and u Then saturation S is determined by equations 1 and 2 as QO Xow UyApy F 4 Frictional Pressure Drop Full Hadgu Hadgu also uses these flow regimes to calculate different forms for the friction term in the momentum equation These are incorporated into SwelF lo as detailed here To use Hadgu s momentum terms as well as Hadgu s phase velocities it is necessary to choose the Hadgu full correlation in the Input Tab options Hadgu s thesis still has a slightly different treatment to that of SwelFlo of heat flow between wellbore and surrounding reservoir Different heat transfer coefficients are allowed for the different regimes So some small differences might still be expected between SwelFlo with the full Hadgu correlation chosen and another simulator based on Hadgu s thesis F 4 1 Bubble Flow In this regime described in th
41. 25 but since we start out by knowing a it is simpler to check whether it is greater than 0 3 or not F 3 3 Slug Flow Slug flow is considered to hold for S between 0 3 and 0 52 Hadgu considers the rise velocity of a Taylor bubble with wall drag near the bubble being the dominant factor and uses the treatment of Fernandes et al 16 developed from a combination of analytical work and experimental data Many of the variables which appear in this subsection are illustrated in Fig 24 In summary the rise velocity of a Taylor bubble in an infinite medium is calculated as Par 0 5 Pl 42 Uso 0 345 where D is the casing inner diameter the translational velocity Uy of a Taylor bubble is related to the superficial vapour and liquid velocities by Un 1 29Usa Ust Uco the rise velocity of bubbles in the liquid slug is as was used in the bubble phase also gApo as 0 5 Up 1 53 1 ALS oS pi where azs 0 3 the velocity of the liquid film around the Taylor bubble is Urrg 9 916 gD 1 yarg 9 9498yr 72 TAYLOR BUBBLE FALLING FILM 20 Po o9 O Go LIQUID SLUG O O E c 000 bo O Figure 24 A slug unit taken from Hadgu s PhD thesis his Fig 7 5 where D is wellbore inner diameter the average velocity of the liquid in the slug is l a Uns Uy Uy Urre gt ALS the average velocity of the bubbles in the liquid slug is UcLs Urrs Uo the average
42. 4 a 1 4 y gt Figure 29 The values of Fi F3 F3 and Fy used in the bubble region and their dependency on dimensionless liquid viscosity Nj according to Duns and Ros From the uppermost curve downwards they are F3 F F4 and Fy F 5 2 Slug Flow Region The dimensionless slip velocity for slug flow under the Duns and Ros correlation is given by the relationship N0982 7 Ev Vsp 1 Fs 1 Fy Nw oa where F 0 029N Fe and where the parameters F Fg and Fy depend on dimensionless liquid viscosity as illustrated in Fig and listed in Table 5 F 5 3 Mist Flow Region Slip velocity is assumed negligible in mist flow so that Vsp 0 and as in eqn 96 S Use Use Us Note that zero slip flow implies that AAA A D UDE IDAS ON S SA Ue me ee Gaa A R Daa AE E E R E Se ee man EE IS eT ee bol LL pap SSE ma eo ware mati o i ae RE A ma el E E rt ep A bh Figure 30 The values of F5 Fs and Fy used in the slug region and their dependency on dimensionless liquid viscosity Nj according to Duns and Ros F 5 4 Transition Region In the transition region between slug and mist flow regions Duns and Ros calculate a total pressure gradient that is the linearly weighted average of the mist and slug total pressure gradients using N as the measure of distance Since our focus is on computing slip Vsp and since transition is relatively small and not often encountered we will instea
43. Az 49 A Bargl A A243 aLsAr 50 1 arg Az 51 s oe 51 LD Le 52 QTB A 1 A3 53 B us 54 QTB ALS Xp 55 280 Xt Xa dd q G 4 56 Urre 9 916 gD 1 varg 57 gApo 0 5 Us 1 53 E 1 ars 58 i DAA o Us 0 345 a e 59 Pl QTB 0 89 60 ALS 0 3 61 62 and then the gas velocity can be calculated from G Usa Us 74 F 3 4 Churn flow The churn flow region is a relatively narrow one sometimes called transition The boundary between slug and churn is given by a 0 52 and 0 4 E 93 0113 5 0 092 2D Mid 2 5 Apg o D Xp which can be solved to find a value for Ur where Ur Usa Us The boundary between dispersed bubble and churn is simply a 0 52 Hadgu has mixed notations here with ay being a flowing void fraction which I call 8 consistent with 6 while 0 52 in 6 is the usual static void fraction Yet Hadgu sets ay a in his criteria for transition Us 0 725 5 15085 This does not matter much since no correlations are given in Hadgu s thesis for churn flow He uses slug flow parameters for pressure drop calculations so this is what we use for calculating relative speeds of liquid and vapour we use the slug flow correlations F 3 5 Annular or mist flow The transition from churn to annular or mist flow is given by U gt 1 0 where 0 5 Pv In the annular flow regime Hadgu ch
44. Bara 0 000000E 00 Total well length m 1000 00 El Pre Process Figure 7 Screen snap of the SETUP tab in the main SwelFlo window Item Max Number Separately defined sections of wellbore 500 Number of subdivisions of the wellbore 5000 Feedpoints 10 Points where reservoir temperature is specified 20 Data files 20 Data in each file 10000 length of file names 50 Points in simulated output curve 400 Table 2 List of maximum dimensions for variables in SwelF lo 16 4 2 Topdown and Bottomup modes The simulator will run either in topdown or in bottom up mode chosen in the Downhole Profile section of the SETUP tab Either production positive total flowrates or injection negative total flowrates may be simulated in either mode Topdown means that the wellhead fluid flow properties are used to initiate the solution of the steady state conservation equations and the simulator steps from the wellhead down the wellbore adding or subtracting fluid at intermediate feeds on the way down until at bottomhole there is a final set of values of flowrate enthalpy carbon dioxide and pressure computed for the fluid flow there Each bottom up simulation begins with some value for wellbore pressure at bottomhole as well as other properties like flowrate mass fraction of carbon dioxide and temperature or enthalpy or dryness and iterates upwards to solve the conservation equations adding or subtracting fluid at feeds on t
45. Figure 4 SwelFlo first steps showing simulated pressure temperature and spinner velocity profiles from a bottom up run matching values saved from a topdown run 12 using the resulting values to solve back in the other direction and seeing that the numerical solutions match to within graphical accuracy provides direct evidence of that accuracy More commonly the correct bottomhole pressure that matches wellhead values may not be known a priori In this case the user can set a range of bottomhole pressures to be automatically searched upon or let the simulator choose a range by setting the bottomhole P range values both to 0 0 Try this set both to zero and RUN Note that the search for the correct bottomhole pressure takes a lot more time many simulations are run until a match is found to the wellhead values provided The resulting best bottomup run should be indistinguishable from the previous one Note also that if the feed enthalpies are not carefully chosen there may not be an ideal bottomhole pressure that gives a match to wellhead values The process this tutorial has just gone through of performing a topdown run using wellhead values to infer bottomhole properties then performing a bottomup run using inferred values is an excellent process for setting up wellbore properties given wellhead values It also sets up feedpoints for computation of an output curve which does require bottomup runs with varying bottomhole pressures
46. Gas For the gas mixture the same expression used by Sutton 1976 and Pritchett et al 1981 is used Py P Ss Pco2 5 17 where Ps the density of steam Pcoz the density of CO Both values of density are evaluated at the given temperature and corresponding partial pressures 5 1 4 ENTHALPY 5 1 4 1 Carbon Dioxide CO The enthalpy of CO is evaluated using the formula given by Sweigert et al 1946 4 135E 04 H 1688 1 542T 794 8L0G Ty co2 3 571E 04 Poo oh 1 7 576E 08Po02 5 18 T 100 3 co2 where the enthalpy of CO in kJ kg the partial pressure of CO in Pa the temperature in K co2 P co2 J ka tl 26 5 1 4 2 Heat of Solution The heat of solution of a gas is the change in enthalpy brought about by the dissolution of the gas in water The equation of the heat of solution is calculated using a polynomial fit to the experimental data obtained by Ellis and Golding 1963 and is valid for temperatures in the range 100 300 C Holin 771 33 6 0198T 0 07438T 2 9244E 04T 4 4522E 077 5 19 where H oln the heat of solution in kJ kg T temperature in C 5 1 4 3 Enthalpy of the Mixture The enthalpy of the mixture is evaluated using Ha XH 1 x H 5 20 where x mass fraction of gas phase Hy enthalpy of the gas phase in kJ kg Hy enthalpy of the liquid phase in kJ kg The liquid and gas phase enthalpies are evaluated as
47. Law formulation 44 YC XC XC id 44 YC 18 1 YC where P vo CO2 P or the user can choose a combination that lies between these two options XC rXC AS pos by setting the parameter r to a value between zero and one on a slider control 108 For pure vapor both choices have XC XCO Then for the GWELL mass fraction the partial pressure of carbon dioxide is given by P XC P For the Dalton s Law formulation the partial pressure of carbon dioxide is given by P YC P The combination option uses P rYC 1 r XC P Given the careful comparisons made by Pritchett et al and by Sutton we recommend using the mass fraction option not the Dalton s Law or the combination option H 4 Computation of the Phase State The logic used to decide the phase state is the same as that used by Pritchett et al and in GWELL The approach is first to check if P lt Ps in which case everything is in the vapor phase Otherwise solubility XC is checked at the partial pressure of carbon dioxide Poor P Ps If solubility exceeds the total mass fraction of carbon dioxide present the phase is all liquid Otherwise the mass fraction XC in the vapor phase is computed as in the previous subsection and the formula XCO XC lt 1 a XC XO am is used to determine a steam dryness X This dryness is in turn used to check again finally that the state really is two phase or if X is effectively z
48. Presently Swelflo is designed to run under Windows XP Vista 7 and the exe file is about 10MB in size SwelFlo may run under other Windows operating systems but is not intended to It is recommended to run it on a computer or laptop with at least 1Gb of RAM Topdown simulations are almost instantaneous while output curves and bottomup runs which use an output curve to help match wellhead pressure take about 10 15 seconds to run under Parallels on a 2010 MacBook 2 26 GHz Intel Core 2 Duo 3 First Steps SwelFlo is usually installed in its own directory at the installer s choice When SwelFlo is first started 1t will ask the user for a project directory to work in This directory will be used to save input and output files Subsequent startups of SwelFlo will open the last project directory that was used on that computer Once SwelFlo is installed start it by double clicking on the exe file icon or on the desktop shortcut If necessary ensure you are starting it as administrator In Windows 7 this is done by selecting the exe file then choosing Organise then Properties then Compatibility then tick the run as administrator box Or you can right click the exe file and choose the run as administrator option You will be prompted to choose a folder to work in which should be somewhere in your documents folder You can create subfolders as part of the selection process If SwelFlo has not been run in the working folder befor
49. SwelFlo User Manual Version 2 02 Dr Mark J McGuinness Marsan Consulting Ltd Wellington NZ December 2015 Abstract This is a User Manual for the geothermal wellbore simulator SwelFlo The graphical user interface to the simulator is compiled to run under Windows XP Vista and Windows 7 and is designed to give easy and fast setup and to allow quick plotting of data together with simulated results while retaining functionality It is intended as a tool to model flow in a geothermal well to be used by reservoir engineers and modellers who wish to run simulations predicting the effects of changing the setup and operability of a geothermal well or the effects of changes in reservoir properties over time on production from a geothermal well Both production and injection flows may be modelled Flow of vapor and liquid phases of water together or alone together with carbon dioxide is modelled Simulations may proceed from the bottom up or from the top down with corresponding differences in which boundary condition is important Downhole values of variables may be examined for a simulation or multiple bottom up runs may be combined to produce an output curve Multiple feedpoints are allowed with variable wellbore geometry heat flow between wellbore and the surrounding country and production or injection Mark McGuinness 2015 Contents 1 Version notes current version 2 Introduction 2 1 System Requirements kad 6
50. TAF P Q T H Hflo X Xflo XCO2 floating point The second line of wellhead values if another simulation is wanted line 7 NDATAF NOTOPS WHPTOLERANCE PMINSET PMAXSET floating point values WHPTOLERANCE is the desired relative tolerance used to decide when wellhead pressure is matched for bottom up simulations Recommended 1073 If this is set too small the simulation will still find a nearby WHP as best it can PMINSET PMAXSET are the maximum and minimum allowed values for bottom hole pressures to use when computing an output curve or to restrict the range when searching for a bottomup match to desired wellhead conditions line 8 NDATAF NOTOPS LENGTH a floating point number the total length of the well mMD metres measured depth not vertical depth but is distance down the casing If it does not match the depth of the bottom feed it is reset internally to that depth Then if it does not match the result of adding up all of the section lengths the program stops It is used exclusively together with angles to the vertical to calculate the changes in potential energy of the flow line 9 NDATAF NOTOPS THCON floating point number the thermal conductivity of rock usually 2 0 W m deg C line 10 NDATAF NOTOPS RHOR floating point number the density of rock usually 2800 0 kg cubic m line 114NDATAF NOTOPS HCAP floating point number the heat capacity of rock usually 1000 0 J kg line 12 NDATAF NOTOPS TIME
51. The simulator will attempt to further subdivide a subsection where there is no convergence of solution to the conservation equations down to a size of about two centimetres This happens particularly if a flash point is encountered It is recommended that the user set sizes of the order of five to twenty metres for subdivision lengths mesh size Casing Identifier values will override any radius values entered and any roughness values entered Preprocessing will cause the over riding values to appear in the CASING GEOMETRY window radius roughness and angle If radius is set to zero or roughness to a negative value they are over written by the default values radius 0 1m and 21 roughness 0m The Casing Identifier drop down menu is read from the file called Casing_Liner_Database csv in the directory that the exe file is stored in A user may add lines to this database as desired using a text editor in the Winows environment The format of this csv file is that the first line is ignored so can be a header and the following lines take the form Casing ID text which can include the character inner diameter inches outer diameter inches and roughness inches All of these should be separated only by a comma Quote characters should not be used to surround each entry the commas are sufficient to separate values Here is a sample from the first four lines of this file CASINGLIN INNERDIAM OUTERDIAM ROUGHNESS CASING BASIC 0
52. Tk R Tk where Peo density of CO in kg m R the gas constant 1 88919E6 erg g K Tx temperature in K Pr pressure in bars z Ph Tx gas compressibility factor evaluated using an analytical fit of the data by Vargaftik 1975 For pressures less than 300 bars Z PprTx A B P 300 C P 300 D P 300 E P 300 4 5 8 24 For pressures greater than 300 bars Z P Tx A B P 300 F P 300 5 9 where Po the pressure in bars Tk the temperature in K The temperature dependent coefficients are evaluated from 2 3 4 A A A T aaTK aah ATE 5 10 B By BT BaT B T2 BaTk 5 11 3 4 4 D Dy D T DT D T DATE 5 13 1 A 300B 3002c 300 D E 7 5 14 300 2 3 4 F Fo F T F T FIT FIT 5 15 The values of the coefficients given in Table 5 2 give a satisfactory fit to the exper imental data between 77 to 350 C TABLE 5 2 VALUES OF COEFFICIENTS FOR CALCULATION OF CO DENSITY 1 72976E 10 1 06781E 12 5 36712E 14 Be 05175E 16 95665E 15 1 3 2 Mixtures 1 3 2 1 Liquid Since the amount of CO dissolved in the liquid phase is small low solubility of CO in water it is assumed that the density of the liquid is equal to the density of pure water at the same temperature So 25 Pi Pio 5 16 where P the density of mixture Puzo the density of pure water 5 1 3 2 2
53. a velocity file generated in the current project folder The names are set to PDataSim001 tzt TDataSim001 txt and VDataSim001 txt with numbers that match the corresponding SWFoutput file The format of these files is the same as the required format for data input files a first line with Q whp and H flowing enthalpy at wellhead subsequent lines with two columns the first being depth and the second being the variable all in units chosen by the user Then a simulation that is subsequently run in another project folder or in the same project folder can use these data files and display them for comparison against perhaps differently set up simulated results The spacing in m between depth values chosen for sampling the simulated results may be set by the user under the INPUT SETUP tab in the box labelled Spacing for data output m 5 ProgressReport tat in the project folder holds the text that was most recently displayed in SwelFlo s Progress grid window Each time the simulator is restarted or visits a project folder from elsewhere this file is overwritten 6 SGWlog txt which is a general purpose text file for logging purposes and is stored in the same directory as the exe file It is overwritten whenever SwelFlo starts 7 ErrorReports txt is the file where error reports are kept in the same directory as the SwelFlo exe file It is overwritten whenever the simulator exe is started up 10 Command Line Execution SwelFlo can no
54. adius flow regime mass concentration of CO2 total mass flowrate flowing enthalpy flowing steam quality that match most closely the wellhead pressure and flowrate given depth of flash point if a transition from liquid to two phase occurs in the wellbore and a summary of feed and wellhead values for successful output curve simulations 7 Graphs Tab Plots may be viewed under the GRAPHS tab Plots of simulated results are produced combined with any data provided for immediate feedback on how closely the simulated 31 wellbore setup matches available data Presently pressure temperature velocity and if requested output curve plots are presented The pressure plot shows total pressure solid line and saturation pressure of pure water Psat dashed line together with any pressure data provided The temperature plot shows the fluid temperature versus depth together with any temperature data provided The velocity plot shows the simulated mass weighted velocity together with any spinner velocity data provided The mass weighted velocity is the combination Xu 1 X u 17 where X is static steam dryness u is vapor phase velocity and w is liquid phase velocity These plots are only intended to provide quick feedback so cannot be altered or adjusted Separate csv files are output to the project folder that allow the user to tailor plots for presentation purposes using an application like Grapher or Excel The units for th
55. adius to 0 1m 66 where Lae 1 8 1 2 which as illustrated in Fig is a very close match to the detailed Table results over a wide range of flowrates and wellbore conditions for 6 lt 0 9 It is also seen to provide a good match to Armand s data Since this excellent fit also behaves smoothly near 8 0 9 this is what is adopted in SwelFlo the simple formula 25 F 2 Orkiszewski Correlation The Orkiszewski correlation provides a different way to get saturation S given P or Afi Orkiszewski s correlation starts out assuming that P and temperature is known and from that value it finds liquid and vapour phase velocities Then saturation can be computed S 6 is given by the following Orkiszewski procedure 1 start with a value for Xgoy and we also know the temperature and pressure 2 compute XfowVG _ 26 XfowVG 1 X flow UL 1 Xpow Vz j XfowVG XfowG pi iji _ PL 28 oe Pu C Pa G Ni b n 6 fov 1 Xa E 29 pi LCP pv Po pl 0 7275v2 De max L071 0 13 30 Ls 50 0 36 0u p0 31 Lu 75 0 84 0 0 p0 J2 32 0 25 Us 153 eee 33 Pi Us 0 35 9D 1 po pi 34 where o is surface tension D is wellbore diameter vz is specific volume one over density of liquid and vg is specific volume of the gaseous phase 67 3 Then if 6 lt Lp the regime is bubble flow set u vr Ugi and GX fow ga eii 2 Uvpu 4 Oth
56. al arrays and the simulator is now based in the newly chosen Project directory The user can select a new path involving multiple new folders and the simulator will create each new folder as required to create the new path 8 1 2 Import Existing SwelFlo setup This item allows the user to select and import the setup from another existing SwelFlo directory by reading its SWFinput file without changing the current working directory Project folder The new setup is written to the current and unchanged Project directory over writing any previously set up simulation state in the current directory Hence the user is warned before any action is taken here 8 1 3 Load Template This menu item loads the internal template detailed in Section 3 1 for a topdown simulation of the downhole profile of a producing vertical well of length 1000m with wellhead pressure P 8 bara flowrate Q 20 kg s temperature T 169 75 C and mass fraction of carbon dioxide XCO 0 001 providing a simple ready to run simulation that may be used as a starting point for modification 8 1 4 Switch to Existing Project Selecting Switch to Existing Project also changes the current working project directory to the one containing the selected input file effectively switching to another existing project Usually the user will select the file named SWFinput txt which is the file used by SwelFlo to store the current setup and is the file read by default when SwelFlo visi
57. alue Ny Ly plus a small amount A Then if S is less than S A we need to correct it to linearly interpolate between S and S A based on the current value of Usa Note that L actually depends on Ng They are equal when 0 25 0 25 A v go pl go 0 25 pla 2 42 pit Lop that is when Usa Usa We then need an explicit value Ug for Usg that corresponds to the saturation value S A 0 25 in the slug flow region by solving simultaneously setting V 2 Vsp and using eqn in the equations Use UsL S A 1 S A NO982 4 FY 1 FN Va Vsp 14 Fs Since Ug varies with Usa keeping G fixed and since also N depends on Usz this is a nonlinear relationship that needs to be solved using a Newton method to find the value Usa Then the interpolated value for S in the bubble slug transition range lies in the range S to S A and is calculated for the current value of Usa as Sa 2 A Use U a F 5 6 Low Flowrates Note that for low flowrates Q less than 15 kg s for a well of radius 0 1m even with the bubble slug transition region introduced here a plot of S against Usa reveals that the Duns and Ros correlation is not invertible there are S values with more than one Usa value 88 generating them and that some S values cannot be realized at very low flowrates This is illustrated in Fig 31 where flowrates of 1 5 10 and 15 kg s are used with a wellbor
58. and are designed to help identify a particular simulation Each line may be up to 80 characters long line 4 ANS VELT NSIMVALS WANTOUTC DOUBLESIM GWELLFLAG SIMRANOK GOTOUTC READYTOPLOT GETPROFILE ALLOWHFLO QMAXOC QMINOC NSIMOUTCPTS SPACNG These variables should be separated by blanks This line is read in as a large character max length 300 characters then parsed by looking for the blanks ANS is 11 VELT is I1 and NSIMVALS and WANTOUTC are I of any size They are integers DOUBLESIM is a real F3 1 GWELLFLAG is an integer and SIMRANOK GOTOUTC READYTOPLOT GETPROFILE ALLOWHFLO are logicals T or F QMAXOC and QMINOC are real numbers read in free format NSIMOUTCPTS is an integer and SPACNG is an integer read in free format ANS determines whether simulation is topdown ANS lt 1 or bottomup ANS gt 1 say 2 VELT determines which 2 phase flow correlation will be used Currently VELT 0 gives Armand VELT 1 gives Orkiszewski VELT 2 gives Hadgu for phase velocities only VELT 3 gives a full Hadgu correlation with Hadgu s choice of momentum terms as well as phase velocities VELT 4 gives the Duns and Ros correlation for phase velocities only VELT 5 gives the full Duns and Ros correlation including its treatment of momentum conservation and VELT 6 gives the Hagedorn and Brown correlation NSIMVALS is the number of points desired to be computed in the output curve DOUBLESIM is redundant GWELLFLAG is redundant
59. ation can be 100 8000 Chart Area 7000 6000 0 8 5000 0 6 4000 PGRAD s dP dz 3000 0 4 S homog 2000 0 2 1000 0 0 0 10 20 30 40 50 60 70 USG Figure 38 A plot of the pressure gradient obtained Pa m versus Usa value input to the Hagedorn and Brown correlation for P 40 bara S 0 44 Q 50 kg s showing more detail near zero Usa values Also shown are the pseudo saturation 5 and the homogeneous saturation which forms an upper limit on what the value of true saturation can be value from eqn 96 Su Usa Usa Usu This is unphysical if S is interpreted as actual saturation since the homogeneous saturation corresponds to equal gas and liquid phase velocities and for upflow vapor can not flow slower than the liquid phase since the vapour phase is more buoyant Hence it is necessary to set S equal to the smaller of S and Sy as illustrated in figure 39 In this figure we see a shortcoming of the correlation it has S values greater than the maximum possible homogeneous flow value that is at a given saturation the gas flow velocity is less than the liquid phase flow velocity there which is unphysical So in SwelFlo if this happens homogeneous flow is used In homogeneous flow there is no slip between gas and liquid which is also not very realistic for bubble flow Steam fractions are very small here which may help reduce the impact of the homogeneous assumptio
60. be one of four possible types with the codes 1 2 3 4 specified productivity index together with a quadratic Forcheimer term A which may be zero OT specified flowrates or dry steam power law parameters or linear and or quadratic drawdown parameters Productivity Index Feeds The extra parameters that must be provided in columns eight and nine of the FEEDPOINTS grid window for this type of feed are the productivity index Y m of the feedpoint and the nonlinear term A floating point in Forcheimer s equation for drawdown at a feedpoint which should be set to zero for linear drawdown This case means that the following Forcheimer equation relating pressure drop Pa at a feed to flowrate Q kg s is used 1 A Pg Ewe SY e n 70 00 where P 2 v is the usual productivity converted from t hr bar to kg s Pa if A is zero gt is the productivity index which is needed if drawdown information is to be independent of viscosity changes at a feedpoint for output curve simulations v m s is the effective kinematic viscosity for flow between the feed and the reservoir and A is a nonlinear Forcheimer term If A is set to zero this equation is just Darcy s law for linear drawdown at a feed Forcheimer s equation is an extension of Darcy s law to 24 situations where turbulent flow near the wellbore leads to lower flowrates and nonlinear dependence on pressure drop 28 Please see Appendix for details about ef
61. ber requested but will not usually be less than half A range of flowrates may be specified if it is desired to tidy up the appearance of the output curve A range of bottomhole pressures may be specified and often needs to be specified to obtain a reasonable number of points in the output curve The detailed RESULTS file will indicate the bottomhole pressures that were successful in producing steady state solutions and hence points on the output curve If just one or two points are achieved narrowing the range of bottomhole pressures to nearby values will improve the simulator s chances of producing a good looking output curve The order does not matter for the range of flowrates or bottomhole pressures Heat flow to or from the country may be turned off or on and rock properties set The time since flow began box reflects the use of a transient solution to the diffusive temperature equation and should relate to the time that the well has been constantly discharging stabilization time if comparing with measured downhole data It is possible to STOP a simulation if desired during a simulation while a progress bar is showing if computation of an output curve or a bottomup match to wellhead values is taking too long It is possible to change the total length of the well arc length which should match the total casing length An incorrect length entered here will be over written by the sum of the casing lengths on preprocessing 10 T
62. ble Flow aoaaa a 79 F 4 2 Slug and Churn Flow gues doo eek a eS 80 F 4 3 Annular or Mist Flow sopas 80 F 4 4 Gravitational Pressure Drop o o eee 81 F 4 5 Momentum Flux Termal F 5 Duns and Ros Correlation 20002 200 0 ee ee Beh ASD a Pe ee ene Se a ee Ge da Laka eeee enna en ee ee ee we Ee A oS ee ee Gee Se oes bee ek oe a Ea ee a ee ee a Pee ee bee Pee Pee be ee ee ee F 5 6 Low Flowratesl aaa a a nie deny ue Ge eet ae pee ee oe a ee a e ee ee A F 5 9 Momentum Flux Terms 4 24 54 644 64 Red ee Ee eee ee F 5 10 Gravitational Pressure Drop 2 2 24 405 6 5 be a eee ee Laos pa oe Ob et me Down HG A oe Lea died vik ee bate Vee bee ee Oe Oe BH ae ed F 7 2 Implementing the Hagedorn and Brown correlation G Choking H_Non condensable Gases H 1 Phase states pee a os es es be oe bee be oes ee a eee pees ad oo Eas Se tani od UE pln aes Gul neces pons ea ek eon nine ok ied nes E ed errada aos rr II hawk ene eee cena ee ess AS I Corrected Area Treatment un Productivity Index at Feedpoints K A Sample Input File L GWELL coding errors M Version Notes 103 104 104 104 106 106 106 106 109 110 111 118 121 124 124 125 1 Version notes current version This manual is copyrighted to Mark McGuinness Except as provided by the New Zealand Copyright Act 1994 no part of this publication may be reproduced or stored in a retrieval
63. cessary to be very careful about the rate at which the ratio of liquid holdup to 4 approaches the value 1 as Usa gt 0 that is as Nnoia 00 since a too rapid approach to zero can lead to negative saturations at very low vapour content This can affect simulations particularly near the flash point when saturations are very small 1 0 8 R 0 6 Y 04 0 2 u 1 20 1 2 3 4 10 10 10 10 10 10 N hold Figure 33 A graph of the correlation presented here of R against Nnoia a slight modification of the Hagedorn and Brown correlation Nw P o1 106 Nhola Cn n ii 5 a x r 97 Furthermore where C exp 4 895 1 0775 In N 0 80822 In N 0 1597 In N 0 01019 In N provided that also C 0 011440529 N gt 0 4 0 0018693617 N lt 0 002 and where y 1 exp 6 6598 8 8173 In N 3 7693 In Ne 0 5359 In Nu with the limits 1 8200679697 Nsec gt 0 0857938 where the limits on 4 and C n have been refined to make them monotonic and smoother to seven decimal places at the limits The decay of y to 1 0 for small N rather than setting Y exactly to 1 0 for small N is set so that small values of saturation can be explored near the flash point The formulae presented above for Ca and w are graphed in figs and 35 and look the same as the graphs used by to deduce their formulae for the Hagedorn and Bro
64. ch gives dP Zn 1 Pas Pa IG 20S la G lt v gt f lt gt gsin 8 bro ES dz A Z Jro So far this is exact for a constant area Evaluating the integrals of the gravity term and of the friction term is where the approximation is made and the trapezoidal rule applied to these integrals gives simple averages of the endpoint values Hence what is implemented in the SwelFlo code is Zn 1 f lt p gt gsin0 dz es de iffa Tar gt aP 7 lah era pial J ee fE amp The momentum flux term needs more careful treatment when area changes A simple discrete approximation to the momentum flux term is obtained by using the trapezoidal rule on G as lt p gt nti lt p gt n gsin Az NI and 56 a function of lt v gt tus tor ES i Gt dz G d lt v gt S Sii lt gt n Then the momentum conservation equation becomes in discrete form allowing for variable cross sectional area A and hence variable G di ot Gn 5 lt U gt n4i1 lt U gt n 14 tE lt P gt n lt p gt n gsin9 Az 15 lee BEN Similarly the energy equation becomes upon integrating with respect to z lt A gt 4 lt h gt Erny LE 5 ene a Q4 Az g sin 0Az a 0 where the flowing enthalpy is lt h gt Xfowhu 1 Xow ht and the flowing kinetic energy is Dee L Ao Ey i de ala 2 Uy 9 l D 1 Newton Step Given values of pressure
65. ction lengths 4 If bottom feed depth does not match sum of section lengths within 0 01 simulator stops with detailed warnings 5 If any feed depth does not match up with a subsection DELZ within 0 01 the subsection is split at the feed depth line 14 NDATAF NOTOPS NUSEC NUPO 45 integer the number of points in the reservoir at which temperatures will be provided to calculate heat loss or gain through the casing to from the country No more than 20 Can be Zero lines 15 NDATAF NOTOPS NUSEC to 14 NUSEC NUPO T_DEPTH TEMP floating point numbers separated by space s depth and reservoir temperature at that depth There should be NUPO of these lines Note that depth here as elsewhere in the input file is distance along the wellbore MD not the vertical depth in the reservoir If the bottom value is recorded as a negative value this is interpreted to mean that the absolute values of all entered depths here are in fact true vertical depths rather than distance along the wellbore next line NUFEED integer the number of feedpoints that will be modelled No more than 10 can be handled This is the last set of lines read in from the input file next NUFEED lines FDEPTH I FICODE I other values The program first reads the depth floating point and the code integer Then it backspaces and what it then tries to read again from the same line depends on the value of FICODE it just read F_DEPTH I feed depth
66. d C are the temperature dependent coefficients given above in subsection H 2 The above rearrangement is only valid if the mass based empirical law XC 4 S92 is used to relate steam mass fraction to partial pressure If Dalton s Law is used the quadratic to be solved for the partial pressure Pco2 in bars is aaPoz baPcoz ca 0 where Ad 44 B X CO Xto 1 Xto ba 44A XCO X10 18P B XCO Xgo 1 Ca 18A PX fo If a combination of the mass and Dalton s laws is chosen by the user then an iterative method is used to find pressure given flowing steam fraction and temperature from the method in subsection H 4 110 H 6 Density Viscosity Enthalpy Surface Tension of Carbon Dioxide The next six pages are lifted from the GWELL manual as the code is taken directly from GWELL code 111 where 002 Mass fraction CO in gas phase Poo2 partial pressure of CO as expressed by Equation 5 2 P the total pressure For cases of dry gas all gas state the above relation becomes P Cco2 c02 a 5 6 where co total mass fraction of CO Equations 5 5 and 5 6 above fit the experimental data better than Dalton s Law which states that the mole fraction of the component gas is proportional to its partial pressure 5 1 3 DENSITY 5 1 3 1 Carbon Dioxide CO The density of CO is calculated from the expression obtained from Pritchett et al 1981 P A A 5 7 dia z P
67. d use a linear transition from the slug flow value of Vsp to zero F 5 5 Bubble Slug Transition Numerical results indicate that the saturation S typically increases monotonically as Usa increases but not when there is a transition from bubble to slug regions where S can decrease temporarily when going from bubble to slug especially at low flowrates This makes the correlation non invertible S as a function of Ng or Usa then has no inverse function In initial versions of the simulator we started out knowing saturation and wanted to find Use and we cannot solve uniquely for Usa as a function of S when the correlation is noninvertible This issue is further complicated by the fact that the dependence of S on V is not simple as Ny increases It is observed numerically that V increases monotonically with Ng through the non invertible S region but with a jump in value there It is not clear that the dependence should be monotonic given that we start out knowing Usa and we want a saturation value But jumps in saturation prevent the Newton method from 87 solving the conservation equations so in SwelFlo we choose to avoid such jumps by making an invertible transition region between the bubble and slug regions To make the correlation invertible and smooth we introduce a Bubble Slug transition region as follows If in the slug region we check the value of S against the value S obtained by using the bubble correlation with maximal v
68. ded It is recommended that the user chooses the default input filename SWFinput txt If no input file is provided in silent mode a warning is sent to the log file and the simulator attempts to run on the last saved state in PrevSimState txt in the execution folder if any If an output path is provided it is ignored in this case References 1 Aunzo Z P 1990 GWELL a multi component multi feedzone geothermal wellbore simulator M S Thesis University of California Berkeley 250 pp 2 Battistelli A Calore C amp Pruess K 1993 A Fluid Property Module for the TOUGH2 Simulator for Saline Brines with Non Condensible Gas in Proceedings Eighteenth Workshop on Geothermal Reservoir Engineering Stanford University Stanford California Jan 26 28 1993 SGP TR 145 pp 249 259 38 3 10 11 12 13 14 15 Bilicki Z amp Kestin J 1980 FLOW IN GEOTHERMAL WELLS PART IV Transition Criteria For Two Phase Flow Patterns Report GEOFLO 6 DOE ET 27225 9 Division of Engineering Brown University Providence R I 26 pages Bilicki Z amp Kestin J 1983 Two phase flow in a vertical pipe and the phenomenon of choking Homogeneous diffusion model I Homogeneous flow models Int J Multiphase Flow Vol 9 Issue 3 pp 269 288 Aunzo Z P Bjornsson G Bodvarsson G 1991 Wellbore models GWELL GWNACL and HOLA User s Guide Lawrence Berkeley Lab University of Ca
69. dgu correlation 19 the Duns and Ros correlation 14 B6 and the Hagedorn and Brown correlation 7 Dissolved solids like NaCl are not accounted for Noncondensable gases are modelled as carbon dioxide The thermodynamic subroutines require that some carbon dioxide is present in the fluid in order to compute fluid properties Specifying pressure temperature and a nonzero carbon dioxide concentration fixes the in place dryness and enthalpy Alternatively specifying pressure nonzero carbon dioxide and an in place dryness fixes the temperature and in place enthalpy Specifying any three of the properties pressure temperature nonzero carbon dioxide concentration static dryness that is between one and zero or in place two phase enthalpy then sets required consistent values for the other two properties As a consequence the simulator internally sets a minimum mass fraction of XCO2 0 001 Note that this is a Mark McGuinness 2013 little less than the equilibrium value for water in contact with air at one atmosphere and at room temperature so is not an unreasonable minimum value Multiple feeds of four different types are allowed drawdown linear or nonlinear productivity index linear or nonlinear fixed flowrate or dry steam The user may specify multiple variations in casing geometry radius roughness inclination or a casing section may be chosen from an editable database file of casing types 2 1 System Requirements
70. dient The way forward is to use saturation as a measure of pseudo saturation This is uniquely given by the value of Ugg While Hagedorn and Brown acknowledge that pseudo saturation is not necessarily the same as actual saturation they also note that they take some trouble to use realistic values for the contribution of friction to pressure drop and this leads to pseudo saturation values that should be close to the true saturation The pseudo saturation S is seen in fig to be sometimes greater than the homogeneous 99 5500 5000 4500 4000 dP dz 3500 3000 2500 2000 10 20 P 40 150 Bara Q 50 kg s 30 40 50 60 70 80 USG 90 Figure 36 A plot of the pressure gradient obtained Pa m versus Usa value input to the Hagedorn and Brown correlation for each of five different pressure values Parameter values set are S 0 44 pressure range is 40 150 bars Q 50 kg s 5000 4500 4000 3500 3000 dP dz 2500 2000 1500 1000 500 40 50 60 70 gt PGRAD i FRIC i GRAV s S homog Figure 37 A plot of the pressure gradient obtained Pa m versus Usa value input to the Hagedorn and Brown correlation for P 40 bara S 0 8 Q 50 kg s Also shown are the contributions from friction and gravity as pressure gradients Also plotted are the pseudo saturation S and the homogeneous saturation which forms an upper limit on what the value of true satur
71. e it will automatically create a basic simulation from the built in template Then the setup can be modified or another setup can be imported from another folder to replace the current setup if there is another one available An important file in the SwelFlo project folder is the file called SWFinput txt It saves all of the latest details of the wellbore simulation and is the file that the simulator looks for when it is directed to a project folder More information on this main input file may be found later in this manual Fig shows the resulting windows that should be visible at this stage The main window contains tabs for STARTUP DATA INPUT RESULTS display and GRAPHS display The smaller Progress window is to provide information to the user The resizable Error Reporting region at the bottom of the main window is to report more serious errors found in the setup or encountered while running Project Help Pre Process Run Stop Edit Search Q snae ProjectName project Last Modified 09 09 2013 19 23 36 Downhole Profile requested Output curve not requested Topdown simulation completed Change or Create Project Directory parent A template used in absence of an input file Prodn flash in well roject Descipicn One feed topdown flow Orkiszewski correlation No outputfile asked for No reservoir T profile one casing section 1000m vertical Ready to run SETUP DATA INPUT RESULTS GRAPHS Flow Correlation Orkiszewski
72. e bottom of the INPUT SETUP tab Preprocess saves input values from the input windows to internal arrays does some checking of the setup to ensure it is ready to run and if it is ready to run then the input file SWFinput txt is saved to the current project directory as displayed in the Status Bar at the bottom of the main SwelFlo window Run will run the simulation it will also Preprocess if necessary first Note that hence it is always possible to ONLY press Run without first pressing Preprocess Bottomup runs and or output curves requested will cause a progress bar to appear the simulation progress bar show progress during a bottomup search for a match to whp and the output curve progress bar monitors progress during the generation of an output curve by running bottomup with various bottomhole pressures as starting points The stop menu and stop button will be active during simulation and cause an immediate stop in the middle of simulation if pressed returning simulator control to the user Note that there is no need for a Close menu item or button opening an input file is really just an import process which reads in the setup then closes the input file If it is desired to wipe clean all existing input values the user should select mouse drag or shift click and Clear values inside each grid window using the grid window menu Clear An alternative is to import the setup from another project or from an input file that the user has saved es
73. e digits or to save the csv file with any other name if desired This is to give some flexibility in saving wellhead setups for example if the user wishes to import read open them from another project from within the corresponding grid window These saved csv files are what may be opened by the Feeds window menu File gt Open Trying to Open an incorrectly shaped csv file usually results in the exe file crashing 4 9 Pre Process Run Stop Re Plot The first three of these buttons at the bottom of the INPUT SETUP Tab have the same effect as the corresponding menu items Pre Process saves input values from the input windows to internal arrays does some checking of the setup to ensure it is ready to run and if it is ready to run then the input file SWFinput txt is saved to the current project folder RUN will run the simulation it will also Preprocess if necessary first Note that hence it is always possible to ONLY press Run without first pressing Preprocess Bottomup runs and or output curves requested will cause a progress bar to appear the simulation progress bar show progress during a bottomup search for a match to whp and the output curve progress bar monitors progress during the generation of an output curve The STOP button and stop menu will be active during simulation and cause an immediate but non catastrophic stop in the middle of simulation if pressed returning simulator control to the user The Re Plot button allows the u
74. e in place It is related to the in place dryness or steam quality X by Xp t 1 X py S a 3 and to the flowing dryness by XfowPl X flow Pl K 1 Als Flowing dryness can be expressed in terms of static dryness X as Xu XUy 1 X u X fow and in terms of saturation S and flow ratio K as SpK X ow gt SpyK 1 S pi and flowing enthalpy can be written in the form R tow Ale 1 E X tow ht Other useful rearrangements of these equations are for static dryness X Spy X Spy 1 pe S5 pi and X flow AE and for flowing steam quality S Puty X fow G C 2 Mass Conservation of Water and Carbon Dioxide Note that mass flowrate input and output is always the total mass water plus carbon dioxide The relative mass of carbon dioxide is always very small though The equation for conservation of mass of water is simply in the absence of feedpoints Q Alpr 1 S u pySuy c 51 where c is some constant The same equation applies if M is the total mass of water and carbon dioxide in the vapor phase and M is the total mass of liquid water and dissolved carbon dioxide This is what is coded Further details on conservation of carbon dioxide and its consequences for computing flowing steam fraction may be found in appendix Note for now that Xgow is given by conservation of carbon dioxide once the constant total mass flowrate is known tog
75. e independent of depth An exception is that Taitel s horizontal flow modelling development 84 has variable area correctly included The work by Yadigaroglu et al 37 eqn 11 has the same form of momentum conservation equation as presented here in equation 111 and also presents the AG v form of the momentum flux eqn 14 The sign convention used above that the velocity terms G Um and v are positive for flow in the direction of increasing z that is downwards in the well is the opposite of the usual geothermal convention of positive velocities for upflow Then the geothermal convention if used would reverse the sign on the friction term only 120 J Productivity Index at Feedpoints The development in this section is explained more fully in The usual definition of feedpoint productivity P is in t hr bar Q where Q is the flowrate at the feed t hr Pes is reservoir pressure in bars and Pasen is wellbore pressure in bars This applies at a fixed downhole pressure and drawdown and hence at given wellbore and reservoir viscosities When increased production rates are simulated for an output curve the wellbore viscosity can reduce appreciably especially if the flashpoint reaches the feedpoint To account for this change in viscosity the above equation is usually rewritten guided by Darcy s law for flow as Q gt Poes Pye 115 where the number is called a productivity index units m3 v is the kinematic
76. e is no choice the header line The subdivision length DELZ should divide exactly into the section length and determines the distance between computational nodes If it does not divide exactly SwelFlo calculates a subdivision length that does with the nearest integer number of subdivisions It is unusual to go much less than 4m in subdivision length except internally at the flash point The total number of all of the subdivisions in the well cannot exceed 5000 in total Any feed points listed later in the input file should be exactly at the join between two subdivisions If they do not SwelFlo ensures that a subdivision is split into two at the feed depth The total well length should match exactly the sum of section lengths The deepest feedpoint should appear exactly at the bottom of the well but if it does not the simulator tries to place it there The maximum total number of subdivisions overall is 5000 There are a number of built in checks and fixes of subsection lengths that do not exactly divide into section lengths and fixes for overall LENGTH Checks have been added in the following order 1 if DELZ does not divide evenly within 1E 07 into SECL it is altered so that with the nearest integer number of sections it does 2 If bottom feed depth does not match LENGTH of well within 0 01 length is set to bottom feed depth 3 If sum of section lengths does not match LENGTH within 0 01 LENGTH is adjusted to match the sum of se
77. e plots are as chosen by the user for depth pressure and flowrate Otherwise temperature is in C and velocity is the mass weighted velocity expected to be obtained from a spinner in m s Note that the data files provided by the user are assumed to use the units chosen by the user for that well Simulated results are converted from the internally used SI units to the user units before displaying together with any data provided 8 Menu Items The menu in the main window has five major items Project Help Pre process Run and Stop Help Pre Process Run Stop New Project Folder Import Existing SwelFlo setup Load Template Switch to existing project Save ce of an inputfile Change Units rkiszewski correlz Exit SETUP DATA INPUT RESULTS GRAPHS casing section 1 Figure 15 The Project menu in the main SwelFlo window 32 8 1 Project This menu item has five submenus New Project Folder Import Existing SwelFlo setup Load Template Switch to existing project Save Change Units and Exit 8 1 1 New Project Folder This menu item has the same action as the Change or Create Project Directory button at the top of the main window They allow the selection of different existing parent directories to work within or the creation of new ones Selection of a different existing Project directory with a SWFinput file already present results in that file being read by the simulator and loaded to intern
78. e previous subsection dP Ge fmt 69 where r is wellbore radius and fy is a Fanning friction factor that depends on roughness and on a Reynolds number Reg which is given by at Reg ea Li where Pm Spy 1 Ni 5 pi Note that this frictional pressure gradient has been derived from the usual theoretical formula 47 D There is a slight difference in GWELL and SwelFlo convention so that fp 4fy 79 where fp is the Darcy friction factor due to White 1974 that has value 64 Re for laminar flow and fp is the friction factor used in SwelFlo F 4 2 Slug and Churn Flow In this regime Hadgu uses the results of Dukler et al 13 summarised here as dP Imp Ving Z EE aa 70 where Pre arspy 1 ars pi 0 3py 0 71 and fm depends on the Reynolds number _ DUzrspre Mi l To compute these we need Uzzs and 6 Referring to the previous subsection F 3 3 we compute given liquid and vapour phase velocities and saturation S Rem Usa Sty 71 Usz 1 Sju 72 DgAp Uso 0 345 73 Pl Un 1 29 Use Us Us 74 Urre 9 916 gD 1 arg 9 9498yr 75 1 ULLS Un Un F UrrB G azz 0 857Uy 0 148U rrp 76 QLS using arg 0 89 azs 0 3 We also have Xp a 77 su Kar l Aja 77 asu QLS asu 0 3 B 22 2 0 59 F 4 3 Annular or Mist Flow Hadgu notes that the flow of the liquid phase on the walls dominates frictional losses in this
79. e pure water saturation value On the other hand if the temperature of vapor phase is lowered it begins to form liquid droplets at almost the same temperature as the pure water saturation value Hence the single line that is a saturation curve for pure water opens out into a two phase region when carbon dioxide is present The more carbon dioxide the wider the two phase region opens out An illustration for a percentage by mass of carbon dioxide XCO2 0 01 is in Fig 40 The conservation equations for pure water still hold provided the liquid water phase is replaced with the liquid phase water and dissolved carbon dioxide and the steam phase is replaced with the vapor phase with partial pressures of the water steam component and of the carbon dioxide component adding to produce the total vapor pressure The presence of even just a little carbon dioxide then means that the two phase line is replaced by a two phase region so that in principle the static steam dryness X is now dependent on pressure and temperature and XCO2 Given any three values in the set X P T XCO the other value is then determined for a two phase fluid Usually XCO is given or known and P is retained as a primary variable SwelFlo retains T as the third primary variable and determines X and other secondary variables like enthalpy as a consequence H 1 Phase states A sequence of checks provides the state of the fluid liquid water plus carbon dioxide in s
80. e radius of 0 1m a pressure of 70 bara and a static steam fraction of 0 5 The regions for Q 5 10 kg s where the correlation becomes multivalued are at the upper end of the slug region and the lower end of the transition region When Q 15kg s and higher the correlation is single valued and invertible and covers the entire S range of 0 1 When Q is less than 5kg s the correlation is single valued but does not reach to S 1 1 0 9 0 8 0 7 0 6 0 2 4 6 8 10 12 14 USG Figure 31 An illustration of issues with the Duns and Ros correlation at low flowrates 1 black 5 red 10 green and 15 kg s blue a wellbore radius of 0 1m a pressure of 70 bara and a static steam fraction of 0 5 Usa values are input from zero to maximum for the given flowrate and S is computed according to the Duns and Ros correlation as discussed in the text We observe that higher saturation values are not realizable at low flowrates and that sometimes several Usg values correspond to one saturation value making the correlation noninvertible there F 5 7 Duns and Ros Momentum Equation If the full Duns and Ros correlation is chosen the treatment of momentum terms is different to the usual SwelFlo treatment and is as detailed here F 5 8 Friction Pressure Changes The usual SwelFlo treatment is to use a two phase multiplier based on Martinelli and Nelson and modified by Chisholm 10 Duns and Ros take a slightly different approach dependi
81. egion or the vector T P in a single phase region then the method is 23 Newton Method with Backtracking and the Armijo Rule Set To J L o f x Momtol Ento Entol and Momtol are used to normalize or weight the momentum and energy differences so they have equal influence on the size of f Hence the critical value it is desired to reduce f to is V2 While f x gt 7 70 Ta where 7 is a relative error and 7 is an absolute error we will choose Ta V2 and 7 0 we do the following to take one modified Newton step where f is the gradient of f If f x lt 107 exit and try reducing subsection size Choose search direction d 4 the Newton step direction Fe Set A 1 Set x x Ad a trial point in the Newton direction If f x4 lt 1 aA f x then set x x accepting this step else set A A 2 and go to j again Once the step is accepted this entire process is repeated from Typically a is set to 1074 The process will stop when f is small enough in practice a limit of 50 is placed on the total 59 number of iterations allowed The process is linear until it steps close enough to the root then it is quadratic as Newton s method is E Conservation equations at Feedzones The sign convention is for positive flow values for production hence Q is positive for flow up a wellbore and for flow into a well from the reservoir at a feed At a
82. elations is that there are no flow regimes and no true value is computed for saturation S The emphasis of Hagedorn and Brown is on matching observed pressure changes for a water air mixture where presumably energy conservation is not required to be considered This makes the Hagedorn and Brown correlation of limited interest in the geothermal context where temperature changes are also important Futhermore implementation is challenging for a geothermal wellbore simulator as the primary output of the correlation is the pressure gradient in a situation where temperature is irrelevant The output that is needed is saturation given flowing steam quality that is given U SG and U sL It is not possible to use pressure gradient to obtain saturation since as noted below the relationship is not invertible one pressure gradient corresponds to two different saturations The way forward is to use the pseudo saturation that the Hagedorn and Brown correlation gives as a measure of steam saturation In fact since the pseudo saturation is sometimes larger than the physically realisable maximum saturation obtained for homogeneous flow the minimum of pseudo saturation and homogeneous saturation is used as actual saturation S Some of the following discussion relates to the computation of pressure gradient which is how the Hagedorn and Brown correlation was originally set up However the saturation we extract from the correlation does not requi
83. elocity profiles as depicted in Fig 3 To change the wellhead conditions click the SETUP tab and then click the Wellhead button The Wellhead Conditions grid window will open with two rows one row for each desired Project Help Pre Process Run Stop Project Name project i LastModified 09 09 2013 22 18 24 Change or Create Project Directory Poeci A template used in absence of an inputfile Prodn flash in well rojec SS PERET A ATAR Desciiption One feed topdown flow Orkiszewski correlation No outputfile asked for No reservoir T profile one casing section 1000m vertical Ready to run SETUP DATA INPUT RESULTS GRAPHS P Tee Mama _OC P Bara 25 30 10 15 20 35 40 45 Simulation 100 200 300 400 500 z mMU 600 700 800 900 1000 Topdown simulation completed a E Current Directory C Users Markm Documents SwelFlo GenericWell project Figure 2 SwelFlo first steps after running showing simulated pressure versus depth profile PaT ly Poel V 170 175 180 185 EEES 200 205 210 215 A x P A 4 v m s 6 3 o 100 igo 200 200 300 300 _ 400 400 E 500 E 500 600 600 700 700 Simulation Simlation 800 800 900 900 1000 1000 Figure 3 SwelFlo first steps showing simulated temperature and spinner velocity profiles
84. eothermal brine can be treated as a solution of NaCl in water to evaluate its fluid properties The principal effects of dissolved solids are boiling point elevation increased viscosity increased density increased surface tension and decreased specific heat 29 I Corrected Area Treatment A careful derivation of the equations for conservation of momentum reveals that the usual form taken for the acceleration term that is the momentum flux term is only correct for constant area and is incorrect at places where the casing size changes Correct versions can be found in 251 841 37 and are reviewed briefly here The momentum conservation equation after careful averaging of the three dimensional form across the wellbore cross sectional area is usually written as dy 10 Oz Ade where lt pu v gt is an appropriate average value pm Sp 1 S p is the usual two phase average fluid density and f takes account of two phase correlations but is related to the Fanning friction factor which is equal to the Darcy or Moody friction factor divided by four The correlation coefficients Cy and Cy are less than but of the order of one and account for the fact that the average of the square of velocity is not the square of the average velocity over the cross section of the well 2 AC ow S pou ACy 1 S piu El lt pv v gt pmgsin 111 The first term on the right hand side the momentum flux term is often rewritten
85. equations but starts with the given wellhead setup in particular with a known flowrate and wellhead pressure and predicts the resulting flow and pressure at bottomhole Injection of single phase fluid is handled correctly For two phase production the rise velocity of vapour bubbles imposes a minimum upwards mass flow velocity in the correlations typically about 1 kg s 4 2 1 A single bottom up simulation If a negative reservoir pressure is provided at the bottom of the well for a bottom up simulation this is taken to be a signal that the absolute value of this pressure is to be used as 17 a single bottomhole flowing pressure wellbore pressure for a single bottom up simulation There is no attempt to match whp in this case Feedpoint flowrate must be specified if this option is used An alternative way to ensure this is to enter the same value of bottom hole maximum pressure and minimum pressure in the Input Setup tab then reservoir pressure can also be separately specified and the bottom feed can be of any type This can be a useful option for comparing bottomup and topdown simulations a topdown simulation can give bottomhole flowrate enthalpy and wellbore pressure Setting these and running a bottomup simulation should give the same downhole profile The advantage of this over the alternative option of setting a negative reservoir pressure at bottom hole is that here a different positive reservoir pressure may be set so that
86. ero or one If fluid is not moving the dryness X determined by the above equation is the static steam fraction since the mass of carbon dioxide in Mwater kg of water gas mixture consisting of Miiguia Of liquid phase and Myapor of vapour phase is given by XCO M water XC Miiquia P XC Myvapor which on dividing through by Mater gives XCO XC 1 X XC X In the above equation X Myapor Mwater is the static steam fraction It rearranges to give eqn 109 If fluid is moving as in a producing well conservation of carbon dioxide gives XCO Q XC M XC M which on dividing through by Q gives XCO XC 1 Xg XC Xoo 109 since Xq M Q This rearranges to give eqn 109 with flowing steam quality that is _ XCO XC R i XC XC 110 Hence flowing steam quality is determined by mass conservation P T and gas fraction XCO when the two phase mixture is flowing in a wellbore H 5 Finding Pressure The equation 110 may be rearranged together with the law XC Y 52 to obtain an equation that gives a quadratic equation for pressure in terms of Xao T and gas fraction XCO as follows P P P Ps XCO Xa P Xe aap where P is the saturation pressure for pure water in bars absolute at temperature T This rearranges to give the quadratic in P in bars not Pa aP bP c 0 where DOCO 1X b A XOOz Xa BPO XCO PUA C P Xgol A E BP where A B an
87. erwise if 8 gt Lg and vp lt Ls the regime is slug To get a smooth and monotonic function S of 8 there needs to be a transition region between this and the bubble regime called bubble slug It is not enough to simply set the vapour velocity to transit linearly from bubble to slug values say as 8 goes from 0 9Lpg to 1 1Lpg since although this gives a continuous transition in vapour relative velocity it does not always give a monotonic S 68 This can lead to discontinuous jumps in saturation S and consequential failure to converge when searching for values of T and X that conserve mass momentum and energy The resolution of this issue in this regime is to find the value P for which the implicit relationship 6 Lg holds this relationship is illustrated as the intersection of y P and y Lg in Fig 21 and to use this value 6 to compute Xfow Ur and Uy Ur Ug and use these to calculate S P As parameters like flowrate are varied the piecewise defined curve for Lg changes keeping the level part at 0 13 but the parabolic part moves back and forth Hence solutions to 6 Lg vary between 0 13 and 1 0 Then we set a transition region with a slightly positive slope SL using the equation Sbu sl S 6 ES SL P B i Sslug in the slug regime is also determined by u vr Us and Sstug oe The extent of the transition region is determined by the inequality Sbu sl gt Sue G That is while this inequality holds
88. es are RESV 10 D C and RESV 11 D n both floating point where C and n are the power law parameters in the dry steam power law Q C P RN h cel and usually n 0 75 and C ranges from 0 4 for typical producers to 0 75 for good producers to 1 5 for the very best producers Note that in this formula however pressure in the reservoir and in the well are first converted to bars before squaring etc and flowrate is taken to be in t hr so is then converted to kg s afterwards So C and n must be chosen accordingly for pressure in bars and flowrate in t hr Dryness may be specified not equal to one but this is not very appropriate for a dry steam feed 47 If FICODE is 4 linear or quadratic drawdown the whole line should contain F_DEPTH I FICODE I RESV 1 D RESV 3 1 RESV 12 D FDRY D RESV 6 RESV 8 D RESV 9 1 The new variables are RESV 8 I c the floating point linear drawdown in bars t hr which is in fact one over the usual productivity in t hr bar and RESV 9 I d the nonlinear drawdown term for flowrate in t hr and pressures in bara in the drawdown equation o Pwel cQ ag dQ Q To get just linear drawdown you can set d 0 in the input file Note that here as for dry steam pressures are first converted to bara and Q is assumed to be t hr so is then converted to kg s afterwards THE END The input file is read until end of file looking for all feeds provided Nothing further is sought in the
89. es of CO Proz The results are shown in Figure 5 3 28 TABLE 5 3 VALUES OF COEFFICIENTS FOR CALCULATION OF CO VISCOSITY 2 96610E 03 2 85290E 06 2 182905 09 2 58250E 01 7 11780E 04 6 95780E 07 9 00870E 01 2 47270E 03 2 41560E 06 1 12474 2 98864E 03 2 85911E 06 8 50257E 01 1 99076E 03 1 73423E 06 6 73731E 01 1 419905 03 1 13548E 06 5 00750E 01 9 04721E 04 6 19087E 07 4 08927E 01 6 35032E 04 3 53981E 07 Although the results showed that interfacial tension decreases with increasing Poop partial pressures greater than 10 bars at the wellhead rarely occur in geothermal well discharge fluids At Poo lower than 10 bars the decrease in surface tension is less than 15 Therefore in this study the interfacial tension of H O CO is assumed to be approximately the same as that for pure water On Ouo 5 27 where Sm surface tension of mixture Ta surface tension of water 5 2 WATER SODIUM CHLORIDE SYSTEM H O NACL The total dissolved solids in geothermal brines varies from that of ordinary well water up to concentrated solutions as high as 40 by weight Sodium chloride NaCl is typically 70 to 80 of the total dissolved solids The other most abundant components are potassium chloride KCI calcium chloride CaCL and silica SiO The silica concentration in geothermal brines is usually between 200 and 600 ppm Wahl 1977 Since NaCl is the major component of the total dissolved solids the g
90. es required There are a number of typos in GWELL Table 4 3 as well as in Chisholm Table 4 4 A corrected version of these tables is presented in Table 4 for clarity D is wellbore diameter and negative signs are chosen for downflow positive signs for upflow The terms v and v represent specific volumes v 1 p and v 1 p The drift velocity is 1 Xpow uwo ty un E o Jan A Pl GWELL mistakenly uses the slip velocity in place of the drift velocity in the Fortran code and replaces the specific velocity Ug in the last line of the table with the liquid viscosity u in the Fortran code and in the GWELL manual Note that if 8 lt 0 9 we need the drift velocity to apply the correlation and we need K to find drift velocity but K depends on saturation S That is C4 depends on both S and 2 So the formula S BCA B S defines saturation implicitly in terms of 8 In this range of 6 we want to find saturation but we need saturation to apply the correlation If 8 gt 0 9 we need the specific velocity Us 1 Xgow Q Ap to compute C4 This can be done without knowing saturation GWELL deals with the implicit nature of the Armand correlation by first estimating saturation by using C4 Can then computing S Cap Then GWELL calculates a value for K K using this value of saturation This value for K is then used to compute C4 from the table Chisholm deals with the implicit nature of the correlation by using Labu
91. ether with pressure temperature and mass fraction of carbon dioxide The following conservation laws are presented for pure water only C 3 Momentum conservation The momentum conservation equation is 5 lt u gt ae lt p gt gsind dis dP dz dz do dz where lt u gt 1 Xgow U Xgowuy is an average flowing velocity lt p gt 1 a p apy is the usual average density z is distance down the well and G Q A is taken to be positive for upflow production The first term on the RHS is the divergence of momentum flux in the case of constant A when G can be put inside the derivative the second is due to gravitational force and the third is due to friction forces between phases and wall parameterised as a wall loss term The average density lt p gt can be written in terms of flowing steam fraction by using equation 4 to obtain X flow K 1 X flow X flowPl K 1 A pv l The original code in GWELL is missing the term K in the denominator of the previous equation so it does not conserve momentum lt P gt Pipo Note that the momentum flux term is not always written correctly for varying wellbore cross section area A 25 The term G is often put inside the z derivative which is only correct for constant A Uy X flow l a pi 7 Ul 1 X flow Q Pv it is straightforward to rewrite the momentum flux term in the case of constant A in the alternative forms Note that using SS
92. ex approach is not appropriate Hence in both topdown and bottomup cases a simulation that attempts to match a given well head pressure is used except for the dry steam power law case to produce a productivity index at feeds where required especially at bottomhole and these indices are used to generate an output curve Generation of an output curve requires that a reservoir pressure be provided at least at bottomhole If a fixed flowrate is specified at a feedpoint and a positive reservoir pressure is provided there the output curve will use the internally computed productivity index to calculate flowrates there not the fixed flowrate The matched bottom up simulation however will use the fixed flowrate This can lead to differences between output curve and simulation run if the reservoir pressure is inconsistent with the sign of the fixed flowrate specified for example if reservoir pressure is below wellbore pressure and a positive flowrate is specified If the user wishes to retain the fixed flowrate at a feedpoint during computation of the output curve they can impose this by setting the reservoir pressure to a negative value at that feed 4 11 Rock Properties Heat flow between casing and country may be turned on and off by selecting Yes or No in the Allow Heat Flow To From Country radio buttons in the Rock Properties section of the SETUP tab in the main SwelFlo window Values for thermal conductivity density hea
93. expression is a good fit to the experimental data of Takenouchi amp Kennedy 33 Furthermore Sutton notes that the above expression relating the masses of carbon dioxide and water in the vapor phase rather than the mole fractions Daltons Law says the mole fraction should be used gives a better fit to the data of Malanin 24 Todheide amp Franck 35 and Takenouchi amp Kennedy 33 for pressures between 100 and 300 bars The fit is within 5 for temperatures in the range 110 280 C and within 25 for temperatures in the range 300 330 C An alternative approach for ideal gases is based on Dalton s law Dalton s law states that the sum of the partial pressures of carbon dioxide and of water is the total pressure Combining Dalton s law with the ideal gas law gives P P Pecos nu RT V n RT V where P is the partial pressure of water component or saturation pressure P is total pressure and n is the number of moles of water present in the vapor phase and ne is the number of moles of carbon dioxide present in the vapor phase both in the same volume V of vapor and at the same temperature T Then it follows that Poo2 Me P netnw where YC n N Nw is the mole fraction of carbon dioxide in the vapor phase Since the molar mass of carbon dioxide is 44 and of water is 18 the mass fraction of CO in the vapor phase can be written as YC 44 YC XC Co 44 YC 18 1 YC so that mole fractio
94. fective viscosity and productivity index and for a discussion on their relationship to the more usual productivity in t hr bar Note that the usual productivity in kg s Pa is the productivity index divided by effective viscosity The productivity index and perhaps also A for quadratic drawdown is needed for correct output curve calculations since the effective viscosity changes with pressure and especially also if flashing occurs at a feedpoint at higher flowrates and this changes the flowrate Note that if the feed is specified to be a fixed flowrate feed or a drawdown feed and if an output curve is requested and a positive reservoir pressure is provided at the feed depth then a productivity index is computed from the topdown or bottomup run that produces a downhole profile This productivity index is then used in computing the output curve and the feed is temporarily treated as a productivity index feed for the purpose of computing the output curve Fixed Flowrate Feeds The only extra parameter required here is the fixed flowrate desired at the feedpoint set in column eight Positive flow is for production that is flow from the reservoir towards the wellbore Reservoir pressure is only required for this feed type if the user wants a productivity index calculated Dry Steam Feeds The extra parameters to be provided for dry steam feedpoints are C and n in columns eight and nine the power law parameters in the dry steam power law m
95. fication is that occasionally the Newton search for a solution to the conservation equations gets stuck in a two cycle stepping over the region where the solution is especially for the Hadgu correlation at the change from churn flow to transition where there may be large changes in slope of the curve of vapour phase velocity versus saturation This behaviour is analogous to the one dimensional behaviour of Newton s method when trying to find the root of tanh x 0 if one starts with too large an x value it steps further and further away from the solution x 0 The idea is to use Newton s method to find the direction to step say in T Xgow space for two phase flow but then to backtrack check the size f of the error vector with components being the change in energy and the change in momentum after stepping and ensure this size reduces from the current value before accepting the Newton step If the size does not reduce the length of the step is halved but the Newton direction is retained we backtrack and the size f checked again This halving and checking is repeated until a smaller size is achieved This works for the tanh function and fixes the issue with rapid changes in the slope of the energy difference 23 pp 169 ff If Fg and Fy are the errors in energy conservation and momentum conservation respectively and it is desired to have these less than Entol and Momtol respectively and if x is the vector T Xgow in a two phase r
96. file crashing 4 8 Feed Points The FEEDPOINTS window allows the user to set or modify the measured depth of up to ten feedpoints the type see below the reservoir pressure the reservoir temperature or enthalpy XCO2 in the reservoir and optionally the dryness of the reservoir fluid near the feedpoint Other parameters are also set depending on the feed type or code Dryness steam mass fraction is optional and if entered and nonnegative it will cause the enthalpy or temperature to be overridden to be consistent with the desired dryness at the given reservoir pressure The 22 resulting enthalpy is then taken to be the flowing enthalpy of any reservoir fluid that enters the feed fe FEEDPOINTS File Edit Search Save SaveAs EditRow Clear DERA ana Dapth mMU Type Res P Bara Res T oC Res H kJ kg ResX ResXCO2 sigma 1 Q 2 kg s C 3 4 A 1 n 3 d 4 1000 0 prod index 55 000000 212556 o 9100000 03 0 0000000E 0c 0 0010 0 30000E 11 0 00000 Pi Row 1 Column 1 Figure 11 A view of the Feedpoints grid window Any value of enthalpy entered for a feed point is interpreted as flowing enthalpy since it is flowing enthalpy that is required for applying energy conservation laws If a temperature is set for the reservoir at a feed point and no enthalpy or dryness values are set then this temperature is used to determine a static enthalpy and this is then taken to be the flowing e
97. files are assumed to be the same as the units the 29 SETUP DATA INPUT RESULTS GRAPHS Add Data File Type Select loaded file to modify Edit Remove Re Plot D P __ C SwelFlo testruns p 1 p2 PDataSim003 b i ischarge Test C SwelHlo testruns p 1 p2 T DataSim003 bt Mass Flow kg s WHP Bara Enthalpy kJ kg on first line of data file Existing New _ Pressure mMU Bara New P file Temperature mmu oC New Tile Fluid Velocity mMU m s Output Test kg s Bara New OC file Figure 13 Screen snap of the DATA INPUT tab in the main SwelFlo window user has most recently chosen in the menu item Change Units and they default to Pa and m For the output curve the default flowrate units are kg s Note that if the user changes units in a project the simulator makes no attempt to change the units in the data files in the project directory The data files are always assumed to contain data that is in the units currently chosen by the user for the current well Hence the user needs to ensure that data units match their chosen set of units by directly editing data files from inside the simulator or by using an editor like Notepad from outside the simulator Enthalpy data at the wellhead is simply passed through to csv files if provided in P T or V data files so the units are up to the user Enthalpy is simply passed through to csv files not plotted or used Temperature is always in degrees C
98. flow of gas and liquid mixtures in wells Proc 6th World Pet Congress Frankfurt Section II Paper 22 PD6 pp 451 465 1963 June 19 26 Ellis A J amp Golding R M The Solubility of Carbon Dioxide Above 100 C in Water and in Sodium Chloride Solutions American Journal of Science Vol 261 1963 pp 47 60 39 16 17 18 19 20 21 27 28 29 30 31 Fernandes R C Semiat R amp Dukler A E Hydrodynamic Model for Gas Liquid Slug Flow in Vertical Tubes AIChE Journal Vol 29 No 6 1983 pp 981 989 Garg S K Pritchett J W amp Alexander J H A new liquid hold up correlation for geothermal wells Geothermics 33 2004 pp 795 817 Grant M A Donaldson I G and Bixley P F Geothermal Reservoir Engineering Academic Press 1982 Hadgu T Vertical Two Phase Flow Studies and the Modelling of Flow in Geothermal Wells PhD thesis University of Auckland 1989 Hadgu T amp Freeson D H A Multipurpose Wellbore Simulator GRC Transactions Vol 14 Part II pp 1279 1286 Aug 1990 Hagedorn A R amp Brown K E 1965 Experimental Study of Pressure Gradients Occurring During Continuous Two Phase Flow in Small Diameter Vertical Conduit J Pet Tech Volume 17 Number 4 pp 475 484 Kjaran S P and Eliasson J 1983 Geothermal reservoir engineering lecture notes University of Iceland UNU Geothermal Training Programme Iceland Report 983 2 http www
99. he Hadgu correlation with transition regions included when G 1718 T 246 C The solid circles delineate the bubble slug around S 0 3 and transition near S 0 6 regions Note that as G decreases to zero C2 oo and C3 becomes complex valued as C2 passes through the value one As this happens S increases to equal the value one and is not real for values of Cy greater than one This signals the situation mentioned earlier that there is then no mist regime For a typical wellbore this occurs as flowrates are reduced through values of about 5 kg s As for other correlations things begin to fail at small flowrates where counterflow becomes possible Counterflow usually vapor flowing upwards and liquid flowing downwards is not simulated at this stage in SwelF lo As G increases Co vanishes C3 increases and S tends to zero When S reduces below 0 32 this is taken to mean there is no slug or transition region between bubble and mist regions only the bubble slug now bubble mist transition F 3 8 Implementation logic Forward Problem If S lt 0 28 bubble flow is used to find Us Otherwise Compute S where mist flow Ugg 1 0 If S lt 0 32 bubble slug transition flow is used to find Us If S lt 0 32 this will be bubble mist transition Otherwise Compute S where slug flow Ug is a little greater than mist flow at S If S lt S slug flow is used to find Us Note that this will not arise if S lt 0 3
100. he three lines of Project Description text may be used as desired to describe the current project They appear in the output files generated for each simulation Error information is offered in the bottom section of the main window while progress reports appear in the small resizable window in the top right part of the screen To illustrate some of the features of the simulator the next section shows how to compute a bottomup simulation that matches a topdown simulation and how to compute an output curve 3 3 Bottomup Match To start from a common baseline click Project Load Template and agree to the subsequent over writing of the current project Click RUN button or menu We will now set up a plot of the results just obtained from this Topdown run for later comparison with a bottomup simulation We will treat them as data Click DATA INPUT existing Pressure Choose the latest PDataSim file offered Similarly choose T and V files generated by the most recent simulation performed in the previous paragraph Re Plot top right in current window or bottom of SETUP window Note the plots now also have symbols from the most recent simulation spaced apart every 50m This spacing can be altered using the box in the bottom right corner of the SETUP window then RUN and reload the data files and rep lot as above The symbols will trivially lie exactly on the lines since they came from the same simulation Now set up the feed at bott
101. he way Then the simulation arrives at the wellhead with values for flowrate enthalpy mass fraction of carbon dioxide and pressure that are determined by the computation Many such simulations are performed in bottomup mode and the value of bottomhole wellbore pressure is varied in an attempt to match the wellhead pressure To perform a bottomup simulation of a downhole profile and to compute an output curve the user must provide reservoir pressure feed properties and either temperature enthalpy or vapor mass fraction at bottomhole and at any other feeds The range of bottomhole pressures that are used to perform this search for a match to wellhead values defaults to a maximum of the reservoir pressure at bottomhole and a minimum of 1 bara The user may help the simulator by refining the pressure range to be searched on at bottomhole by entering two different pressure values into the boxes labelled Range of bottomhole P s to use at the lower end of the Output Curve section of the main SwelFlow window This range serves two purposes one to help the downhole profile to match wellhead values and the other to set the range used to compute an output curve The two modes have different uses Bottom up is useful for predicting the flowrate for a target wellhead pressure given feed and reservoir properties Output curves wellhead flowrate versus wellhead pressure require a series of bottom up simulations Topdown mode solves the same
102. his is ignored A different approach is used by Brill and Beggs 7 who use the Duns and Ros mist flow correlation to compute the acceleration term rather than the Hagedorn and Brown correlation This will make the pseudo holdup value incorrect when the acceleration term is significant since it was set by Hagedorn and Brown by fitting to data and using a different acceleration calculation Since Hagedorn and Brown essentially fit a pseudo holdup to measured pressure changes it is important to use the same method as the one they have used to compute the derivative of V2 irrespective of what actually happens Any deviation from the formulae used by Hagedorn and Brown will invalidate the fitted correlation for pseudo holdup A close examination of their original paper gives a hint of what has been done to compute the difference of V from one point to another point that is a pressure and temperature increment away Hagedorn and Brown comment 21 p 481 item 4 that they determine p at an average pressure and temperature and calculate the velocity of the mixture here presumed to mean Vm at both ends of the increment using the pressures at those two points and the ratio of the flow rates as measured at the out flow end of the tube presumably the ratio referred to is of Usa and Usz This suggests that the ratio Ug Usa is taken to be fixed at a value say R over the increment so that since pus pUse G and then us
103. in one metre of casing length with units J s m The first two terms are a flowing average enthalpy and the next two terms are an average kinetic energy of the two phase flow 54 C 5 Heat Flow between country and casing Following 5 a transient radial heat flow model is used for conductive heat flow between casing and country at depth z oT D oT t r r Z where D E is the thermal diffusivity of the surrounding rock Boundary conditions are assumed to be that the wellbore temperature is fixed at T while the temperature far from the well is the reservoir temperature Tes A transient solution valid for large times t gt r D gives _ Akt Tres Tw Qi A In ES 2n a where 7 0 577216 is Euler s constant This solution assumes diffusive heat transport between a constant wellbore temperature Tu and surrounding reservoir rock that starts out at the temperature Tes at time zero It does not take into account any changes in wellbore temperature over time D Discretisation and numerical method The wellbore is specified by the user as a number of straight sections with constant geometry joined end to end at various angles to the horizontal A section is subdivided into smaller subsections for numerical accuracy as illustrated in Fig 16 The user provides the length and properties of each section The first subsection is assumed to start at z 0 and the node at the top is labelled n 0
104. ing the ratio we have G Usa _ n Pv piRy At node 1 specific gas velocity is given by G AS py e Ru and at node 2 since G is the constant mass flowrate per unit area and R is assumed constant G es ey dd pe a Ry These can be rearranged to get 1 1 v P Ry ye yp 2 l 106 SG SG p di p R 95 which tends to Torrens formula in the limit that p R is much less than py Since this is the original formulation in 21 it provides the best way to calculate dV dP The change in V is then given by Ap RuApi AV Va Pu F Rap 1 maA UsLAp 107 where A means the value in node 2 minus the value in node 1 Since pressure gradient is not known a priori in the forwards correlation the implicit nature of the accelerational pressure derivative is accommodated by using dp __ y Wim _ _ y Wm aP a p y ave dP dP dz me dP az AP Vaz az where AV is as given in the previous equation so that computation of e in SwelFlo is given by AV PsVin dp dp Vig map D Uso Us gt 0 The acceleration pressure gradient is taken across to the left side of eqn 103 and the equation is made explicit in total pressure gradient to get 2 dt a gas dz l e 108 Note that a comparison between this and Torrens approach is possible ra one considers the physically reasonable limit that dor is relatively small so that equation
105. input file 48 B Nomenclature A c Ca D D fp fu g G h k K M Nig Nov Na N P Q Qi r Re S T u Usa Ust v lt v gt X X fow Subscripts wellbore cross section area m thermal capacity of rock J kg C Armand coefficient diameter of wellbore m thermal diffusivity of rock m s Darcy friction factor Fanning friction factor fp 4 gravitational acceleration me mass flux kg s m specific enthalpy J kg thermal conductivity J m s C ratio of phase velocities mass flow rate of a phase kg s dimensionless liquid velocity dimensionless gas velocity dimensionless casing diam dimensionless liquid viscosity pressure Pa total mass flow rate kg s heat flow rock to casing J m radius of casing m Reynolds number vapor saturation temperature O velocity m s superficial gas velocity m s superficial liquid velocity m s specific volume m kg average flowing velocity m s steam dryness static mass fr steam quality flowing mass fr homogeneous liquid phase reservoir vapor phase well 49 vapor saturation vol gas flow ratio roughness m density kg m ave two phase density kg m ave two phase density kg m surface tension water N m dynamic viscosity Pa s C Governing Equations The steady state equations that SwelFlo solves in sections of the well between feedpoints are outlined here C 1 Two Phase Flow Terminology The notation introduced here follows 10 and is also broad
106. is described in more detail in the next subsection Speed of computation is generally excellent due in part to initial estimated values of pressure and temperature for the Newton iteration being chosen as linear extrapolations of the previous step so that often only one or two Newton steps are required The exception to this is near phase changes It was found that stepping across a phase change can lead to larger errors than tolerances allow unless the code automatically reduces subsection size there So the code reduces step size down to subsections of length just 2cm at a flash point The momenta and energies on either side of the flash point are required to match across the flash point with single phase values on one side and two phase values on the other For the SwelFlo code tests show that flash depths are reproducible and are a match when the same simulation is run topdown and bottom up which provides a good test of accuracy and reliability A good match is also seen between pressures temperatures and velocities when comparing topdown and bottomup simulations Accurate flash point determination is important for pressure matching and 2cm gives an accuracy of better than about 200 Pa D 2 Modified Newton with Armijo Rule This is a summary of the modified Newton method that has been developed to solve for T and Xfow that satisfy the conservation equations in the next node for two phase flow in the 58 well The reason for the modi
107. kenouchi and Kennedy and by Ellis and Golding is used Pritchett et al fitted the following relationship to the data gathered by 83 XC Pcoz A BPco2 where XC is the mass fraction of carbon dioxide dissolved in liquid water and Pco2 P Ps T is the partial pressure of carbon dioxide in the coexisting gas phase The values of A and B depend on temperature as A a aT aT B bi ht T bT where the values of a and b are given in Table Qi bi 1 03549E 03 2 04465E 01 1 60369E 01 1 07449E 01 4 83594E 02 1 44701E 04 gt w N e Table 7 The coefficients in the quadratic dependence of A and B on temperature for calculation of the solubility of carbon dioxide in water H 3 Mass Fraction of Carbon Dioxide in Gas The mass fraction of carbon dioxide in the vapor phase which is in equilibrium with the liquid phase can be calculated in different ways One way is described by the following fit to the 106 results of Takenouchi and Kennedy 33 Poo XC gt P where XC is the mass fraction of carbon dioxide in the vapor phase This is used in GWELL with a statement that it is better than Dalton s law for matching downhole data A reference is also made to Pritchett et al 30 and in this manuscript a very detailed set of instructions for determining the state of geothermal fluid in the presence of carbon dioxide is presented In particular Pritchett et al note that for states of geothermal interest the above
108. length 0 01 1000 m casing deviation angle 0 90 S vertical extent of casing section 1 5000 m feed depth 0 5000 m vertical depth to reservoir T 1 5000 m reservoir T for profile 1 374 C reservoir pressure at feed 230 0 1000 0 Bara reservoir H at feed 4 3 3000 kJ kg A n or d at feed 0 0 3x 10 XCO at feed or wellhead 0 001 1 0 X static dryness at feed or wellhead 2 0 1 0 wellhead pressure 1 0 230 Bara wellhead flowrate 500 0 500 0 kg s wellhead or feed H 4 3 3000 kJ kg bottomhole pressure range 0 0 1000 Bara thermal conductivity 0 10 W K m rock density 0 4000 kg m heat capacity 0 4000 J kg K time since discharge began for heat transport to country calcs 0 10 s spacing for data plot output from simulation 5 1000 m relative tolerance for matching WHP 1076 0 5 max and min wellhead flowrates for o curve plot 0 0 2000 0 kg s Table 1 Ranges imposed by the simulator on input values of parameters Input values outside these ranges will not be accepted by the input windows The user can avoid entering a valid value once started by clearing the value using the window menu 15 SETUP DATA INPUT RESULTS GRAPHS Downhole Profile Flow Correlation Orkiszewski Yes Produce a profile O No Output Curve Simulation Direction E ii Produce an Yes Topdown Bottomup Output Curve No No of points in Wethead No panisi 50 Range of flowrates 0 0 kg s Range of 0 000000E 00 bottomhole P s to use
109. lifornia Earth Sciences Division Report LBL 31428 UC 251 103 pp Barnea D A unified model for predicting flow pattern transitions for the whole range of pipe inclinations Int J Multiphase Flow Vol 13 No 1 1987 pp 1 12 Brill J P and Beggs H D 1991 Two Phase Flow in Pipes University of Tulsa sixth edition Carslaw H S and Jaeger J C 1959 Conduction of Heat in Solids Oxford University Press 2nd edition Colebrook C F and White C M 1937 Experiments with Fluid Friction in Roughened Pipes Proceedings of the Royal Society of London Series A Mathematical and Physical Sciences 161 906 367 381 Bibcode 1937RSPSA 161 367C doi 10 1098 rspa 1937 0150 Chisholm D Two phase flow in pipelines and heat exchangers Pitman Press Ltd Bath 1983 DeGance A E and Atherton R W 1970 Chemical Engineering aspects of two phase flow Part 6 Vertical and inclined flow correlations Chem Eng 77 pp 87 94 October 1970 Drew D A amp Wood R T Overview and taxonomy of models and methods for workshop on two phase flow fundamentals in Proc 1st Int Worshop on Two Phase Flow Fundamentals National Bureau of Standards Gaithersburg Maryland Sep 22 27 1985 pp 1 88 Dukler A E Moye Wicks III amp Cleveland R G Frictional pressure drop in two phase flow B An approach through similarity analysis AIChE Journal Vol 10 No 1 1964 pp 44 51 Duns H Jr amp Ros N C J Vertical
110. lots unless it is Run again New output files do not overwrite old ones due to the numbering system used for the filenames see below A number of properties of the flow may be set at wellhead the pressure and XCO are considered fundamental while temperature T enthalpy H and steam mass fraction X are considered interdependent Note that a distinction is also drawn between static and flowing values of enthalpy and steam mass fraction which coincide for a single phase flow or for homogeneous two phase flow Any one of these five variables is sufficient to determine the other four given P and XCO 3 1 1 Consistent Wellhead Values There is a hierarchy in the case that more than one of T H Hs X X 510 are set by the user at the wellhead to ensure that the values are consistent The hierarchy is that X fio is most dominant a value between 0 and 1 will override any value that is set for T H X or H pio Next is H sio which will override any value that is set for T H or X Next is X which will override any value that is set for T or H Next is H which will override any value that is set for T Hence the T value that is set will only have effect if none of H Hyio X X fio are set To clear a cell in a grid window use the menu item Clear current cell or Clear current selection after selecting the cells you wish emptied An exceptional case is when X or Xj is set to zero or one AND the other values T H and H fio are emp
111. low and the friction becomes independent of Reynolds number ATT eater OM aa 53 fo C 3 2 Two Phase Correction Factor The two phase multiplier 1 o friction liquid only is taken from Chisholm for the case of a rough casing with Blasius exponent zero when it takes the form rro 1 T 1 BrXaow 1 X ow Xgow gt where T is the ratio of pressure drop if fluid is single phase vapor to the pressure drop if the fluid is single phase liquid and is expressed as e P Hi Pv with the Blasius exponent n set to the value 1 4 for smooth flow at this stage since I will be used to determine the smooth parameter B The semi empirical coefficient Br for rough well casings is given in terms of the smooth B by 2 Bpr 0 5B p m poe Hi The value of B depends on and G as in Table 3 which is Table 7 7 from Chisholm 10 T G B lt 9 5 lt 500 4 8 500 lt G lt 1900 2400 G gt 1900 55 V G 95 lt TP lt 28 lt 600 520 T VG gt 600 2T gt 28 15000 I2VG Table 3 Vaues of B for smooth pipes C 4 Energy conservation Conservation of energy for a producing well takes the form see for example Chisholm X fow ds A gw Q dz Xfowhy 1 Xtow hi 5 uZ 7 u O gsind 0 where Q is the heat flow rate per metre of casing at depth z from the country to the well This is the total heat flowing per second through the area of the casing surface contained
112. lso read import input files from other SwelFlo projects or directories When the Save or Preprocess or Run buttons or menus are pressed the simulator saves the current simulation setup to the default input file The default name for the input file is SWFinput tat The Save command in the main menu creates this file in the project directory containing the input setup ready for reading if the simulator is directed to this project at another time The simulator can be directed to read any text file as an input file if desired The default name SWFinput txt is where all changes made are saved to The project structure associates each such input file uniquely with the Project folder it is inside There is no Save As command in the main menu If the user desires to save the current setup in the correct format to a different file name which will not be automatically overwritten they may copy that file within the main Windows desktop environment Having another file name 41 for an original input file is one way to protect that input file against changes being written to it by the simulator or to set up what is effectively a template since the simulator will not overwrite such a file unless it is called SWFinput txt Here is a line by line review of what is in the input file The line number is also indicated but should not appear in the actual file line 1 TEXT1 line 2 TEXT2 line 3 TEXT3 these lines of text are echoed in output files
113. ly consistent with the GWELL manual and Hadgu s thesis Note that the two phase correlations discussed in this manual have all been developed only for upflow of both phases and are incorrect for counterflow or downflow M is the mass flow rate kg s of vapour phase up the well M is the liquid mass flow rate Uy is the actual speed m s of vapour up the well M aApy Uy 1 with a similar eqn for liquid where a A A is the static void fraction or saturation A is the cross section area occupied by vapour p is vapour phase density and A is the total cross section area of the well in X tow gt 2 Q where Q M M is the total mass flow rate The velocity ratio is Ka Ul Slip velocity is u u and drift velocity is u uy where uy is the homogeneous velocity which both phase velocities are equal to if K 1 X flow LA un 6 E E E Pv Pl X gow K 1 A o uy 6 i 2 Pu Pl where G 2 Superficial gas velocity in m s is saturation times gas velocity and is M Usa D las SG pA u This is the steam velocity in the absence of the liquid phase Similarly superficial liquid velocity m s is M Usp 1 S SL pA yur gt 50 where M is the mass flowrate of the liquid phase and p is liquid density Void fraction a is used in Hadgu s thesis since he often deals with air water mixtures Void fraction is saturation S that is volume fraction of vapour phas
114. m dp dp oP MF 102 2 P dz dz Poe 102 where P is pressure Pa and 2 is the total pressure divergence This leads to dp dp dp dp dp dp dp MF E E 7 CL 2 N 2 a M where the subscript acc means momentum flux or acclerational pressure drop f means friction and grav means gravitational pressure drop so that MOMOM F 5 10 Gravitational Pressure Drop Duns and Ros use the same form as SwelFlo lt p gt gsin 6 F 6 Hagedorn and Brown Correlation Hagedorn and Brown made pressure measurements on a 457m test well for nearly 3 000 pressure points and a range of gas liquid flow rates viscosities and tubing diameters The key idea is to use a homogeneous equal phase velocity model to compute frictional and acceleration pressure changes and then to match the measured total pressure drop by choosing a value for the pseudo liquid holdup or equivalently a pseudo saturation that make the gravitational pressure change correct The dimensionless numbers listed in the Duns and Ros correlation are used to correlate the pseudo holdup Hence pseudo saturation is used as a parameter to correlate the equations with the measured pressure drop and saturation itself is never computed in this correlation It is noted in that this correlation is generally adequate for vapour saturations less than 0 99 0 95 but may not be accurate for mist flow 92 A major difference here to most other corr
115. m MD floating point Note that all depths in the simulation are actually distances along the casing that is measured depth not true depths The casing angles are used to find true depth changes when needed Feeds should be listed in order from shallowest to deepest The other values depend on what number is entered for FICODE I an integer They can be different for each feedpoint If FICODE is 1 Forcheimers equation the whole line should contain F_DEPTH I FICODE I RESV 1 1 RESV 3 D RESV 12 FDRY I RESV 6 RESV 4 1 RESV 7 I where RESV 1 1 is reservoir pressure Pa floating point RESV 3 D is reservoir temperature C RESV 12 1 is flowing fluid enthalpy J kg in the reservoir FDRY 1 is the floating point flowing dryness of the reservoir fluid at this feed RESV 6 I is the mass fraction of CO in the reservoir fluid RESV 4 I is the productivity index m of the feedpoint floating point RESV 7 I is the nonlinear term A floating point in Forcheimer s equation for drawdown at a feedpoint which should be set to zero for linear drawdown Note that enthalpy is interpreted as a flowing enthalpy from reservoir to well if provided This case means that the following Forcheimer equation relating pressure drop Pa at a feed to flowrate Q kg s is used 46 E ea 1 D a ll pet where P 2 v is the usual productivity converted from t hr bar to kg s Pa if A is zero
116. n can be written in terms of mass fraction by rearranging this to give 18XC YOu 18xG 44 1 XC A plot of XC against YC Fig reveals that XC is larger so that switching from the GWELL mass fraction formulation to Dalton s law will lead to smaller partial pressures of carbon dioxide given XC or to larger values of XC given a value for Poo 107 Note that Battistelli et al 2 describe a modified Equation of State module for the geothermal reservoir simulator TOUGH2 designed to handle three component mixtures of water sodium chloride and a non condensible gas Their algorithm for linking the partial pressure of the gas with the mass fraction of gas in the vapor phase is to assume densities combine as for an ideal gas so that total vapor phase density is the sum of the water vapor density at saturation pressure for the brine and the density of the non condensible gas at the partial pressure of that gas Then the mass fraction of the vapor that is non condensible gas is the ratio of the density of the gas to the total density of the vapor phase This leads to exactly the same expression for mass fraction as that derived above using Dalton s Law 0 2 04 0 6 0 8 1 YC vV Figure 41 Mass fraction of CO in the vapor phase XC versus the mole fraction in the vapor phase YC SwelFlo allows the user to choose to use the GWELL mass fraction formulation P XC XC 2 v P El or the Dalton s
117. n on downhole profiles 101 os S homog 0 10 20 30 40 50 60 70 80 USG Figure 39 A plot of the saturation obtained versus Usg value input to the Hagedorn and Brown correlation for P 40 bara S 0 9 Q 50 kg s showing more detail near zero Usa values Also shown are the pseudo saturation S and the homogeneous saturation which forms an upper limit on what the value of true saturation can be 102 G Choking As noted in for example 3 4 when the speed of a fluid flow in a pipe reaches the speed of sound any further increase in the pressure difference driving that flow has no observed effect on the flow rate and the flow is called choked Choking corresponds mathematically to a steady solution to the conservation equations ceasing to exist if one begins at bottomhole and integrates upwards to a depth where the speed of sound is exceeded by the velocity of a component Physically the cause is that in a time varying flow pressure differences propagate at the speed of sound so that if a flow is going at the speed of sound then any information about a change in pressure downstream cannot propagate upstream so cannot change the flow upstream Another view is that explained in Chisholm 10 that choking occurs when the entire pressure gradient is accounted for by the momentum flux term leaving no possibility of accounting for friction or gravitational terms The method used in SwelFlo is the same as
118. nce v is always positive so the area under the curve 1 v increases monotonically as the range of integration increases What is currently implemented in SwelFlo is to approximate the integral by using a twenty one point trapezoidal rule that gives good results and monotonic behaviour for flowrate versus pressure drop across the feedpoint Note that when flow is single phase and or the kinematic viscosity is practically constant this formulation drops back to the simple case 115 The pseudopressure approach provides a smooth combination of both single and two phase flow conditions 122 The left hand side of this equation may be understood by considering the concept of pseudopressure The quasi steady equation for isenthalpic flow in a porous medium is L8 p 292 k v Eve 0 121 V where v is in general the two phase kinematic viscosity 1 kri kro V VI Vy where k and kr are relative permeabilities for liquid and steam respectively and 1 and 1 are kinematic viscosities for liquid and steam respectively Equation 121 can be rewritten as Vim 0 122 where the pseudopressure m P is defined to be 1 m f Zap V Equation 122 then says that in general pseudopressure varies linearly with r in the steady state while pressure does not unless kinematic viscosity is constant so that Pros dP where C is some constant and R is the distance from the wellbore to the location that reservoir pressure is
119. nd Ros correlation Flow Region Limits bubble 0 lt Nyu lt Lo slug Li lt Nyo lt Ls transition Ls lt Ngu lt Lm mist Ngu Z Lm Table 6 Duns and Ros flow regions 84 After doing this we can find saturation S as follows using the definition of slip velocity in the form iy iy y 286 _ Use S 1 8 we rearrange this to obtain a quadratic equation for the static vapour saturation S in terms of known values Ugg Usz and V 93 VS V Usa Us1 S Usg 0 94 The solution of this that lies in the range 0 1 is given explicitly by AV5SU g gt 95 When V is zero this takes the limiting value Use s _ 96 Usa Use oe That is given Usa and Usz we can use the Duns and Ros correlation to calculate a value for the slip velocity V and we substitute these three values into eqn to find saturation S This defines the forward problem Finally we obtain liquid and vapour phase velocities by using Uy Usa S w Usz 1 S F 5 1 Bubble Flow Region Dimensionless slip velocity is given by N a Ven Fi FaNn Fl a 97 SD 1 alNiy Ge gt 97 where 7 F F A Na and the parameters Fi F2 F3 and Fy depend on the dimensionless liquid viscosity as determined by Duns and Ros and as illustrated in Fig 29 Digitisations of these graphs are given in Table 5 85 AAA X H ie Lt gut EA US A E 8 DI O E A PS EE D A E A PL ID a a a 9
120. ng on the flow region 89 Bubble Flow and Slug Flow For bubble flow and for slug flow Duns and Ros assume the continuous liquid phase is the main contributor and they correct the single phase pressure drop with their own two phase multiplier fm to obtain dp Vin LI t oU e 2 frapi D where Vm Usa UsL and their two phase multiplier takes the form p fm LA where f is the usual liquid phase Darcy friction factor that solves the Colebrook equation using the Reynolds number Re P iUsr Dh SY Mi and where f corrects for two phase effects and is a function of fiUs N A 3 4UsL plotted as the vertical correlation in Fig 32 and listed in Table B The parameter f is a further correction term for saturation and viscosity effects and is given by um am EE 4 a A mm i Figure 32 The value of f versus ETAF 4Us1 used in computing frictional pressure changes in the bubble region according to the Duns and Ros correlation 90 Note that even their liquid phase pressure drop term Vm Usa UsL We Depa ROT ale Pl SLop PIU SL 2D takes a different form to the usual SwelFlo term based on 10 which emphasizes the liquid velocity by weighting it with density G PUse pU sr PU se pUsz 2DhPm 2Dn Spv 1 S pr Transition Regime A weighted average depending on the value of N of the slug and mist flow values for frictional pressure d
121. ng steam quality to be determined by using phase velocities Given the criteria on saturation that Hadgu uses to distinguish flow regimes the forward problem can be easily set up with no iteration required to solve Then since in fact we have flowing steam quality and we need saturation values we call the direct problem iteratively to solve the inverse problem of finding saturation given flowing steam quality F 3 2 Bubble Flow Hadgu s equation 7 36 is obtained by using bubble rise velocities where buoyancy is balanced by a nonlinear drag term then modified to allow for variations across the radius of the wellbore and it reads 1 1 53 gApoy o Us Use z ET a 7 38 where is surface tension and Ap p pv In the forward problem we know at each step what the in place dryness or the static saturation is so from eqn 3 we know what the saturation a is and we know what Q is 71 total mass flowrate kg s so we rearrange eqn using G pw Usa where G 39 to obtain an explicit formula for the liquid velocity p 1 G 1 1 53 TA AV 0 A Oi oe ea 40 Us Po m2 Po a Tap e N Transition from elongated bubble to slug flow is taken to occur when a 0 3 that is for smaller void fractions we have bubble flow transitioning to slug as increases through the value 0 3 This gives the approximate condition 5 41 Us 1 67U sg 0 72 gt l sh 0
122. nthalpy of any fluid entering from the reservoir If a dryness is set at a feed point it is used to compute a static enthalpy which over rides any temperature or enthalpy value entered for this feed and is then taken to be the flowing enthalpy at this feed The effective two phase viscosity at a feedpoint is computed as an average by averaging the inverse viscosity 1 kri krv V n Vy over values taken between the reservoir pressure and the wellbore pressure where k 1 S and k S are relative permeabilities for liquid and steam respectively S is vapor saturation and 14 and v are kinematic viscosities for liquid and vapor phases respectively The computation of this viscosity is discussed in Appendix Jp It essentially assumes sub horizontal flow from the reservoir to the feed For a bottom up simulation there must be at least one feedpoint at a depth corresponding to the well length For topdown simulation a feed is not required to be specified at bottom hole or anywhere else since flow there is already predetermined by the rest of the well above it The Edit Row menu item in this grid window allows the user to insert or delete entire rows or clear their contents A maximum of ten rows or feeds is currently allowed The Clear menu item allows the user to clear the contents of the current cell the selected cells or the current row Feed depths are not required to match a subsection depth on input It they do not
123. ntsov s correlation to find the drift velocity without needing to know saturation Chisholm summarises it on his p 49 Uwp wWp 63 B 0 Equations lt 0 9 horiz 0 Can 4 59 is correct i E a Chisholm eqn vert 90 If Ug gt eTa then Cay Can USO 0 2 5 H Un lt ea then w 14 3 1 2 and g 1 v vv Uy i l If D 19 R then Cay f 1 53w Bu ae aar 17 1 else Cay E gD 1 v1 v v All Ca Can Cav Can ease patos Uy gt 09 Al Cy 1 22 ete 1 29 Table 4 Equations for the Armand Coefficient a corrected version of Table 4 3 in the GWELL manual and of Table 4 4 in Chisholm 10 9 is the angle to the horizontal in degrees non negative Where there is a sign the negative sign is to be used for downflow 0 2 5 m2 8 Pu Pl and the bubble velocity is computed as a slug velocity for small enough wellbore diameter or as the usual bubble velocity for larger diameters 0 25 1 53 enter DSD Up Pi gt 0 35 9D 1 p p D lt De where where the critical diameter is 0 5 D 19 10938775 55 g pr Pv and the value of 19 is refined here to 19 10938775 to give better continuity in up when D passes through the value D Typos in Chisholm include that the first line in his Table 4 4 should have P instead of the value 0 7 and eqn 4 72 should read 3 333 sin 1 80 0 333 sin 1 89 he has the sign
124. of T H Hilo X X 510 The simplest way to set one is to ensure the other four are cleared using the Clear menu item in this window There is a hierarchy used in the event that more than one of temperature enthalpy or steam mass fraction are specified and the most important term in the hierarchy is used to determine consistent values of all of the others and replace them with the consistent values if they are present otherwise show the consistent values when the user clicks the Pre Process or Run button 4 5 1 Wellhead Thermodynamic Setup Priorities Note that a distinction is drawn between static and flowing values of enthalpy and steam mass fraction which coincide for a single phase flow or for homogeneous two phase flow Any one of the five variables T H H fio X X fio is sufficient to determine the other four given P and XCO z The hierarchy is that 1 Xyio is most dominant a value between 0 and 1 will override any value that is set for T H X or H fio Next is Hfio which will override any value that is set for T H or X Next is X which will override any value that is set for T or H Next is H which will override any value that is set for T oa AeA Ww N The T value that is set will only have effect if none of H Hfi X X 510 are set 19 An exceptional case is when X or X yy is set to zero or one and there are no values set for T H and Hypo This is taken as a signal that the user wants to se
125. olution two phase liquid water with carbon dioxide in solution steam and carbon dioxide in gas form and vapor water and carbon dioxide both purely in vapor phase H 1 1 Liquid The mixture is all liquid if the total pressure is greater than the saturation pressure of pure water at the given temperature and the solubility of CO in water is greater than the given value of XCO 104 LIQUID bubble point curve Oo e e dew point curve 100 150 200 250 300 1 Figure 40 Phase diagram for water with 0 01 XCO2 The solid curves show bubble point boiling point and dew point condensation point while the dashed curves are contours of constant static steam mass fraction inside the two phase region between bubble point and dew point The dew point curve is only slightly higher than the saturation curve for pure water 105 H 1 2 Two phase If the total pressure is greater than the saturation pressure of pure water at the given temperature but the solubility of carbon dioxide in pure water is less than the given value of XCO then a gas phase will also exist and the fluid is two phase H 1 3 Vapor If the total pressure is less than or equal to the saturation pressure of pure water at the given temperature Ps T the fluid is all in the gas state H 2 Solubility of Carbon Dioxide The relative amounts of carbon dioxide dissolved in the liquid phase and present in the gas phase is determined by solubility Work by Ta
126. omhole for a bottomup simulation The best feed type for an output curve is Productivity Index Set this by clicking on the Feed Points button To get a match with the Topdown simulation it is necessary to set up the correct flowing enthalpy for fluid entering from bottomhole and the correct productivity index value These may be read from the RESULTS tab but should already match the values that have been read in from the template That is sigma should already be set to the value computed from the topdown run 0 45052E 1 and reservoir enthalpy should already match the value resulting from the topdown simulation 929 86 kJ kg To obtain a single bottomup simulation set the range of bottomhole pressures in the Output Curve section of the SETUP tab both to the single value obtained from the topdown simulation 48 505 bara Click RUN The graphs should show that the bottomup run matches the symbols from the previous topdown run as in Fig 4 If not check the work you have just done for mistakes It is a good test of numerical accuracy that topdown matches bottomup That is solving the differential equations describing conservation of momentum and energy in one direction then 11 SETUP DATA INPUT RESULTS GRAPHS ol lala P L v oc T deg C oc 170 175 181 185 190 Tos P Bara 10 15 2 5 30 35 40 45 Simulation O PDataSim006 txt Simulation O TDataSim006 txt Simulation O TDataSim006 txt
127. on 2 00 corrects the interpretation of the steam mass fraction that is returned by the subroutine CO2 based on the one in GWELL It is in fact a flowing steam quality that is returned when there is two phase flow in the wellbore not the static steam fraction Previous versions of SwelFlo which used the incorrect interpretation of a static steam mass fraction led 125 to enthalpies that were higher by up to about 5 than they should be Now the simulator works with flowing steam quality as a fundamental variable alongside pressure and temperature allowing correlations to be more simply computed then previously by starting with flowing steam quality and using the correlation to find static steam volume fraction saturation S Version 1 11 corrects a coding error in Version 1 10 that led to all feeds being treated as dry steam during output curve computation if any feed was set to dry steam type Version 1 10 adds a facility to invoke SwelFlo from the command line and to provide optional argument that cause it to run silently no windows are put up so an input file must then be provided that give an input file and that give a path to a folder that will accept output Version 1 00 is the first released version of SwelFlo It used some parts of routines originally from the Fortran 77 programme called GWELL converted to Fortran 90 but is also substantially improved and corrected and augmented Various errors in the coding of GWELL have
128. ooses the Baroczy correlation lt a 1 X ow 0 74 E 0 65 0 13 ee a a 9 X flow Pl Ly which we rearrange to find flowing dryness in terms of saturation a 1 C 1 CH J 65 ON o The way to implement this is to first compute Xgow from the above equations then combine equations and to get u from QX fow Uy aan 67 QU Y 64 B Then use Ugg au Then compute Ug and if it is less than or equal to one it is slug flow If Ugg gt 1 it is mist flow and Us G p Usapu m Implementation will run fast since it is explicit in saturation and does not need iteration to solve 75 F 3 6 Implementation of Hadgu Correlation When the results of the above correlation are plotted in the form of Usz versus S as in Figs and 26 a couple of difficulties are apparent The different regimes do not connect smoothly and the dependence of Us on S is not monotonic so that this function cannot be inverted unless the dependence is forced to be monotonic by introducing transitions between flow regimes 0 1 Usl 0 1 0 2 l 0 0 2 0 4 0 6 0 8 1 S Figure 25 An example of Us versus saturation for the Hadgu correlation with very small flowrate when G 54 T 246 C The leftmost brown curve is bubble flow the middle blue curve is slug and churn and the rightmost red curve is for mist flow There is a jump at the transition between bubble and slug flow which may cause issues when
129. or these different regions are discussed further below In each of these regions as detailed in the next subsection Duns and Ros offer a correlation that gives a value for the dimensionless slip velocity 0 25 Vas 1 2 92 go where the slip velocity is Vs Uy U To summarise then given the Duns and Ros correlation and given values for the superficial velocities Usa and Ug we are given a way to calculate a value for the slip velocity V 83 2 3 L La Na A R R Fi Ni Fs Fs Fi H f 20 05 10 0 1 25 0 24 0 83 5 0 0 003 0 218 0 48 0 125 0 001 1 0 2 0 0 55 20 0 1 25 0 24 0 86 9 0 0 005 0 205 0 2 0 11 0 01 1 02 1 97 0 74 30 0 1 25 0 24 0 98 16 0 0 007 0 195 0 01 0 098 0 1 1 05 16 0 9 40 0 1 25 0 24 1 25 220 0 01 O18 0 1 0 088 0 3 1 07 13 1 0 50 0 1 26 0 24 1 64 32 0 0 015 0 177 0 12 0 078 0 5 1 0 1 05 1 06 60 0 13 0 27 19 360 0 02 0 17 0 1 0 07 1 0 0 78 10 11 70 0 1 36 038 23 430 0 03 0 158 0 34 0 062 3 0 0 51 10 1 1 80 00 164 059 28 51 0 0 05 0 13 0 96 0 051 10 0 0 34 10 11 900 19 083 31 55 0 0 08 0 08 1 7 0 043 70 0 0 22 10 1 1 1000 21 0 95 3 25 55 0 0 1 0 058 2 14 0 039 200 0 0 208 10 1 1 150 0 2 05 1 03 34 550 015 0 045 2 16 0 034 500 0 0 208 1 0 1 1 200 07 19 0 98 3 7 55 0 0 3 0 053 1 86 0 03 1 65 0 93 3 8 55 0 0 5 0 067 1 76 0 027 13 085 40 550 0 8 0 085 1 74 0 026 0 95 0 75 4 1 55 0 1 5 0 105 1 74 0 024 Table 5 Digitized versions of numbers used in the Duns a
130. osen to give a match to total pressure gradient as observed in the Hagedorn and Brown experiments and as detailed here The Darcy friction factor fp depends on the Reynolds number Pn VmDh GD Ls hs Re where the two phase viscosity is defined in terms of the pseudo saturation as 1 8 Ss ksn The acceleration pressure gradient term used by Hagedorn and Brown is based on homogeneous flow dp _ ps d V dz J ace 22 dz Note that the minus sign is unclear in 86 section 5 3 due to the use of a A symbol but is correctly present in his coding 36 Appendix 8 p 149 Determination of the derivative of V2 in the above equation is tackled in different ways in different interpretations of the correlation One approach is that of Torrens A look at the coding in Torrens thesis reveals that there the derivative term above is evaluated by using 1 2 ya Pe 2 _ 77 Usa Usa Usi Usi where the superscript 1 means to use P to evaluate and the superscript 2 means to use P to evaluate Then the new value of Vm is given by V U a U B and the derivative is approximated as 2 1 y p vO MEVA yo AVe 2 dz Sem P ES Pi SAR AP l 94 Rearranging the Torrens calculation gives AV e aa Use Ap 105 p ANT py That is the Torrens approach is to consider that changes in vapour phase density dominate changes in Vm as pressure changes However Usa will also change and t
131. pecially with a different name so that it serves as a template Or simply change the project directory to a new one for a clean slate start with the template These changes are then written to the input file if the user presses Save Preprocess or Run 9 Output Files Generated Output files generated by SwelFlo are 1 The text file SWFoutput001 dat is the main output file always created in the current project folder It contains a a record of the input used to create it and b a summary of flowrates and properties at the wellhead and at the feedpoints the depth of the flash point and c detailed downhole simulated properties of the flow if a downhole profile was requested and d a summary of wellhead and feedpoint pressures and flowrates and properties and flash depth for each successful output curve point calculated if an output curve 35 was requested If SWFoutput001 dat already exists in the project folder then the name SWFoutput002 dat will be tried and so on up to SWFoutput999 dat before giving up and prompting for a filename from the keyboard This file is displayed if the user presses the Show Detailed Output button in the RESULTS tab in the main simulator window Columns output in part c of this output file are in order measured depth along the wellstring centreline vertical depth pressure temperature flowing steam quality the enthalpy of the liquid phase water plus carbon dioxide the en
132. r the Hadgu correlation with a larger flowrate when G 1718 T 246 C The leftmost brown curve is bubble flow the middle blue curve is slug and churn and the rightmost red curve is mist flow Note that calculations indicate that the S value where the slug correlation gives zero Usg does not vary much with G and is about 0 65 0 66 In the event that the mist value of U is always less than one in value there is no mist flow region and it is necessary to also ensure that no slug flow values are considered with S greater than the value which gives a zero value for Us S say We do this by using zero for Ugg and S 1 If the mist value of Ug is always greater than one in value for S in 0 32 1 there is no slug or transition region and we use the mist value of Us for the right hand end of the bubble slug transition region which will now be named the bubble mist transition region in fact F 3 7 Slug to Mist transition It is possible to compute an explicit value for S where Ug 1 0 for mist flow by combining and rearranging eqns 23 24 and 63 66 together with Usg au This leads after some algebra to the following explicit formula 1 S _ 1 0 68 where 1 C 0 74 C B a 2 Ca yg9DApp G The result of using these transition regions to remove jumps in Us is illustrated in Fig 27 TT bubble bu slug Figure 27 An example of Us versus saturation for t
133. raction solved for along with temperature in two phase regions and the flowing steam quality that appears in correlations momentum conservation energy conservation friction calculations and choking computations The following GWELL subroutines are affected by the confusion between flowing steam quality and static steam fraction ARMAND CHOKED ENERGY FEED1 FEED2 FEED3 FRIC1 MOMENTUM REGIME RESMOM and RESMOM2 which also has the wrong sign on the momentum flux divergence term A1 The change required to correct can be nontrivial For example in ARMAND the correlation essentially gives saturation in terms of flowing steam quality which when written correctly has to be iterated to get flowing properties given the current static properties In GWELL this subroutine is written in a contradictory and confused way by starting as if X is flowing not static steam mass fraction then computing saturation from that through the correlation but not properly iterating just repeating twice and taking an average if necessary and calculating a slip velocity where it should compute a drift velocity So correction of this subroutine involves a complete rewrite More details on this issue with GWELL and any code that is based on GWELL are to be found in 27 M Version Notes Version 2 01 has a minor fix to correct the output of flowing wellhead enthalpies to the file SW Foutcrv csv that contains output curve simulated results if any Versi
134. re any computation of pressure gradients within the correlation Such computation is left to the numerical computation of conservation equations elsewhere in the code F 7 The Forwards Hagedorn amp Brown Correlation The Brill and Beggs 1991 7 interpretation of the Hagedorn and Brown correlation is considered as described in 7 36 I together with the original publication 21 A homogeneous flow model is used to compute frictional and acceleration divergence of momentum terms and the pseudo holdup or saturation value to be used in these terms and in the gravitational term is determined by the correlation itself This pseudo holdup is essentially a correlation parameter that Hagedorn and Brown have fitted to their data that gives a match to experimental pressure gradients for a given pair of values of steam and liquid flow Usa and Us The conservation of momentum equation is written in the following form in terms of the total pressure gradient dp dp dp dp aP O E les of 103 7 E A E grav i where subscript acc means accelerational pressure gradient associated with the momentum flux divergence subscript f means pressure gradient due to friction and subscript grav 93 means gravitational pressure gradient The friction pressure gradient is dp V2 iG m_ 104 where 2 Ph PiUst psa G j j c Spt i s Vm Us Use Pa Ps y UsL Use Vin p Je ds did and S is a pseudo saturation ch
135. rmand s correlations for various values of Q The line S is included for comparison in the S vs P plots The data points are Armand s experimental data scanned from Chisholm s book The solid line is the correlation given in Table 4 The dashed line is the simpler correlation C A Cn Pressure is set to one bar and wellbore radius to 26 mm as in Armand s experiments However the most straightforward fix for the jump in the Armand correlation suggested by the good fits to the Armand data is to use the simple formula C4 Ca given in equation 4 59 in Chisholm which gives S a 25 65 f S 0 0 02 04 06 08 1 0 1 03 05 0 7 0 9 0 02 04 06 08 1 01 03 05 0 7 0 9 b Q 5 c Q 10 4 Q 10 0 5 0 0 5 0 02 04 06 08 1 0 1 03 05 0 7 0 9 0 02 04 06 08 1 01 03 05 0 7 09 B B B B e Q 20 1 Q 20 8 Q 50 h Q 50 o o a 3 02 04 06 08 1 T ol 03 05 0 7 0 9 03 T 0 4 06 0 8 1 0 1 03 05 0 7 09 i Q 100 j Q 100 k Q 200 1 Q 200 Figure 20 Plots of saturation S BC4 Sp 6 versus volume flowrate 6 and of the Armand coefficient C4 versus 8 for Armand s correlations for various values of Q The line S is included for comparison in the S vs P plots The data points are Armand s experimental data scanned from Chisholm s book The solid line is the correlation given in Table 4 The dashed line is the simpler correlation C4 Cap Pressure is set to 10 4 bar and wellbore r
136. rop is used here That is if DP is the slug term and DP the mist term then A pP gu Hs DP DP A Mist Flow Regime Friction in the Duns and Ros mist flow regime is assumed to be due to the drag of the vapour on the wavy liquid flow that clings to the walls of the wellbore so that dp z Usa a fibazp gt 99 where the Darcy friction factor f depends on the vapour phase Reynolds number _ PyUscDr Hv Re The roughness this flow experiences is associated with ripples on the liquid film on the walls They are correlated to the number N defined as 2 Ny Po o Pl so that 100 E ED Na lt 0 005 D 17482 N gt 0 005 pgl2cgDn If these relative roughness values are found to be greater than 0 05 the friction factor is found using the extended Colebrook equation 1 EN LT 4 E log10 0 270 D J 0 067 5 101 91 fi F 5 9 Momentum Flux Terms Duns and Ros take the pressure drop due to changes in momentum flux to be zero unless the flow region is mist flow SwelFlo adopts this approach if the user selects the full Duns and Ros correlation option However our experience with other correlations suggests that momentum flux changes can be significant in accelerating slug flow regions Then in the mist regime the pressure gradient due to divergence of momentum flux in the mist flow region is taken to be a given fraction of the total pressure derivative dp Usa Ust Usap
137. s Run Stop 4 4 2 2 a 46 A A eS Ge 35 9 Output Files Generated 35 10 Command Line Execution 37 38 A Input File Format 41 B Nomenclature 49 C Governing Equations 50 Col Two Phase Flow Terminology 2004 84 4 s 2 BS a en B Hd HEY 50 C 2 Mass Conservation of Water and Carbon Dioxide 51 E A Bae AR O en ees See Pe ee 52 C31 A AAA 53 C 3 2 Two Phase Correction Factor a oe Bee eR Owe e e eS 54 pA Energy conservation be ta ee eb Re Oe ee EEE OS ESE oH ewe 54 C 5 Heat Flow between country and casidgl o o 55 D _Discretisation and numerical method 55 D 1 Newton Step rre ar eee saes 57 D 2 Modified Newton with Armijo Rule o 58 E Conservation equations at Feedzones 60 F Correlations 62 EJ Armand correlation daa es wa Glad Be OER AEE Ow RE ME RES 62 F 2 Orkiszewski Correlation f2 d e 08eeee eane bteaeebhbdde ses 67 F 3 Hadgu Correlation hh ee he ee ee ee A BES 71 F31 Phase Velocitie o s roes so oe DRE Oe oe a OS be 71 ERE OW 2 4 8 aae Se He hk HR oH AG ee i 71 Pad Slug NOW cae eke Bad WE Dae WP eee ee ee Oe ee a x 72 F 3 4 Churn II 75 F 3 5 Annular or A a 75 F 3 6 Implementation of Hadgu Correlation 76 F 3 7 Slug to Mist transition 2 24 44 s arre sens ee Ae eae TT F 3 8 Implementation logic a o a a o o 78 F 4 Frictional Pressure Drop Full Hadguj 79 F 4 1 Bub
138. s for varying v P with a pressure gradient dP dr across a narrow region of practically constant pressure in the form of Darcy s law kA dP v P dr 116 where k is an effective permeability to flow W is the mass flow positive if flowing outwards from the well and A is the cross sectional area the flow is passing through Rearranging this and using Q W so that positive Q corresponds to production and integrating from undisturbed reservoir to wellbore gives dP Q dr 117 l v P kA oot In a linear one dimensional geometry equation 117 becomes P res dP R f a 118 Pao UCP kA and putting Y kA gives P res dP Q yl 119 Pwen v P and comparison with 115 indicates that the effective viscosity should be written as 1 1 Pes qP f ss 120 Var Pres well J Por v P This is in the integrated sense the average value of 1 v when plotted against pressure The same integral average formula for productivity index arises if a radially symmetric cylindrical or spherical geometry is used instead of linear The simplest approximation to this would be the usual arithmetic average of 1 v using reservoir and wellbore values This is however not very accurate in some situations it can give non monotonic flow dependence on the pressure difference Pres Pwen In particular the form of equation 119 makes it clear that flowrate Q should be monotonic in the pressure difference Pres Pye si
139. ser to re plot any existing graphs of simulator results and data files taking into account any changes recently made since the last simulation run to the ranges of flowrates desired in the output curve and taking into account any changes made to the data files the user wants displayed with the simulation results or any new data files added to the list the user wants displayed alongside the simulated results Re plot uses the most recently generated output file for simulated values hence if the user has changed units without re running the simulation plots will show incorrect units upon replotting 4 10 Output Curves An output curve is produced if Yes is selected in the radio buttons under the Output Curve label in the SETUP tab of the main SwelFlo window see Fig 7 The number of points desired may be set just below these radio buttons This number is 26 indicative only If it was zero and the user selects Yes for an output curve this number defaults to 50 It may be adjusted by the user The range of flowrates desired at the wellhead is used if nonzero to restrict the plotted points presented in the output curve The range of bottomhole pressures set by the user are to help the simulator find those bottomhole pressures that give steady state solutions for the current setup If these pressures are set to zero the simulator defaults to searching between one bara and the reservoir pressure at bottomhole This range is
140. servoir is consistent with the fluid before the feed even if it is not set up that way A warning is issued Cases One and Two are OK since the feed injects the wellbore fluid into the reservoir and there is no physical contradiction if the injected fluid properties do not match the reservoir fluid properties F Correlations The algorithms used in SwelFlo to implement the correlations are described in more detail in this section The original correlation publications were consulted and coded from scratch The simulator needs values for vapor or liquid phase velocity given pressure temperature or enthalpy gas fraction and flowing steam quality Saturation or equivalently vapour phase velocity is provided by the correlation Correlations are a short cut accounting for the effect of interfacial terms in the momentum and energy conservation equations when they are written in terms of the individual phases then added together to get simpler and fewer equations for total amounts F 1 Armand correlation The Armand correlation requires a knowledge of P to compute void fraction a from the Armand coefficient C4 using a C43 where 8 is the volumetric gas flow ratio volume of 62 vapour per second divided by the total volume flow rate which can be written as j 1 Xfow Pv aa P o S 0 Then the velocity ratio is calculated from NETO 0 Q e K 1 in A py pl Then Uy and u Ue E give the phase velociti
141. simulation Typically the simulator is run with just the first row active a pressure value less than one bar in a row will cause the simulator to ignore it For example you can change the wellhead pressure to 7 bara and Run again There is no need to Save any changes they are automatically saved every time the simulator is Pre Processed or Run The Save commands on the grid windows are really to allow the user to save for example the wellbore geometry to a file that can be opened later by the same kind of grid window in a different project To check that correct conditions are set at the wellhead and at feedpoints and that the simulator is ready to run click the Pre Process button in the SETUP tab window Clicking the Run button also performs a Pre Process to check that the simulator is ready to run before running it Besides the quick plots of simulation results text values of simulated variables versus distance down the well are available under the RESULTS tab A summary of wellhead values flash depth and feedpoint values is presented plus a button to browse the latest output file produced by SwelFlo The first parts of that file describe the input setup used so that this output file provides a complete summary of the simulation performed This output file is read by SwelFlo to produce the RESULTS summary and the plots displayed so that if a user edits the latest output file it is likely that SwelFlo will be unable to show RESULTS and p
142. sity may be written down 0 25 Ni Ust 2 gt 85 go pi 0 25 Ny U 86 w Uso L 86 WN Dig 87 O g 0 25 N n 88 a i 88 where D 4A P is hydraulic pipe diameter which matches the actual diameter for a circular cross section A is cross section area and P is the length of the pipe perimeter Note that a variety of sources 29 9 36 11 7 confirm that the density in N and in the other dimensionless numbers is indeed liquid density p Then the flow regime is determined as bubble slug mist or transition according to the values of the variables Ly Ls Lm Ly Lit Lays 89 82 L 50 36N 1 90 Lm T5 84N 2 91 v The numbers L and La are correlated to the nondimensional pipe diameter as illustrated in Fig 28 and this correlation is evaluated in SwelFlo by using interpolation on a digitized version of this figure as given in Table 5 Figure 28 The values of L and Ly and their dependency on diameter number N4 according to Duns and Ros The flow regions bubble where the liquid phase is continuous slug where neither phase is continuous with churn flow at low liquid flow rates and dispersed bubble flow at high liquid flow rates transition which is where there is a transition to mist flow and mist where the gas phase is continuous and annular liquid flow occurs are delineated by the values of Lp L and L as listed in Table 6 The correlations f
143. system in any form or by any means without the prior written permission of the copyright owner Version 2 02 provides for the display of saturation pressure of pure water dashed line on the same graph as pressure solid line This gives an indication of how far above the bubble point of pure water the fluid is which echoes the temperature profile and the carbon dioxide concentration Known Bugs When the user was allowed to direct the simulator to open and read from a saved Feedpoints comma separated file into the Feedpoints Grid window it crashed Hence Open and New have been removed from the File menu for this grid window for the time being In a future release an Open command will be added separately to its menu 2 Introduction The geothermal wellbore simulator SwelFlo is intended as a tool for geothermal reservoir engineers and modellers and for academics investigating two phase upflow in pipes to allow them to rapidly set up and solve a model of a geothermal wellbore and to test that model against data from production tests injection tests and output tests The governing equations are as described in Appendix Ch They solve coupled steady state mass momentum and energy conservation equations for liquid vapour or two phase water with carbon dioxide gas present in both phases The flow correlations available to determine the relative velocities of vapour and liquid phases are the Armand correlation the Orkizewski correlation the Ha
144. t capacity and time since flow began may be set or altered here also The choice of time since flow began reflects the transient heat flow model that is used to compute heat flow It would usually be the time in seconds since well discharge commenced starting from an equilibrated zero heat flow initial state with the temperature in the wellbore assumed equal to the temperature in the country at each depth value A typical value is one week or 604800 s 4 12 Depth Spacing for data output In the SETUP tab in the bottom right hand corner is a box labelled Spacing for data output m This allows the user to choose the depth spacing between simulated results output to text files in the project directory for pressure temperature and velocity These files are in the correct format for input as data files to SwelFlo say from later runs in the same project folder The files have names that include the project folder name together with a three digit generation number that matches the number on the main output file and this allows the user to compare the results of simulations in different projects or between different simulations in 28 Rock Properties Allow Heat Flow Yes To From Country No Thermal Cond 20000 W K m Rock Density 2800 0 kg cubic m Heat Capacity 1000 0 J kg K Time Since 0 60480E 06 Flow began s Spacing for gt data output m 50 g Figure 12 A view of the Rock Properties section of the SETUP
145. t high temperatures and pressures AM J Sci Vol 262 1964 p 1055 1074 34 Taitel Y Advances in Two Phase Flow Modeling SPE27959 1994 presented at the Tulsa Centennial Petroleum Engineering Symposium held in Tulsa OK USA 29 31 August 1994 35 Todheide K amp Franck E U 1963 Das Zweiphasengebeit und die Kritische Kurve im System Kohlendioxid Wasser bis zu Drucken von 3500 bar Zeitschrift fur Physikalische Chemie Neue Folge 37 pp 387 401 36 Torrens T S Vertically upward two phase flow in an annulus a thesis submitted in 1993 for the degree of Master of Engineering Dept of Mechanical Engineering University of Auckland New Zealand 37 Yadigaroglu G amp Lahey R T Jr On the various forms of the conservation equations in two phase flow Int J Multiphase Flow Vol 2 1976 pp 477 494 A Input File Format Here is a full description of the file SWFinput txt It can be generated directly by SwelFlo so that the user does not need to create an input file It contains all of the wellbore parameters and simulation settings that SwelFlo needs to determine a simulation and is the main file read in by SwelFlo when it first enters a project directory It is where SwelFlo saves its current setup It is a text file so can be modified directly using a text editor if desired The advantage of using the simulator to change the input file is that the format is then guaranteed correct The simulator can a
146. t the fluid to the bubble point zero or liquidus one but only if all of the other values T H and Hy are Cleared If both X and X fio are set to zero this means the same thing If there is a conflict X fio wins It is never necessary to Save the Wellhead Conditions grid window or any of the other grid windows in SwelFlo All simulation settings are automatically saved to the current project folder each time Pre Process or Run or Save in the main window menu is clicked and if a user tries to exit after altering a simulation setting SwelFlo will prompt for an overall Save before exiting The overall save writes to the file SWFinput txt over writing any previous setup in the current project directory This file is the main place that a simulation setup is saved Since it is over written each time the user modifies the project it is useful to copy the SWFinput txt file to another name using the Windows environment if a user wishes to preserve that setup say as a template for other simulations Such a template may be imported to any current project over writing it after a warning using the menu in the main SwelFlow window Project Import Existing SwelFlo setup It is possible to Save the Wellhead Conditions setup from the Wellhead grid window menu Then a csv file is created with a default name SWFWHcon csv that can be Opened by another Wellhead Conditions grid window providing the user with a way to share parts of the simulator
147. tab in the main SwelFlo window the same project for example different correlations or different feed parameters used The units used when saving are the currently chosen units set by the user for pressure and depth and flowrate 5 Data Input Tab The DATA INPUT tab pictured in Fig allows the user to add data files for display on the same plots as the simulated P T V versus depth profiles and the output curve if asked for It allows the user to choose files generated by previous simulations in the same folder or elsewhere or to create and then edit completely new data files from scratch Up to twenty data files may be read edited saved or created to go into the well folder containing either pressure temperature or mass averaged fluid velocity values versus depth downhole or output curve values wellhead pressure vs flowrate The format that these files must have is summarised in the Data Input tab and is described in detail here Downhole profile data files for P T or V have either pressure temperature or velocity data versus depth MD and the first line should be Mass flowrate at wellhead wellhead pressure flowing enthalpy and the remaining lines up to 10 000 of them should each contain depth MD data value The mass weighted velocity is the combination Xu 1 X u 17 where X is static steam dryness u is vapor phase velocity and w is liquid phase velocity The units for pressure and depth in data
148. taken from 123 K A Sample Input File This is an example of the file SWFinput txt that is used by the program to save the simulation input setup between simulations It is kept in the project directory and is read by the simulator on startup and whenever the project directory is changed It is saved to whenever a simulation is run It is not intended to be edited directly by the user but indirectly by SwelFlo However it is just a text file and can be edited by the user if desired A template used in absence of an input file Prodn flash in well One feed topdown flow Orkiszewski correlation No output file asked for No reservoir T profile one casing section 1000m vertical Ready to run 1 1 0 O 0 000000 OTFTTF 0 00 0 00 O 50 O No data files 2 wellhead conditions 800000 0 20 0 169 75 837557 42 920000 0 5 955E 002 0 10 1 0E 03 0 0 0 0 0 0 0 00 0 00 0 00 0 0 0 0 1 00000000E 003 0 0 0 0 1000 000 2 00 2800 00 1000 00 604800 0 1 casing section 1000 000 0 100E 00 0 00E 00 0 200E 02 90 0 0 00 1 O no reservoir profile 1 feed point 1000 0 1 5500000 1 216 916 929860 0 0 0 1 0E 003 4 5052E 012 0 0 L GWELL coding errors This is just a short note on some issues with the existing coding of the Armand and Orkiszewski correlations in GWELL which also affect the computation of energy conservation in GWELL and in codes based on GWELL These issues arise around the meaning and use of dryness In the GWELL code dryness X is used in
149. terchangeably in 1 thermodynamic calculations where it reflects the void fraction and is dryness for a static fluid or for a homogeneous flow with equal velocities for vapour and liquid phases and 2 in correlations and in the energy conservation equation where it must be the flowing steam quality or the ratio of the mass flow rate of vapour to the total mass flow rate 124 The correct dryness to use in the correlations and in the Orkiszewski correlation as noted explicitly in the GWELL manual is Xgow the ratio of steam mass flow rate to total mass flow rate the flowing steam quality This mass fraction is a flowing mass fraction different to the in place mass fraction whenever vapour is flowing at a different velocity to liquid Dryness and enthalpy measured at wellhead will in fact reflect this flowing mass fraction But the thermodynamics of vapour liquid carbon dioxide equilibrium depend on saturation or mass fraction in place GWELL s correlations and energy balance assume that the flowing dryness is the same as the dryness in place which is generally incorrect This manual describes corrected Orkiszewski and Armand correlations as well as a more careful energy balance that reflects the flux of energy into and out of a computational element rather than the existing GWELL energy conservation equation which uses energy in place multiplied by total mass flow In GWELL the variable X is used to denote both the static steam mass f
150. thalpy of the vapor phase water plus carbon dioxide the velocity of the liquid phase and of the vapor phase the density of the liquid phase and of the vapor phase water plus carbon dioxide total static enthalpy of the fluid radius flow regime name the total mass fraction of carbon dioxide the total mass flowrate the flowing enthalpy of the mixture the enthalpy of the water component in the vapor phase the steam the enthalpy of the water in the liquid phase the enthalpy of the carbon dioxide component in the vapor phase the enthalpy of solution of the carbon dioxide dissolved in the liquid phase the partial pressure of carbon dioxide the saturation pressure of pure water at the current temperature and the word CHOKED if the flow is choked in this section the comma separated file SWFbplot001 csv is created in the working folder if a downhole profile was requested with numbering in the name being the same number as that currently applying to SWFoutput dat It contains comma separated simulated results depth MD downhole P T C and mass weighted velocity m s positive for production each in its own column versus depth MD plus any data provided in data files This file is intended for graphing when the user wants total control over the appearance of the graphs and is ready to be used as input to the Grapher package or opens very nicely in Excel Each column has a header describing what it is depth
151. the bottom feed may be of any type productivity index fixed flowrate dry steam or drawdown 4 3 Total well length There is a place to enter the total well length arc length along the casing centerline just above the button for Pre Process Run Stop and Re Plot The reason for this is to provide a check on the casing geometry setup If the sum of the casing sections provided by the user in the Casing grid window is different that sum over rides the value entered by the user and is displayed there In the event that a top down simulation has no feed information provided the total well length is used as the default location of the bottomhole feed A bottom up simulation requires at least one feed to be provided in the Feed Points grid window 4 4 Correlations The conservation equations that are solved in the simulator are simplifications of the conservation equations for the individual phases liquid and vapor phases of water obtained by adding them together to simplify interfacial interaction terms in the original equations As a result extra empirical conditions are needed to determine the relative velocities of liquid and vapor phases depending on saturations They form a correlation between static steam saturation and flowing steam quality or equivalently between static steam saturation and the flowrate of one of the phases of water A number of these conditions are in the literature and the ones currently on offer in SwelFlo
152. ts a folder When the user makes changes to a simulation setup these changes are written to SWFinput txt whenever a Pre Process Run or Save is performed 33 The user is only prompted if there are unsaved changes in the original Project directory which will be lost when switching to the other directory The input file selected in the newly chosen pre existing directory is then read into internal arrays and any graphical data is also read in and displayed 8 1 5 Save This item causes the simulator to save the current setup to the default input file SWFinput txt essentially saving the project ready to be read in again if the simulator is restarted or switches to this Project directory at a later time If there are preprocessing errors the user is prompted before the input file is written to There is no Save As menu item each input file is uniquely identified with a project folder so that output files etc are uniquely identified with the input file that generates them Of course the user may choose to duplicate by Importing a SWFinput txt file from another project And in the Windows environment a user may save a SWFinput txt file under another name to preserve it for use as a template ready to be imported by using the menu item Import Existing SwelFlo setup But a Project folder has only one working input file associated with it and saved inside it and it is called SWFinput txt The Save button is greyed out and unresponsive if
153. tual directions of flow rather than the conventional directions for positive flow values In Cases Three Four and Five conservation of carbon dioxide across a feed when the overall mass fraction is XCO2 takes the form Qn4 1XCO2n41 QnXCO2 BE Qs XCO2 gt where the subscript f refers to the properties of fluid flowing between reservoir and wellbore and the flowing fluid enthalpy takes exactly the same form in these three cases Qa lt h gt n41 Qn lt h gt n Qg RA and all enthalpies are flowing enthalpies rather than static enthalpies No distinction is made between flowing and static enthalpy for any reservoir fluid entering at the feed The fluid in node n 1 is then reverse engineered as necessary to find a new temperature that corresponds to the new values of flowing enthalpy and XCO2 While no account is taken of the kinetic energy of any incoming reservoir fluid the kinetic energy of the up flowing mixture is accounted for in the next energy balance In the other Cases One Two and Six the values of lt h gt and XCO2 must all be the same before and after the feed For a topdown or a bottom up simulation this makes Case Six problematic unless the user ensures that the values before the feed match the reservoir values because there is otherwise a contradiction between the required values from the reservoir and the actual values set The simulator continues past such a case by assuming that the fluid entering from the re
154. ty This is taken as a signal that the user wants to set the fluid to the bubble point zero or liquidus one but only if all of the other values T H and Hy are all empty If both X and X yy are set to zero this means the same thing If there is a conflict X fi wins A simple way to ensure the desired value is being set from the list T H Hyto X X fio is to clear the other four values Then Run Clicking Pre Process will populate the remaining four with consistent values 3 2 SETUP Other settings in this tab include Produce a downhole profile or not turning this off will ensure no downhole values of P T etc are produced or presented Simulation Direction can be Topdown or Bottomup Wellhead values determine the starting point for a topdown simulation and fluid properties at the lowest feed which is also taken to be the bottom of the simulated wellbore follow from the simulation results Bottomhole values reservoir values of P T etc are used in a bottom up simulation with the exception of bottomhole pressure which is varied in an attempt to match wellhead pressure and to a lesser extent to match flowrate Reservoir properties as well as feed properties at a feedpoint including bottomhole may be set using the Feed Points button and entering values into the grid window that appears An output curve may be requested and the number of points desired in it The actual number achieved will often be less than the num
155. velocity of the gas in Taylor bubbles is a Ucrg Un Un Vers QTB the ratio of Taylor bubble length to slug length is j Usa ALsUaLs 43 arBUGTB and the volume average void fraction is asu Parsg 1 Bjazs 44 Equation 43 given that azs 0 3 and arg 0 89 gives P as a function of the superficial liquid and vapour phase velocities Then equations through to 44 provide a recipe for going forwards from superficial liquid and vapour phase velocities to saturation S QSU 73 The forward problem to be solved is to find the phase velocities given static dryness X pressure temperature and total flowrate The flow Q may be used as in eqn to eliminate Usa Then eqn 43 through its dependence on Uy in numerator and denominator gives as a ratio of linear functions of the superficial liquid phase velocity Usz The known in place dryness X gives an independent value for volume averaged void fraction or saturation X pi Xpit 1 X po Asu 45 Combining eqns and allows to express in terms of in place dryness X Then eqn can be rearranged after some algebra to give an explicit expression for Us for slug flow as a function of dryness X or saturation S This result is summarised below BAs Ag h 4 Us L29A7A pr where 46 Ag ALS Pu Un ULA UrrgAs G 1 1 29az s41 47 As 1 29G arg arsAs argUsopulAr A243 48 Upprars UrrBpvaLs
156. w be called with command line arguments from a Windows command window The arguments that have effect are command followed by effect S silent mode no windows f full path and filename input file required for silent mode a path for directory output is to p path name go to if different to directory con taining input file Note that silent mode requires an input file to be provided The name of the input file must contain the full path If spaces are desired in this name enclose the whole thing in double 37 quotes There must be a gap between the f and the name and between the p and the path name if a separate path for the output directory is provided For example the command gt swelflo exe s f c swelflo testruns p1 SWFinput txt will run SwelFlo silently on the input file SWFinput txt in the folder c swelflo testruns p1 and output will go to the same folder The command gt swelflo exe s f c swelflo testruns project one SWFinput txt p c swelflo testruns out1 will run an input file in the directory c swelflo testruns project one but output will go to the folder c swelflo testruns out1 Note that the directory that output is to go to must exist before running SwelFlo If an input file is provided but silent mode is not invoked SwelFlo will run in full interactive windows mode using that input file and running in the same folder as the input file ignoring any output path provi
157. wn correlation j 1 0 0 17536897977Nsec Nsec lt 0 01 0 014 0 006 n 0 0044 0 002 0 0014 T T T T T 0 005 0 01 0 05 0 1 0 5 l N l Figure 34 A graph of the function C versus N used in the Hagedorn and Brown correlation Finally we have Na ies Nea and the dimensionless numbers Nw N Ni Na are as defined in eqns 85 to 8S in the previous section on the Duns and Ros correlation Ns Note that in 7 the definitions of these numbers appear at first sight to be missing the gravitational terms these are absorbed into the coefficients used there The different conversion factors used in 11 7 arise because surface tension is in different units to the other parameters petroleum engineering units have been used and because gravity has been absorbed into these factors also 98 0 Ol 02 03 04 05 06 07 08 09 410 sec Noy N 380 1 No 1 fe 0 0 01 0 02 0 03 0 04 0 05 0 06 0 07 0 08 0 09 ai Lt N Figure 35 A graph of the function w versus Nsec used in the Hagedorn and Brown correlation Plotted alongside of this for comparison is the corresponding figure from the original publication of Hagedorn and Brown F 7 2 Implementing the Hagedorn and Brown correlation Fig shows a graph of the pressure gradients typically obtained versus a range of values of Usa as input to the forwards Hagedorn and Brown correlation as described in
158. wrong on the second term and in his Table 4 4 the formula for C4 for all inclinations when P lt 0 9 should read Ca Can Cay Can a aa 64 that is the first minus sign in Chisholm s formula should be a plus sign or the subscripts on C4 should be reversed in the first set of parentheses However simulations sometimes failed to converge when using the Armand correlation as summarised above that is they sometimes failed to find sufficiently accurate approximations to the conservation equations This is due to the discontinuous nature of the correlation which is illustrated in Fig obtained by solving the Armand correlation to relate S and for a well of radius 0 026m at a pressure of about 1 bara and temperature 100 C for various total mass flow rates It is the jump at 6 0 9 that gives the trouble This is where Chisholm suggests switching to the Zuber and Findlay correlation for annular flow and there is no attempt to ensure this switch is a continuous one A possible fix would be to smooth the region where 3 gt 0 9 with a straight line 0 0 5 al 0 5 0 02 04 06 08 1 01 03 05 07 09 0 02 04 06 08 1 or 03 05 07 09 a Q 5 b Q c Q 10 4 Q 10 07 0 0 5 0 02 04 06 08 1 0 1 0 3 05 0 7 0 9 0 02 04 06 08 1 0 1 03 05 0 7 0 9 B B B B e Q 20 1 Q 20 e Q 50 h Q 50 Figure 19 Plots of saturation S BC4 Sp 8 versus volume flowrate 6 and of the Armand coefficient C4 versus 5 for A
159. x term relies on A being constant since he writes OVe Oz which relies on A or equivalently G being independent of z to bring it outside the derivative P The momentum flux term in GWELL also is only correct if A is constant with depth since GWELL uses the term 2 Gu instead of G2 um Both GWELL and Chisholm s versions valid for constant area are incorrect by about a factor of two at places where wellbore radius changes Such places are also where the momentum flux term can dominate the pressure gradient so the effect on pressure gradient there is significant However the effect of using the incorrect form on overall pressure profile down a well is typically still small 25 The Hagedorn and Brown correlation uses a momentum flux term of the form dP Ps A 2 SE A E 2 dz m gt Vm Usa Us Su 1 Sju Ps S Po 1 3 p1 where 119 This assumes a constant area It is not obvious how to include the effect of area depending on z since the correlation uses homogeneous flow to compute Vm then corrects using S Note that the correlation is not used to compute momentum conservation directly only to calculate vapor and liquid phase velocities given temperature and saturation Hence no changes have been made to it The full Duns and Ros correlation assumes zero momentum flux except in the mist regime where flow is taken to be homogeneous so that liquid and vapor phase velocities are the

SwelFlo User Manual Version 2.01

Contents

Download Pdf Manuals

Related Search

Related Contents