Programmer’s Manual .pt



Programmer’s ManualThis document describes the VBA code to solve the 1D evolution equation using the case of a river as an example. A river is a 1D problem that includes a free surface through which heat and gas (mass) can be exchanged. A river is also a system where mass and heat can be released through lateral discharges and finally a river is a system that every student knows and thus where transport processes are easy to understand.The code is written to allow the simulation of several properties and to be easily extended to simulate other properties. The code follows a rational identical to that used to simulate heat conduction (diffusion) along a fin or across a solid body used in the first part of the “Energy and Mass Transfer” course. It reads data from excel worksheets, writes the results on the same worksheets and plots the results in a specific worksheet.Data is read using keywords to easy input data specification. Data defining the river, the atmosphere and the numerical options are read from worksheets: “River”, “Atmosphere” and “Main”. The “Generic” Property is transported and can decay (or grow) at a constant first order rate. The computing parameters are specified in the worksheet “Generic”. Other properties are specified one per worksheet, Temperature, DO, BOD, etc. The graphs are plotted in the “Graph” worksheet.The theory supporting the global transport model, including the upwind algorithm is provided in the document “Evolution Equation.docx”. Other documents describe additional methods and issues specific of each property (temperature, DO, BOD).Requirements for the ApplicationThe application will be used as a numerical laboratory. For that reason, it must:Allow the simulation of several properties and several discretization methods: upwind, Central-Differences and QUICK for space and explicit, implicit or semi-implicit for temporal discretization,Input must be flexible and able to detect errors,Standard graphics must be done automatically to easy results assessment,The structure of the application must be simple and easily upgradable to facilitate the inclusion of new features by students and to give that some information about programing practices.The VBA language was chosen to implement this application because it allows the use of excel worksheets to organize input and output and because it allows the organization of data into data structures (Types) that easy the inclusion of new properties in the code. The data types used in this application are identical to those used in the “Fin Application”.The main data structures are listed in REF _Ref39157682 \h \* MERGEFORMAT Table 1. The first column lists the name of the structures (Types) and the second column gives an overview of the data included in that structure. The content of each structure is listed in the next paragraph.Table SEQ Table \* ARABIC 1: Main data structures used in the “River Application”. On the left column is listed the name of the structure and on the right is a small explanation of the content.Type GridSpecify the geometryType Flow_Properties Includes variables to specify the flow,Type ControlsIncludes the variables used to control execution Type AtmosphereIncludes all data related to the atmospheric propertiesType Discharge_PropertiesDescribe the properties of the water from a lateral dischargeType DischargesLists the Lateral DischargesType Boundary_ConditionIncludes de data used to specify boundary conditionsType PropertyIncludes all the variables necessary to describe a propertyType PropertyListFormed by the list of properties to be simulatedType Equation_CoefIncludes the coefficients used to solve the evolution equation.River geometryA river is a water body where the length is much longer than the width and the depth ( REF _Ref39163635 \h Figure 1). REF _Ref39163833 \h Figure 3 shows the parameters necessary to describe the geometry of a river trench: width, depth, length and the flow velocity. In a real river the width and the average cross-section depth vary along the length. Along a cross-section, the depth and the width also vary in time. According to the river flow. The most intense the flow is, deeper and wider the cross-section becomes. In this course the aim is to introduce the transport processes and we are going to consider simple river with uniform width, depth, and flow. Considering average values, the river geometry is characterized by three values (length, width and depth) and the flow is characterized by the water velocity. To these values one must add the longitudinal diffusivity and the values necessary to describe the numerical grid: spatial step and number of cells and derived from these values one can compute the area of the cross-section and the horizontal area and volume of one cell. Figure SEQ Figure \* ARABIC 1: Image of a river, a water body where length is much longer than depth or width. Therefore, properties globally change much more along the length than across the vertical cross section and thus rivers are one-dimensional systems.Figure SEQ Figure \* ARABIC 2: Schematic representation of a river geometry (adapted from )In real rivers, another important variable is the distance from the center of a cell to a reference in the river (as we are used to see in motorways). REF _Ref39163833 \h Figure 3 shows the list of variables necessary to define the grid and REF _Ref39165439 \h Figure 4 shows the list of variables necessary to characterize the flow. The list also includes the river slope that is necessary to characterize turbulence. All these parameters are entered in the worksheet river.Figure SEQ Figure \* ARABIC 3: List of variables necessary to define the grid.Figure SEQ Figure \* ARABIC 4: List of variables necessary to define the flow. REF _Ref39166268 \h Figure 5 shows one example of river input data. Actual data is entered in column 2 in front of the respective keyword specified in column 1. Figure SEQ Figure \* ARABIC 5: Data provided in worksheet “River” to specify the river Geometry and Flow.Variables necessary to control the program executionThe simulation can use various methods to discretize the spatial and temporal derivatives and different time steps. These variables as well as all the variables necessary to control printing are grouped in the data structure “Controls” listed in REF _Ref39167374 \h Figure 6. REF _Ref39167376 \h Figure 7 shows the corresponding data provided in the worksheet “Main”. Data is read in column 2, in the line identified by the corresponding keyword. The Boolean variables “advection” and “diffusion” specify if the simulation is going to compute these processes. These keywords are in fact redundant because they can be set off specifying respectively “zero velocity” and/or “zero diffusivity” in the “River” worksheet.Figure SEQ Figure \* ARABIC 6: List of variables to control execution.Figure SEQ Figure \* ARABIC 7: Example of Input data to control execution, provided in the worksheet “Main”.Discharges REF _Ref39182576 \h Figure 11 shows a typical anthropogenic discharge in a river. Natural discharges are tributaries of a main river. Being natural or manmade, a discharge is characterized by its location (and eventually by a name) and by a volume discharge rate (m3/s) and a set of specific properties (e.g. temperature, colour, and a set of concentrations) measured in the framework of a monitoring programs or computed by a model. In real cases, the discharge flow rate and the concentrations depend on time. In this application they will be set constant. In the “Discharges Worksheet” the values of every property simulated by this application must be provided ( to be used as a boundary condition). Figure SEQ Figure \* ARABIC 8: Example of a (lateral) discharge in a River.Figure SEQ Figure \* ARABIC 9: Structure to store the list of discharges. REF _Ref39183486 \h Figure 10 shows the variables need to characterize a particular property discharge and REF _Ref39182576 \h Figure 11 show the data provided in the Worksheet Discharges for all discharges and properties.Figure SEQ Figure \* ARABIC 10: List of variables that characterise a discharge of a specific property.Figure SEQ Figure \* ARABIC 11: Example of data to characterise all discharges and all properties, in the worksheet “Discharges”.Discharges are grouped into a “DischargesList” as shown in REF _Ref39184318 \h Figure 9. The list has only 2 variables, the number of discharges and a vector of “Discharges”, each of them having the properties shown in REF _Ref39183486 \h Figure 10.The discharges will be used as boundary conditions for each property. The calculation of a property is not tied to the calculation of other properties through the discharges. Properties can be tied through source/sink terms only. For that reason, the discharges are handled independently for each variable. Therefore, in the definition of the Discharges there is no reference to the number of properties. The rationale for the algorithm will become fully clear below when a property will be described.Definition of a propertyA property is characterized by the list of variables shown on REF _Ref39485617 \h Figure 12. This list of variables includes variables to quantify decaying/growth of the property. The Boundary conditions associated to the property are defined using the structure “BoundaryCondition” listed in REF _Ref39486697 \h Figure 13. This structure includes: (a) the open boundaries at the tops of the channel, (b) the free surface and (c) the list of lateral discharges described in REF _Ref39184318 \h Figure 9. Each discharge has the properties listed in REF _Ref39183486 \h Figure 10. Using this methodology to describe the discharges, each property gets from the worksheet “Discharges” only the information it needs, location, flowrate and the corresponding specify value (concentration of temperature). This allows the discharges of one property to be imposed independently of other properties. Figure SEQ Figure \* ARABIC 12: List of properties used to characterize a property.Figure SEQ Figure \* ARABIC 13: List of variables used to specify the boundary conditions at the left and right sides of channel, at the free surface and lateral discharges.The properties could be used independently in the application or grouped in a properties’ list. In this program a list ( REF _Ref39487636 \h Figure 14) is used although the advantage is marginal. To create a new property, it is enough to add it to the list and to add the subroutines describing its source/sink terms. Figure SEQ Figure \* ARABIC 14: List of properties used in the application.Other variablesThe application uses internal variables that the user does not need to know. The coefficients of the algebraic equation to solve are the most important set. Typically, to solve a transport equation, the 3 coefficients that multiply the property values at the calculation point and at the two neighbors plus one independent term (Left, Center, Right and Ti coefficients) are enough. In this case we have 5 multiplying coefficients (plus the independent coefficient) to allow the use of the QUICK algorithm that considers 2 left and 2 right neighbors for every point. In the application we have two sets of 5 coefficients. This speeds the calculation of several properties. In fact, advection and diffusion coefficients are common to all properties (they depend only on velocity, diffusivity and on spatial and time steps, i.e. on the Courant and Diffusion Numbers and not on the property itself) and thus they are common to all properties. We need two sets of coefficients because implicit and semi-implicit methods require the inversion of the matrix and consequently the modification of the coefficients used to define the matrix. For that reason, the subroutine “CoefTransport” computes the coefficients a, b, … g. The other set of coefficients L_Left, Left,… are made equal to the former set of coefficients prior to the calculation of a property. Figure SEQ Figure \* ARABIC 15: List of coefficients used to solve the transport equations.AtmosphereThe Type Atmosphere groups all the variables necessary to specify fluxes across the free surface: Heat, O2 and CO2. For that it needs air temperature, wind speed, air moisture and Solar Radiation. Solar Radiation need the Solar Constant, latitude and the position of the Sun in the sky. To compute the position of the Sun, a simple algorithm is used assuming that it describes a sinus between sun rise and sun set hours. The seasonal cycle is not included in the application. The exchange of sensible, latent and radiation is described in the atmosphere. The exchange coefficients (convection coefficients) are computed as functions of the windspeed.The exchange of gas is computed using the Henry Constants and the Partial Pressure in the atmosphere, as described in the theorical documents. Figure SEQ Figure \* ARABIC 16: List of variables used to describe the Atmosphere.Code StructureThe code is composed by three major parts (as in the case of the fin): initialization, the calculation cycle and graphing. InitializationInitialization includes (a) the declaration of the variables using the structures listed above, (b) the initialization of those variables reading data from the worksheets and (c) printing of the initial values (at time t=0). Properties are initialized by the Function “CreateProperty” since most of the steps are the same. Differences between properties and the name of the worksheet used to initialize a property are controlled by the “Name” of the variable passed as first argument when the Function is invoke. The Function “CreateProperty” is invoked for every property. The first step inside the function is to verify if the property is to be active. If it is not active, the Function is left without further reading. The property CO2 included in the PropertiesList is not yet implemented in the application. The first step to include this property is to initialize it using the CreateProperty. Other steps will include the source term resulting from the BOD respiration (stoichiometric oxidation equation) and then the activation of discharges, transport and the calculation of fluxes across the free surface.The initial values of the properties are specified as in case of the Fin Application. Three possibilities are included: impose the default value everywhere (Keyword DEFAULT), imposing zero everywhere but in the cells specified (CELL method) or imposing the default value every where but in the cells specified (DEFCELL method). REF _Ref39569009 \h Figure 18 shows an example of input file for property Temperature in the worksheet “Temperature”.Figure SEQ Figure \* ARABIC 17: Activities used to initialise the model. The program defines the variables to be used, initialises them using the functions (Setcontrols, creategrid…) . The properties included in the PropertiesList to be used in the program are initialised using the function createProperty. After creating every property, the code tests the consistency of the calculation verifying if the properties necessary are active.Figure SEQ Figure \* ARABIC 18: Example of input data to initialize temperature. In this case the property is not active and thus nothing is read. If the property was active (yes in front of the keyword “PropertyActive:” then, the temperature would be initiated with the DEFAULT values (15?C), the water would enter in the channel at 20 ?C if the velocity is positive or at 10 ?C if the velocity was negative (flow from the right).The Computation Cycle The calculation cycle includes:Discharge Application,Advective/diffusive Transport,Decaying,Free Surface Fluxes. There is no best sequence for the calculation. However, discharges are usually the major disturbing pressure in rivers and for that reason they are calculated first in the application. Transport is common to all properties as described above. The coefficients are computed in the subroutine “CoefTransport” and are subsequently used for each property. The method for advection is selected in the subroutine Transport (explicit/implicit/semi-implicit, central-Differences/upwind/QUICK). This subroutine is called sequentially for every active property. After the transport, the fluxes across the free surface (heat and gas) are computed as well as decaying of the Generic Property and of OD/BOD. At the end of each iteration the application tests if it is time for a new output. If yes, it prints every property in its worksheet and updates the printing control variables. After completing the computing cycle, the application invokes the subroutine “Graphs” and creates all the graphs. REF _Ref39569693 \h Figure 19 shows part of the computational cycle: time increment, discharges, transport, and free surface fluxes, as described in the Figure’s legend. After the steps shown in the Figure, the source and sink terms are invoked (see the application code for details).Figure SEQ Figure \* ARABIC 19: Main part of the Computation cycle. The cycle is repeated while the actual time is lower than the duration of the simulation specified in the Setcontrols Function and stored in the ctrl variable. At each time step the actual time is incremented by “dt”, then the permanent discharges are applied for every active variable, then the coefficients of the transport equation are computed (CoefTransport) then the new value of each property is updated in the subroutine “Transport”using those coefficients. After these steps common to every variable, processes specific of each variable are computed: fluxes across the free surface and source/sink terms (not shown in this list). Figure SEQ Figure \* ARABIC 20: Subroutine Transport. This subroutine manages the way how the coefficients are used to compute transport. It distinguish QUICK form other methods, because this method deals with a pentadiagonal matrix and then it distinguish between explicit, implicit and semi-implicit. Free Surface fluxesTemperature and Dissolved Oxygen have free surface fluxes. REF _Ref39570960 \h Figure 21 shows the subroutine managing heat fluxes across the free surface. This subroutine uses services of three Functions to compute convective heat exchange (sensible and latent heat) and to compute radiation. Oxygen fluxes across the free surface are also computed by a subroutine that uses the services of simple Functions. Figure SEQ Figure \* ARABIC 21: Subroutine managing heat fluxes across the free surface. This subroutine need the convection coefficient (computed by the Function “h”, the radiative flux computed by the Function “FluxRadiative” an the Solar Radiation Flux computed by the Function Solar_Radiation. With those fluxes it updates Temperature due to free surface heat fluxes.In fact, fluxes across the free surface depend on the detail of turbulence in the atmosphere and in the water column. The formulation used in this application is realistic, but simple. Many formulations can be found in the literature, in general developed for large water masses, as lakes and oceans.Decay calculationDecay depends on a “rate of decay”. The calculation is very simple if the rate is constant. This is the case of the “Generic” property in our application, but not the case of BOD. Biological Oxygen demand depend on the activity of organisms and this depends on the Temperature and on free Oxygen availability. When there is no free Oxygen only anaerobic respiration is possible and this is very slow. Consequently, without free Oxygen BOD decaying rate is very low.The rate of decay of the Generic Property can be computed implicitly or explicitly, both methods being equally simple. BOD decay is computed identically, but before the rate of decay has to be calculated. REF _Ref39575271 \h Figure 22 shows the subroutine used to compute BOD (and DO) decay. The rate of decay is computed as a function of the Temperature, Dissolve Oxygen and of a semi-saturation constant following a Michaelis-Menten approach. Figure SEQ Figure \* ARABIC 22: Subroutine to compute DOD decay. It gets the decayrate from the Function kd_BOD as function of the Dissolved oxygen and of the Temperature and then computed the BOD decay. If decay is larger than existing BOD or than existing DO then updates the value. This procedure avoids negative concentrations of BOD and of DO. Steps to include a new Properties in the codeInclusion of new properties in the code is a quite simple process because all properties share the same grid and the transport processes (when water flows, it carries all the dissolved and suspended materials). Properties differ on the boundary conditions – including the free surface fluxes – and on the source and sink terms. Modifications to include a mixt QUICK/upwind methodThe QUICK method does not fully respect the Transportive property of Advection. For that reason, in cells of very high downstream gradient it can generate negative concentrations. A way to remove this problem is to use upwind to compute those cells. To implement this method, one must create a Boolean variable to decide it is to be used and another variable to choose the cells where gradient is considered too big. REF _Ref39655447 \h Figure 24 shows the new Control Type and an example of the input data for the new method called QUICK_UP. These modifications have been included in the application version “v4”.Figure SEQ Figure \* ARABIC 23: List of variables to control execution, including two variables for QUICK_UP control.Figure SEQ Figure \* ARABIC 24: Example of Input data to control execution, provided in the worksheet “Main” for QUICK_UP control .Figure SEQ Figure \* ARABIC 25: New “upwind” coefficients have been added to the Equation_Coef Type, necessary to implement the QUICK_UP method.If the method QUICK_UP is selected, then the Boolean variable QUICK_UP becomes true and the code looks for the keyword “QUICK_UP_Ratio:” and reads the variable QUICK_UP_Ratio in column 2. REF _Ref39655994 \h Figure 26 shows the code added to the Function Set controls to read the data necessary to implement the method QUICK_UP. The implementation of the method requires the definition of new coefficient computed as in upwind (b_up, e_up and f_up in REF _Ref39656395 \h Figure 25). These coefficient are computed in the subroutine CoefTransport if QUIK_UP is true ( REF _Ref39656559 \h Figure 27) and subsequently used where property gradients are very high. REF _Ref39656763 \h Figure 28 shows the code added to the subroutines dealing with the QUICK transport calculation. Figure SEQ Figure \* ARABIC 26: Code added to the Function Setcontrols to read the data necessary to implement the QUICK_UP method..Figure SEQ Figure \* ARABIC 27: Code added to the subroutine CoefTransport to implement the QUICK_UP method.Figure SEQ Figure \* ARABIC 28: Code added to the subroutine SemiImpQUICKTransport, ImpQUICKTransport and ExpQUICKTransport to implement the QUICK_UP method.Closing notesThis text aims to give an overview of the data structures and of the algorithm to simulate the advective and diffusive transport of Properties, with or without decay. This overview is useful to understand the calculation algorithm and is necessary to modify it due to the inclusion of alternative processes or due to the inclusion of new properties.The code of this application can be merged with the code of the “Fin Application” to create a “quasi-2D) model able to simulate a heating solar panel. In that case, the water pipe receiving the heat would be “the river” and the heat received by the river along the free surface would be replaced by heat delivered by the Fin to the water. The processes along the Fin are those simulated by the Fin Application. The “base temperature” for the fin would be the water temperature in the pipe. For each cell in the pipe, one would have a fin. ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download