diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp index 2d78366e2a..ebfd1cde1e 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp @@ -245,152 +245,6 @@ void RigFractureTransCalc::computeTransmissibilityFromPolygonWithInfiniteConduct m_fracture->setFractureData(fracDataVec); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -void RigFractureTransCalc::calculateStimPlanCellsMatrixTransmissibility( RigStimPlanFracTemplateCell* stimPlanCell, RigStimPlanFractureCell* fracStimPlanCellData) -{ - //TODO: Gjør til egen klasse / kalkulator, som kan holde eclipseCell/StimplanCell og Transm. - - //Not calculating flow into fracture if stimPlan cell cond value is 0 (assumed to be outside the fracture): - if (stimPlanCell->getConductivtyValue() < 1e-7) return; - - RigEclipseCaseData* eclipseCaseData = m_case->eclipseCaseData(); - - RifReaderInterface::PorosityModelResultType porosityModel = RifReaderInterface::MATRIX_RESULTS; - RimReservoirCellResultsStorage* gridCellResults = m_case->results(porosityModel); - - size_t scalarSetIndex; - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DX"); - cvf::ref dataAccessObjectDx = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DY"); - cvf::ref dataAccessObjectDy = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "DZ"); - cvf::ref dataAccessObjectDz = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "DZ"); //assuming 0 time step and main grid (so grid index =0) - - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMX"); - cvf::ref dataAccessObjectPermX = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMX"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMY"); - cvf::ref dataAccessObjectPermY = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMY"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "PERMZ"); - cvf::ref dataAccessObjectPermZ = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "PERMZ"); //assuming 0 time step and main grid (so grid index =0) - scalarSetIndex = gridCellResults->findOrLoadScalarResult(RimDefines::STATIC_NATIVE, "NTG"); - cvf::ref dataAccessObjectNTG = RigResultAccessorFactory::createFromUiResultName(eclipseCaseData, 0, porosityModel, 0, "NTG"); //assuming 0 time step and main grid (so grid index =0) - - RigActiveCellInfo* activeCellInfo = eclipseCaseData->activeCellInfo(porosityModel); - - - std::vector stimPlanPolygon = stimPlanCell->getPolygon(); - std::vector fracCells = m_fracture->getPotentiallyFracturedCells(); - - for (size_t fracCell : fracCells) - { - bool cellIsActive = activeCellInfo->isActive(fracCell); - if (!cellIsActive) continue; - - double permX = dataAccessObjectPermX->cellScalarGlobIdx(fracCell); - double permY = dataAccessObjectPermY->cellScalarGlobIdx(fracCell); - double permZ = dataAccessObjectPermZ->cellScalarGlobIdx(fracCell); - - double dx = dataAccessObjectDx->cellScalarGlobIdx(fracCell); - double dy = dataAccessObjectDy->cellScalarGlobIdx(fracCell); - double dz = dataAccessObjectDz->cellScalarGlobIdx(fracCell); - - double NTG = dataAccessObjectNTG->cellScalarGlobIdx(fracCell); - - cvf::Vec3d localX; - cvf::Vec3d localY; - cvf::Vec3d localZ; - std::vector > planeCellPolygons; - bool isPlanIntersected = planeCellIntersectionPolygons(fracCell, planeCellPolygons, localX, localY, localZ); - if (!isPlanIntersected || planeCellPolygons.size() == 0) continue; - - //Transform planCell polygon(s) and averageZdirection to x/y coordinate system (where fracturePolygon already is located) - cvf::Mat4f invertedTransMatrix = m_fracture->transformMatrix().getInverted(); - for (std::vector & planeCellPolygon : planeCellPolygons) - { - for (cvf::Vec3d& v : planeCellPolygon) - { - v.transformPoint(static_cast(invertedTransMatrix)); - } - } - - cvf::Vec3d localZinFracPlane; - localZinFracPlane = localZ; - localZinFracPlane.transformVector(static_cast(invertedTransMatrix)); - cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO; - directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1)); - directionOfLength.normalize(); - - std::vector > polygonsForStimPlanCellInEclipseCell; - cvf::Vec3d areaVector; - - for (std::vector planeCellPolygon : planeCellPolygons) - { - std::vector >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, stimPlanPolygon); - for (std::vector clippedPolygon : clippedPolygons) - { - polygonsForStimPlanCellInEclipseCell.push_back(clippedPolygon); - } - } - - if (polygonsForStimPlanCellInEclipseCell.size() == 0) continue; - - double area; - std::vector areaOfFractureParts; - double length; - std::vector lengthXareaOfFractureParts; - double Ax = 0.0, Ay = 0.0, Az = 0.0; - - for (std::vector fracturePartPolygon : polygonsForStimPlanCellInEclipseCell) - { - areaVector = cvf::GeometryTools::polygonAreaNormal3D(fracturePartPolygon); - area = areaVector.length(); - areaOfFractureParts.push_back(area); - - //TODO: the l in the sl/pi term in the denominator of the Tmj expression should be the length of the full Eclipse cell - //In the current form the implementation gives correct result only if s=0 (fracture templte skin factor). - length = RigCellGeometryTools::polygonAreaWeightedLength(directionOfLength, fracturePartPolygon); - lengthXareaOfFractureParts.push_back(length * area); - - cvf::Plane fracturePlane; - cvf::Mat4f m = m_fracture->transformMatrix(); - bool isCellIntersected = false; - - fracturePlane.setFromPointAndNormal(static_cast(m.translation()), - static_cast(m.col(2))); - - Ax += abs(area*(fracturePlane.normal().dot(localY))); - Ay += abs(area*(fracturePlane.normal().dot(localX))); - Az += abs(area*(fracturePlane.normal().dot(localZ))); - } - - double fractureArea = 0.0; - for (double area : areaOfFractureParts) fractureArea += area; - - double totalAreaXLength = 0.0; - for (double lengtXarea : lengthXareaOfFractureParts) totalAreaXLength += lengtXarea; - - double fractureAreaWeightedlength = totalAreaXLength / fractureArea; - double skinfactor = m_fracture->attachedFractureDefinition()->skinFactor; - - - - double transmissibility_X = calculateMatrixTransmissibility(permY, NTG, Ay, dx, skinfactor, fractureAreaWeightedlength); - double transmissibility_Y = calculateMatrixTransmissibility(permX, NTG, Ax, dy, skinfactor, fractureAreaWeightedlength); - double transmissibility_Z = calculateMatrixTransmissibility(permZ, 1.0, Az, dz, skinfactor, fractureAreaWeightedlength); - - double transmissibility = sqrt(transmissibility_X * transmissibility_X - + transmissibility_Y * transmissibility_Y - + transmissibility_Z * transmissibility_Z); - - - - - fracStimPlanCellData->addContributingEclipseCell(fracCell, transmissibility); - } - -} //-------------------------------------------------------------------------------------------------- /// diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h index 3eea71beb2..efb31f35ca 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h @@ -48,29 +48,39 @@ public: // Calculations based on fracture polygon and eclipse grid cells void computeTransmissibilityFromPolygonWithInfiniteConductivityInFracture(); - bool planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ); // Functions needed for upscaling from StimPlan grid to Eclipse Grid, for transmissibility calculations on eclipse grid // Obsolete if final calculations will be done on the stimPlan grid void computeUpscaledPropertyFromStimPlan(QString resultName, QString resultUnit, size_t timeStepIndex); + + // Calculations based on StimPlan grid + static double computeStimPlanCellTransmissibilityInFracture(double conductivity, double sideLengthParallellTrans, double sideLengthNormalTrans); + + double computeRadialTransmissibilityToWellinStimPlanCell(const RigStimPlanFracTemplateCell& stimPlanCell); + double computeLinearTransmissibilityToWellinStimPlanCell(const RigStimPlanFracTemplateCell& stimPlanCell, + double perforationLengthVertical, + double perforationLengthHorizontal); + + double cDarcy(); + +private: + + bool planeCellIntersectionPolygons(size_t cellindex, + std::vector > & polygons, + cvf::Vec3d & localX, + cvf::Vec3d & localY, + cvf::Vec3d & localZ); + + // Functions needed for upscaling from StimPlan grid to Eclipse Grid, for transmissibility calculations on eclipse grid + // Obsolete if final calculations will be done on the stimPlan grid std::pair flowAcrossLayersUpscaling(QString resultName, QString resultUnit, size_t timeStepIndex, RimDefines::UnitSystem unitSystem, size_t eclipseCellIndex); double computeHAupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers); double computeAHupscale(RimStimPlanFractureTemplate* fracTemplateStimPlan, std::vector stimPlanCells, std::vector planeCellPolygon, cvf::Vec3d directionAlongLayers, cvf::Vec3d directionAcrossLayers); static double arithmeticAverage(std::vector values); - // Calculations based on StimPlan grid - static double computeStimPlanCellTransmissibilityInFracture(double conductivity, double sideLengthParallellTrans, double sideLengthNormalTrans); - - void calculateStimPlanCellsMatrixTransmissibility(RigStimPlanFracTemplateCell* stimPlanCell, RigStimPlanFractureCell* fracStimPlanCellData); - double computeRadialTransmissibilityToWellinStimPlanCell(const RigStimPlanFracTemplateCell& stimPlanCell); - double computeLinearTransmissibilityToWellinStimPlanCell(const RigStimPlanFracTemplateCell& stimPlanCell, double perforationLengthVertical, double perforationLengthHorizontal); - static std::vector getRowOfStimPlanCells(std::vector& allStimPlanCells, size_t i); static std::vector getColOfStimPlanCells(std::vector& allStimPlanCells, size_t j); - double cDarcy(); - -private: double convertConductivtyValue(double Kw, RimDefines::UnitSystem fromUnit, RimDefines::UnitSystem toUnit); double calculateMatrixTransmissibility(double permX, double NTG, double Ay, double dx, double skinfactor, double fractureAreaWeightedlength);