Buildings.Fluid.Geothermal.Aquifer

Package with aquifer thermal energy storage models

Information

Package with models for aquifer thermal energy storage.

Extends from Modelica.Icons.Package (Icon for standard packages).

Package Content

Name Description
Buildings.Fluid.Geothermal.Aquifer.MultiWell MultiWell Model of a single well for aquifer thermal energy storage
Buildings.Fluid.Geothermal.Aquifer.Data Data Collection of data records for aquifer thermal energy storage
Buildings.Fluid.Geothermal.Aquifer.Examples Examples Example models for Buildings.Fluid.Geothermal.Aquifer
Buildings.Fluid.Geothermal.Aquifer.Validation Validation Validation models for aquifer thermal energy storage

Buildings.Fluid.Geothermal.Aquifer.MultiWell Buildings.Fluid.Geothermal.Aquifer.MultiWell

Model of a single well for aquifer thermal energy storage

Buildings.Fluid.Geothermal.Aquifer.MultiWell

Information

This model simulates aquifer thermal energy storage, using one or multiple pairs of cold and hot wells.

To calculate aquifer temperature at different locations over time, the model applies physical principles of water flow and heat transfer phenomena. The model is based on the partial differential equation (PDE) for 1D conductive-convective transient radial heat transport in porous media

ρ c (∂ T(r,t) ⁄ ∂t) = k (∂² T(r,t) ⁄ ∂r²) - ρw cw u(∂ T(r,t) ⁄ ∂t),

where ρ is the mass density, c is the specific heat capacity per unit mass, T is the temperature at location r and time t, u is water velocity and k is the heat conductivity. The subscript w indicates water. The first term on the right hand side of the equation describes the effect of conduction, while the second term describes the fluid flow.

The pressure losses in the aquifer are calculated using the Darcy's law

Δp = ṁ g ⁄ (2 π K h ln(rMax ⁄ rWB)),

where is the water mass flow rate, g is the gravitational acceleration, K is the hydraulic conductivity, h is the thickness of the aquifer, rMax is the domain radius and rWB is the well radius. The pressure losses in the wells are calculated using Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow.

Spatial discretization

To discretize the conductive-convective equation, the domain is divided into a series of thermal capacitances and thermal resistances along the radial direction. The implementation uses an array of Modelica.Thermal.HeatTransfer.Components.HeatCapacitor and Modelica.Thermal.HeatTransfer.Components.ThermalResistor. Fluid flow was modelled by adding a series of fluid volumes, which are connected to the thermal capacitances via heat ports. The fluid stream was developed using the model Buildings.Fluid.MixingVolumes.MixingVolume. The geometric representation of the model is illustrated in the figure below.

image

Typical use and important parameters

By default, the component consists of a single pair of wells: one cold well and one warm well. The number of paired wells can be increased by modifing the parameters nPai. The effect is a proportional increase of thermal capacity, and the mass flow rate at port_a and port_b is equally distributed to the pairs of well, thus all pairs have the same mass flow rates and temperatures, and the quantities at the fluid ports is for all well combined, as is the electricity consumption PTot for the well pumps.

To ensure conservation of energy, the wells are connected via fluid ports. To avoid thermal interferences, make sure that the aquifer domain radius rMax is large enough for your specific use case.

Circulation pumps are included in the model and they can be controlled by acting on the input connector. The input must vary between [1, -1]. A positive value will circulate water clockwise (from port_Hot to port_Col, thus extraction from the cold well and injection into the warm well). A negative value will circulate water anticlockwise.

The temperature values in the warm and cold aquifers can be accessed using TAquHot and TAquCol. These temperatures correspond to the temperatures of each thermal capacitance in the discretized domain. The location of the thermal capacitance is expressed by rVol.

Implementation

The model computes the flow resistance of the cold well and the cold acquifer. Because of symmetry, the warm side has the same flow resistances. To reduce the number of equations in the model, the flow resistance of the warm side is set to the same flow resistance as on the cold side, but with the sign reversed.

Parameters

TypeNameDefaultDescription
replaceable package MediumModelica.Media.Interfaces.Pa...Medium in the component
Heighth Aquifer thickness [m]
Lengthd Distance between wells [m]
Heightlength Length of one well, used to compute pressure drop [m]
RealnPai1Number of paired wells
TemplateaquDat Aquifer thermal properties
Subsurface
IntegernVol10Number of control volumes used in discretization
RadiusrWB0.1Wellbore radius [m]
RadiusrMaxd/2Domain radius [m]
RealgriFac1.15Grid factor for spacing
TemperatureTCol_start283.15Initial temperature of cold well [K]
TemperatureTHot_startTCol_startInitial temperature of warm well [K]
Properties of ground
TemperatureTGroCol285.15Undisturbed ground temperature (cold well) [K]
TemperatureTGroHotTGroColUndisturbed ground temperature (warm well) [K]
Hydraulic circuit
MassFlowRatem_flow_nominal Nominal mass flow rate [kg/s]
PressureDifferencedpAquifer_nominalm_flow_nominal/nPai*Modelica...Pressure drop at nominal mass flow rate in the aquifer [Pa]
PressureDifferencedpExt_nominal Pressure drop at nominal mass flow rate in the above-surface system (used to size the head of the well pump) [Pa]

Connectors

TypeNameDescription
replaceable package MediumMedium in the component
input RealInputuPump control input (-1: extract from hot well, +1: extract from cold well, 0: off) [1]
output RealOutputPTotTotal power consumed by all circulation pumps in the wells [W]
FluidPort_aport_ColFluid connector
FluidPort_aport_HotFluid connector

Modelica definition

model MultiWell "Model of a single well for aquifer thermal energy storage" replaceable package Medium = Modelica.Media.Interfaces.PartialMedium "Medium in the component"; parameter Integer nVol(min=1)=10 "Number of control volumes used in discretization"; parameter Modelica.Units.SI.Height h "Aquifer thickness"; parameter Modelica.Units.SI.Length d "Distance between wells"; parameter Modelica.Units.SI.Height length "Length of one well, used to compute pressure drop"; parameter Real nPai=1 "Number of paired wells"; parameter Modelica.Units.SI.Radius rWB=0.1 "Wellbore radius"; parameter Modelica.Units.SI.Radius rMax=d/2 "Domain radius"; parameter Real griFac(min=1) = 1.15 "Grid factor for spacing"; parameter Modelica.Units.SI.Temperature TCol_start=283.15 "Initial temperature of cold well"; parameter Modelica.Units.SI.Temperature THot_start=TCol_start "Initial temperature of warm well"; parameter Modelica.Units.SI.Temperature TGroCol=285.15 "Undisturbed ground temperature (cold well)"; parameter Modelica.Units.SI.Temperature TGroHot=TGroCol "Undisturbed ground temperature (warm well)"; parameter Buildings.Fluid.Geothermal.Aquifer.Data.Template aquDat "Aquifer thermal properties"; parameter Modelica.Units.SI.MassFlowRate m_flow_nominal "Nominal mass flow rate"; parameter Modelica.Units.SI.PressureDifference dpAquifer_nominal(displayUnit= "Pa")= m_flow_nominal/nPai*Modelica.Constants.g_n/2/Modelica.Constants.pi/h/aquDat.K*log(d/rWB) "Pressure drop at nominal mass flow rate in the aquifer"; final parameter Modelica.Units.SI.PressureDifference dpWell_nominal( displayUnit="Pa") = 2*resWelCol.dp_nominal "Pressure drop at nominal mass flow rate in the well"; parameter Modelica.Units.SI.PressureDifference dpExt_nominal(displayUnit="Pa") "Pressure drop at nominal mass flow rate in the above-surface system (used to size the head of the well pump)"; final parameter Modelica.Units.SI.Radius rVol[nVol]( each final fixed=false) "Radius to the center of the i-th domain"; Modelica.Blocks.Interfaces.RealInput u( final min=-1, final max=1, final unit="1") "Pump control input (-1: extract from hot well, +1: extract from cold well, 0: off)"; Modelica.Blocks.Interfaces.RealOutput PTot( final unit="W") "Total power consumed by all circulation pumps in the wells"; Modelica.Fluid.Interfaces.FluidPort_a port_Col( redeclare final package Medium = Medium) "Fluid connector"; Modelica.Fluid.Interfaces.FluidPort_a port_Hot( redeclare final package Medium = Medium) "Fluid connector"; Modelica.Units.SI.Temperature TAquHot[nVol] = heaCapHot.T "Temperatures of the hot aquifer"; Modelica.Units.SI.Temperature TAquCol[nVol] = heaCapCol.T "Temperatures of the cold aquifer"; Movers.Preconfigured.SpeedControlled_y pumCol( redeclare final package Medium = Medium, final m_flow_nominal=m_flow_nominal/nPai, final dp_nominal=resAqu.dpMea_nominal + 2* resWelCol.dp_nominal + dpExt_nominal) "Pump to extract from cold well"; Movers.Preconfigured.SpeedControlled_y pumHot( redeclare final package Medium = Medium, final m_flow_nominal=pumCol.m_flow_nominal, final dp_nominal=pumCol.dp_nominal) "Pump to extract from hot well"; FixedResistances.PressureDrop resWelCol( redeclare final package Medium = Medium, m_flow_nominal=m_flow_nominal/nPai, final from_dp=false, dp_nominal= Modelica.Fluid.Pipes.BaseClasses.WallFriction.Detailed.pressureLoss_m_flow( m_flow=m_flow_nominal/nPai, rho_a=rhoWat, rho_b=rhoWat, mu_a=mu, mu_b=mu, length=length, diameter=rWB, roughness=2.5e-5)) "Pressure drop in the cold well"; // The resistances on the hot side are set to be the negative of the cold // side, as they have the same geometry and flow rate. // This reduces the size of the nonlinear system of equations, // and leads to fewer equations in the model. Movers.BaseClasses.IdealSource resWelHot( redeclare final package Medium = Medium, final allowFlowReversal=true, final control_m_flow=false, final control_dp=true, final m_flow_small=1E-4*abs(m_flow_nominal)) "Pressure drop in the warm well"; BaseClasses.MassFlowRateMultiplier mulCol( redeclare final package Medium = Medium, k=nPai) "Mass flow multiplier for the cold well"; BaseClasses.MassFlowRateMultiplier mulHot( redeclare final package Medium = Medium, k=nPai) "Mass flow multiplier for the warm well"; Modelica.Blocks.Math.Add addPum( final k1=1, final k2=1) "Sum of pump electrical power"; Sensors.RelativePressure dpWelCol( redeclare package Medium = Medium) "Pressure drop of cold well"; Airflow.Multizone.Point_m_flow resAqu( redeclare final package Medium = Medium, m=1, final dpMea_nominal=dpAquifer_nominal, final mMea_flow_nominal=m_flow_nominal/nPai) "Pressure drop in the cold and warm side of the aquifer"; protected parameter Modelica.Units.SI.Radius r[nVol + 1](each fixed=false) "Radius to the boundary of the i-th domain"; parameter Modelica.Units.SI.HeatCapacity C[nVol](each fixed=false) "Heat capacity of segment"; parameter Modelica.Units.SI.Volume VWat[nVol](each fixed=false) "Volumes of water"; parameter Modelica.Units.SI.ThermalResistance R[nVol](each fixed=false) "Thermal resistances between nodes"; parameter Real cAqu(fixed=false) "Heat capacity normalized with volume for aquifer"; parameter Real kVol(fixed=false) "Heat conductivity normalized with volume"; parameter Modelica.Units.SI.DynamicViscosity mu=Medium.dynamicViscosity(Medium.setState_pTX( Medium.p_default, Medium.T_default, Medium.X_default)) "Water dynamic viscosity"; parameter Modelica.Units.SI.SpecificHeatCapacity cpWat=Medium.specificHeatCapacityCp(Medium.setState_pTX( Medium.p_default, Medium.T_default, Medium.X_default)) "Water specific heat capacity"; parameter Modelica.Units.SI.Density rhoWat=Medium.density(Medium.setState_pTX( Medium.p_default, Medium.T_default, Medium.X_default)) "Water density"; parameter Modelica.Units.SI.ThermalConductivity kWat=Medium.thermalConductivity(Medium.setState_pTX( Medium.p_default, Medium.T_default, Medium.X_default)) "Water thermal conductivity"; Buildings.Fluid.MixingVolumes.MixingVolume volCol[nVol]( redeclare final package Medium = Medium, each final T_start=TCol_start, each final m_flow_nominal=m_flow_nominal, final V=VWat, each nPorts=2) "Array of fluid volumes representing the fluid flow in the cold side of the aquifer"; Modelica.Thermal.HeatTransfer.Components.HeatCapacitor heaCapCol[nVol](C=C, each T(start=TCol_start, fixed=true)) "Array of thermal capacitor in the cold side of the aquifer"; Modelica.Thermal.HeatTransfer.Components.ThermalResistor theResCol[nVol](R=R) "Array of thermal resistances in the cold side of the aquifer"; MixingVolumes.MixingVolume volHot[nVol]( redeclare final package Medium = Medium, each T_start=THot_start, each m_flow_nominal=m_flow_nominal, V=VWat, each nPorts=2) "Array of fluid volumes representing the fluid flow in the warm side of the aquifer"; Modelica.Thermal.HeatTransfer.Components.HeatCapacitor heaCapHot[nVol](C=C, each T(start=THot_start, fixed=true)) "Array of thermal capacitor in the warm side of the aquifer"; Modelica.Thermal.HeatTransfer.Components.ThermalResistor theResHot[nVol](R=R) "Array of thermal resistances in the warm side of the aquifer"; Buildings.HeatTransfer.Sources.FixedTemperature groTemCol( final T=TGroCol) "Boundary condition ground temperature in the cold side of the aquifer"; Buildings.HeatTransfer.Sources.FixedTemperature groTemHot( final T=TGroHot) "Boundary condition ground temperature in the warm side of the aquifer"; Modelica.Blocks.Nonlinear.Limiter limCol(final uMax=1, final uMin=0) "Limiter for pump signal"; Modelica.Blocks.Math.Gain gaiCon(final k=-1) "Inversion of control signal"; Modelica.Blocks.Nonlinear.Limiter limHot(final uMax=1, final uMin=0) "Limiter for pump signal"; Modelica.Blocks.Math.Gain gaiPum(final k=nPai) "Gain used to scale the pump electrical power"; initial equation assert(rWB < rMax, "In " + getInstanceName() + ": Require rWB < rMax"); assert(0 < rWB, "In " + getInstanceName() + ": Require 0 < rWB"); assert(rMax <= d/2, "In " + getInstanceName() + ": Require rMax <= d/2"); cAqu=aquDat.dSoi*aquDat.cSoi*(1-aquDat.phi); kVol=kWat*aquDat.phi+aquDat.kSoi*(1-aquDat.phi); r[1] = rWB; for i in 2:nVol+1 loop r[i]= r[i-1] + (rMax - rWB) * (1-griFac)/(1-griFac^(nVol)) * griFac^(i-2); end for; for i in 1:nVol loop rVol[i] = (r[i]+r[i+1])/2; end for; for i in 1:nVol loop C[i] = cAqu*h*3.14*(r[i+1]^2-r[i]^2); end for; for i in 1:nVol loop VWat[i] = aquDat.phi*h*3.14*(r[i+1]^2-r[i]^2); end for; R[nVol]=Modelica.Math.log(rMax/rVol[nVol])/(2*Modelica.Constants.pi*kVol*h); for i in 1:nVol-1 loop R[i] = Modelica.Math.log(rVol[i+1]/rVol[i])/(2*Modelica.Constants.pi*kVol*h); end for; equation if nVol > 1 then for i in 1:(nVol - 1) loop connect(volCol[i].ports[2], volCol[i + 1].ports[1]); end for; end if; connect(resAqu.port_b, volCol[nVol].ports[2]); connect(resAqu.port_a, volHot[nVol].ports[2]); connect(pumCol.port_a, volCol[1].ports[1]); connect(pumHot.port_a, volHot[1].ports[1]); if nVol > 1 then for i in 1:(nVol - 1) loop connect(volHot[i].ports[2], volHot[i + 1].ports[1]); end for; end if; if nVol > 1 then for i in 1:(nVol - 1) loop connect(heaCapCol[i + 1].port, theResCol[i].port_b); end for; end if; if nVol > 1 then for i in 1:(nVol - 1) loop connect(heaCapHot[i + 1].port, theResHot[i].port_b); end for; end if; connect(groTemCol.port, theResCol[nVol].port_b); connect(volCol.heatPort, heaCapCol.port); connect(theResCol.port_a, volCol.heatPort); connect(groTemHot.port, theResHot[nVol].port_b); connect(volHot.heatPort, heaCapHot.port); connect(theResHot.port_a, volHot.heatPort); connect(pumCol.port_b, resWelCol.port_a); connect(resWelHot.port_a, pumHot.port_b); connect(limCol.y, pumCol.y); connect(gaiCon.y, limHot.u); connect(limHot.y, pumHot.y); connect(gaiCon.u, u); connect(limCol.u, u); connect(resWelCol.port_b, mulCol.port_a); connect(mulCol.port_b, port_Col); connect(mulHot.port_b, port_Hot); connect(resWelHot.port_b, mulHot.port_a); connect(pumCol.P, addPum.u1); connect(addPum.y, gaiPum.u); connect(pumHot.P, addPum.u2); connect(gaiPum.y, PTot); connect(dpWelCol.port_a, resWelCol.port_b); connect(dpWelCol.port_b, resWelCol.port_a); connect(dpWelCol.p_rel, resWelHot.dp_in); end MultiWell;