199 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
201 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
202 using GridView = GetPropType<TypeTag, Properties::GridView>;
205 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
206 using Indices = GetPropType<TypeTag, Properties::Indices>;
208 numPhases = FluidSystem::numPhases,
211 temperatureIdx = Indices::temperatureIdx,
212 energyEqIdx = Indices::energyEqIdx,
215 H2OIdx = FluidSystem::H2OIdx,
216 AirIdx = FluidSystem::AirIdx,
219 liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
220 gasPhaseIdx = FluidSystem::gasPhaseIdx,
223 conti0EqIdx = Indices::conti0EqIdx,
226 dim = GridView::dimension,
227 dimWorld = GridView::dimensionworld
230 static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
232 using EqVector = GetPropType<TypeTag, Properties::EqVector>;
233 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
234 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
235 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
236 using Constraints = GetPropType<TypeTag, Properties::Constraints>;
237 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
238 using Model = GetPropType<TypeTag, Properties::Model>;
239 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
240 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
241 using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
242 using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
244 using CoordScalar =
typename GridView::ctype;
245 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
247 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
254 : ParentType(simulator)
262 ParentType::finishInit();
267 FluidSystem::init(275, 600, 100,
273 fineK_ = this->toDimMatrix_(1e-13);
274 coarseK_ = this->toDimMatrix_(1e-12);
278 coarsePorosity_ = 0.3;
281 fineMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
282 fineMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
283 coarseMaterialParams_.setResidualSaturation(liquidPhaseIdx, 0.2);
284 coarseMaterialParams_.setResidualSaturation(gasPhaseIdx, 0.0);
287 fineMaterialParams_.setEntryPressure(1e4);
288 coarseMaterialParams_.setEntryPressure(1e4);
289 fineMaterialParams_.setLambda(2.0);
290 coarseMaterialParams_.setLambda(2.0);
292 fineMaterialParams_.finalize();
293 coarseMaterialParams_.finalize();
296 computeThermalCondParams_(fineThermalCondParams_, finePorosity_);
297 computeThermalCondParams_(coarseThermalCondParams_, coarsePorosity_);
300 solidEnergyLawParams_.setSolidHeatCapacity(790.0
302 solidEnergyLawParams_.finalize();
315 std::ostringstream oss;
316 oss <<
"waterair_" << Model::name();
317 if (getPropValue<TypeTag, Properties::EnableEnergy>())
335 this->model().globalStorage(storage);
338 if (this->gridView().comm().rank() == 0) {
339 std::cout <<
"Storage: " << storage << std::endl << std::flush;
350 template <
class Context>
353 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
354 if (isFineMaterial_(pos))
362 template <
class Context>
363 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
365 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
366 if (isFineMaterial_(pos))
367 return finePorosity_;
369 return coarsePorosity_;
375 template <
class Context>
378 unsigned timeIdx)
const
380 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
381 if (isFineMaterial_(pos))
382 return fineMaterialParams_;
384 return coarseMaterialParams_;
392 template <
class Context>
393 const SolidEnergyLawParams&
397 {
return solidEnergyLawParams_; }
402 template <
class Context>
403 const ThermalConductionLawParams&
406 unsigned timeIdx)
const
408 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
409 if (isFineMaterial_(pos))
410 return fineThermalCondParams_;
411 return coarseThermalCondParams_;
429 template <
class Context>
431 const Context& context,
432 unsigned spaceIdx,
unsigned timeIdx)
const
434 const auto& pos = context.cvCenter(spaceIdx, timeIdx);
435 assert(onLeftBoundary_(pos) ||
436 onLowerBoundary_(pos) ||
437 onRightBoundary_(pos) ||
438 onUpperBoundary_(pos));
441 RateVector massRate(0.0);
442 massRate[conti0EqIdx + AirIdx] = -1e-3;
445 values.setMassRate(massRate);
448 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
449 initialFluidState_(fs, context, spaceIdx, timeIdx);
451 Scalar hl = fs.enthalpy(liquidPhaseIdx);
452 Scalar hg = fs.enthalpy(gasPhaseIdx);
453 values.setEnthalpyRate(values[conti0EqIdx + AirIdx] * hg +
454 values[conti0EqIdx + H2OIdx] * hl);
457 else if (onLeftBoundary_(pos) || onRightBoundary_(pos)) {
458 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
459 initialFluidState_(fs, context, spaceIdx, timeIdx);
462 values.setFreeFlow(context, spaceIdx, timeIdx, fs);
482 template <
class Context>
484 const Context& context,
486 unsigned timeIdx)
const
488 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
489 initialFluidState_(fs, context, spaceIdx, timeIdx);
492 values.assignMassConservative(fs, matParams,
true);
501 template <
class Context>
511 bool onLeftBoundary_(
const GlobalPosition& pos)
const
512 {
return pos[0] < eps_; }
514 bool onRightBoundary_(
const GlobalPosition& pos)
const
515 {
return pos[0] > this->boundingBoxMax()[0] - eps_; }
517 bool onLowerBoundary_(
const GlobalPosition& pos)
const
518 {
return pos[1] < eps_; }
520 bool onUpperBoundary_(
const GlobalPosition& pos)
const
521 {
return pos[1] > this->boundingBoxMax()[1] - eps_; }
523 bool onInlet_(
const GlobalPosition& pos)
const
524 {
return onLowerBoundary_(pos) && (15.0 < pos[0]) && (pos[0] < 25.0); }
526 bool inHighTemperatureRegion_(
const GlobalPosition& pos)
const
527 {
return (20 < pos[0]) && (pos[0] < 30) && (pos[1] < 30); }
529 template <
class Context,
class Flu
idState>
530 void initialFluidState_(FluidState& fs,
531 const Context& context,
533 unsigned timeIdx)
const
535 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
537 Scalar densityW = 1000.0;
538 fs.setPressure(liquidPhaseIdx, 1e5 + (maxDepth_ - pos[1])*densityW*9.81);
539 fs.setSaturation(liquidPhaseIdx, 1.0);
540 fs.setMoleFraction(liquidPhaseIdx, H2OIdx, 1.0);
541 fs.setMoleFraction(liquidPhaseIdx, AirIdx, 0.0);
543 if (inHighTemperatureRegion_(pos))
544 fs.setTemperature(380);
546 fs.setTemperature(283.0 + (maxDepth_ - pos[1])*0.03);
549 fs.setSaturation(gasPhaseIdx, 0);
550 Scalar pc[numPhases];
552 MaterialLaw::capillaryPressures(pc, matParams, fs);
553 fs.setPressure(gasPhaseIdx, fs.pressure(liquidPhaseIdx) + (pc[gasPhaseIdx] - pc[liquidPhaseIdx]));
555 typename FluidSystem::template ParameterCache<Scalar> paramCache;
556 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
557 CFRP::solve(fs, paramCache, liquidPhaseIdx,
true,
true);
560 void computeThermalCondParams_(ThermalConductionLawParams& params, Scalar poro)
562 Scalar lambdaGranite = 2.8;
565 Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
566 fs.setTemperature(293.15);
567 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
568 fs.setPressure(phaseIdx, 1.0135e5);
571 typename FluidSystem::template ParameterCache<Scalar> paramCache;
572 paramCache.updateAll(fs);
573 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
574 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
575 fs.setDensity(phaseIdx, rho);
578 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
579 Scalar lambdaSaturated;
580 if (FluidSystem::isLiquid(phaseIdx)) {
582 FluidSystem::thermalConductivity(fs, paramCache, phaseIdx);
583 lambdaSaturated = std::pow(lambdaGranite, (1-poro)) + std::pow(lambdaFluid, poro);
586 lambdaSaturated = std::pow(lambdaGranite, (1-poro));
588 params.setFullySaturatedLambda(phaseIdx, lambdaSaturated);
589 if (!FluidSystem::isLiquid(phaseIdx))
590 params.setVacuumLambda(lambdaSaturated);
594 bool isFineMaterial_(
const GlobalPosition& pos)
const
595 {
return pos[dim-1] > layerBottom_; }
601 Scalar finePorosity_;
602 Scalar coarsePorosity_;
604 MaterialLawParams fineMaterialParams_;
605 MaterialLawParams coarseMaterialParams_;
607 ThermalConductionLawParams fineThermalCondParams_;
608 ThermalConductionLawParams coarseThermalCondParams_;
609 SolidEnergyLawParams solidEnergyLawParams_;