Dear Dumuxers,
I'm simulating a single phase geothermal doublet (a production and an
injection well in a geothermal aquifer where a low temperature is
imposed for the injected water = 293 degree kelvin). The boundary
conditions are the same as initial conditions and they are all Dirichlet
with a geothermal gradient (0.03 Kelvin/m) and a hydrostatic pressure
gradient where they have values of 358 Kelvin and 1.1e6 Pa respectively
at the top of the upper boundary layer (top of domain). The domain is
2500 m x 2500 m x 50 m (30 m aquifer and 10 m upper layer and 10 m lower
layer) with 3 refinement zones (coarser outside to refined inside). The
upper and lower layers are represented by 3 numerical layers each to
ensure undisturbed boundaries. I'm representing the aquifer by 2
numerical layers and each of the injection & production wells by 2 point
sources in the 2 layers (fully penetrating wells) with 100 m distance
between the 2 wells and with a rate of 100 m^3/day for each well. I'm
using the OnePNI model and the horizontal permeability of the aquifer is
250 mD and the vertical one is 100 mD. The permeability of the boundary
layers is negligible (1e-9 mD).
When I set the simulation time for 1 month, at around 14 days, the time
step size kept decreasing from 50000 s till it reached 44 s ! so it
became impractical to proceed and I stopped the simulation. When I
displayed the results on Paraview, I could see what looks like a
numerical error above and below and point sources of the injection well
in the most 2 inner layers for pressure and and the most inner layer for
temperature (pictures attached). I tried the decreasing the vertical
permeability of the aquifer to 1e-9 mD but the error persisted. I also
tried deceasing the convergence criterion from 1e-3 to 1e-7 while
keeping the new low vertical permeability but again the error persisted.
Only when I removed the imposed temperature of 293 Kelvin and I used the
current temperature of the control volume (volVars.temperature()) to
decouple the thermal problem from the flow problem, the simulation
proceeded successfully till tEnd = 20 days with easy convergence and the
code seemed really stable !!
I tried today re-imposing the 293 kelvin for the injected water but with
changing the thickness of the most inner boundary layers to 10 m instead
of 3.33 m and the simulation proceeded successfully (tEnd = 20 days) but
I could still see the numerical error. Also I tried representing each of
the boundaries with just 1 numerical layer of 10 m and again the
simulation proceeded successfully (tEnd = 20 days) but I could still
observe the numerical error above and below the injection well inside
the boundaries.
I'm really confused at this point. Why is the imposed temperature
causing this numerical anomaly ? Does the point source nature have
anything to do with that ? Is it because of the heat transfer by
conduction at the boundary coupled with the low permeability? I
appreciate your thoughts and guidance on the issue.
here are my point source functions and my grid
void addPointSources(std::vector<PointSource>& pointSources) const
{
// The injection well (source term)
pointSources.push_back(PointSource({1205, 1245, 17.5}));
pointSources.push_back(PointSource({1205, 1245, 32.5}));
// The production well (sink term)
pointSources.push_back(PointSource({1305, 1245, 17.5}));
pointSources.push_back(PointSource({1305, 1245, 32.5}));
}
// Using solution-dependent point sources
template<class ElementVolumeVariables>
void pointSource(PointSource& source,
const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
const auto& pos = source.position();
const auto& volVars = elemVolVars[scv];
if (pos[0] < 1250.0) // injection
{
const Scalar volumeSource = 5.787e-4; // injectionRate is
positive and m^3/s
const Scalar massSource =
volumeSource*IapwsH2O::liquidDensity(293, volVars.pressure(0));
const Scalar energySource =
massSource*IapwsH2O::liquidEnthalpy(293, volVars.pressure(0));
source = NumEqVector({ massSource, energySource });
}
else // production
{
const Scalar volumeSource = -5.787e-4; // productionRate is
negative and m^3/s
const Scalar massSource = volumeSource*volVars.density(0); //
using current control volume's water density
const Scalar energySource = massSource*volVars.enthalpy(0);
// using current control volume's water enthalpy
source = NumEqVector({ massSource, energySource });
}
}
---------------------------------------------------------------------------------------------
[Grid]
Positions0 = 0 400 1100 1400 2100 2500
Positions1 = 0 400 1100 1400 2100 2500
Positions2 = 0 10 40 50
Cells0 = 2 35 30 35 2
Cells1 = 2 35 30 35 2
Cells2 = 3 2 3
---------------------------------------------------------------------------------------------
Best Regards,
Mahmoud Atef
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
