From a1a7441c2d4243e3173886f19128a0ff0602daa3 Mon Sep 17 00:00:00 2001 From: astridkbjorke Date: Tue, 16 May 2017 13:21:39 +0200 Subject: [PATCH] #1487 - pre-proto - Adding class EclipseToStimPlanCellTransmissibilityCalculator --- .../RigFractureTransCalc.cpp | 256 +++++++++++++++++- .../ReservoirDataModel/RigFractureTransCalc.h | 27 +- .../RigStimPlanFracTemplateCell.h | 2 +- 3 files changed, 275 insertions(+), 10 deletions(-) diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp index 451c1cf418..1442bc50a9 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.cpp @@ -41,6 +41,8 @@ #include "RimDefines.h" #include "RigStimPlanFracTemplateCell.h" #include //Used for log +#include "RiaApplication.h" +#include "RimEclipseView.h" //-------------------------------------------------------------------------------------------------- /// @@ -789,6 +791,21 @@ double RigFractureTransCalc::computeLinearTransmissibilityToWellinStimPlanCell(c return Tc; } + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigFractureTransCalc::calculateMatrixTransmissibility(double perm, double NTG, double A, double cellSizeLength, double skinfactor, double fractureAreaWeightedlength) +{ + double transmissibility; + + double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; + transmissibility = 8 * cDarcy() * (perm * NTG) * A / (cellSizeLength + slDivPi); + + return transmissibility; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -859,15 +876,248 @@ double RigFractureTransCalc::convertConductivtyValue(double Kw, RimDefines::Unit return cvf::UNDEFINED_DOUBLE; } + + +EclipseToStimPlanCellTransmissibilityCalculator::EclipseToStimPlanCellTransmissibilityCalculator(const RimEclipseCase* caseToApply, + cvf::Mat4f fractureTransform, + double skinFactor, + double cDarcy, + const RigStimPlanFracTemplateCell& stimPlanCell) + + : m_stimPlanCell(stimPlanCell) +{ + m_case = caseToApply; + m_fractureSkinFactor = skinFactor; + m_fractureTransform = fractureTransform; +} + + + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -double RigFractureTransCalc::calculateMatrixTransmissibility(double perm, double NTG, double A, double cellSizeLength, double skinfactor, double fractureAreaWeightedlength) +void EclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsMatrixTransmissibility() { - double transmissibility; + + //Not calculating flow into fracture if stimPlan cell cond value is 0 (assumed to be outside the fracture): + if (m_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 fracCells = getPotentiallyFracturedCellsForPolygon(m_stimPlanCell.getPolygon()); + 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_fractureTransform.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; + std::vector stimPlanPolygon = m_stimPlanCell.getPolygon(); + + 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; + bool isCellIntersected = false; + + fracturePlane.setFromPointAndNormal(static_cast(m_fractureTransform.translation()), + static_cast(m_fractureTransform.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 transmissibility_X = calculateMatrixTransmissibility(permY, NTG, Ay, dx, m_fractureSkinFactor, fractureAreaWeightedlength); + double transmissibility_Y = calculateMatrixTransmissibility(permX, NTG, Ax, dy, m_fractureSkinFactor, fractureAreaWeightedlength); + double transmissibility_Z = calculateMatrixTransmissibility(permZ, 1.0, Az, dz, m_fractureSkinFactor, fractureAreaWeightedlength); + + double transmissibility = sqrt(transmissibility_X * transmissibility_X + + transmissibility_Y * transmissibility_Y + + transmissibility_Z * transmissibility_Z); + + + m_globalIndeciesToContributingEclipseCells.push_back(fracCell); + m_contributingEclipseCellTransmissibilities.push_back(transmissibility); + + } + +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector EclipseToStimPlanCellTransmissibilityCalculator::getPotentiallyFracturedCellsForPolygon(std::vector polygon) +{ + std::vector cellindecies; + + RiaApplication* app = RiaApplication::instance(); + RimView* activeView = RiaApplication::instance()->activeReservoirView(); + if (!activeView) return cellindecies; + + RimEclipseView* activeRiv = dynamic_cast(activeView); + if (!activeRiv) return cellindecies; + + const RigMainGrid* mainGrid = activeRiv->mainGrid(); + if (!mainGrid) return cellindecies; + + cvf::BoundingBox polygonBBox; + for (cvf::Vec3d nodeCoord : polygon) polygonBBox.add(nodeCoord); + + mainGrid->findIntersectingCells(polygonBBox, &cellindecies); + + return cellindecies; +} + + + + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool EclipseToStimPlanCellTransmissibilityCalculator::planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, + cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ) +{ + + cvf::Plane fracturePlane; + bool isCellIntersected = false; + + fracturePlane.setFromPointAndNormal(static_cast(m_fractureTransform.translation()), + static_cast(m_fractureTransform.col(2))); + + const RigMainGrid* mainGrid = m_case->eclipseCaseData()->mainGrid(); + if (!mainGrid) return false; + + RigCell cell = mainGrid->globalCellArray()[cellindex]; + if (cell.isInvalid()) return mainGrid; + + //Copied (and adapted) from RigEclipseWellLogExtractor + cvf::Vec3d hexCorners[8]; + const std::vector& nodeCoords = mainGrid->nodes(); + const caf::SizeTArray8& cornerIndices = cell.cornerIndices(); + + hexCorners[0] = nodeCoords[cornerIndices[0]]; + hexCorners[1] = nodeCoords[cornerIndices[1]]; + hexCorners[2] = nodeCoords[cornerIndices[2]]; + hexCorners[3] = nodeCoords[cornerIndices[3]]; + hexCorners[4] = nodeCoords[cornerIndices[4]]; + hexCorners[5] = nodeCoords[cornerIndices[5]]; + hexCorners[6] = nodeCoords[cornerIndices[6]]; + hexCorners[7] = nodeCoords[cornerIndices[7]]; + + //Find line-segments where cell and fracture plane intersects + std::list > intersectionLineSegments; + + isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments); + + RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons); + + RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ); + + return isCellIntersected; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double EclipseToStimPlanCellTransmissibilityCalculator::calculateMatrixTransmissibility(double perm, double NTG, double A, double cellSizeLength, double skinfactor, double fractureAreaWeightedlength) +{ + double transmissibility; double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D; - transmissibility = 8 * cDarcy() * (perm * NTG) * A / (cellSizeLength + slDivPi); + transmissibility = 8 * m_cDarcy * (perm * NTG) * A / (cellSizeLength + slDivPi); return transmissibility; } diff --git a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h index 25bf52fdeb..e27c52dcd3 100644 --- a/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h +++ b/ApplicationCode/ReservoirDataModel/RigFractureTransCalc.h @@ -69,28 +69,43 @@ public: 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); private: RimEclipseCase* m_case; RimFracture* m_fracture; RimDefines::UnitSystem m_unitForCalculation; - double calculateMatrixTransmissibility(double permX, double NTG, double Ay, double dx, double skinfactor, double fractureAreaWeightedlength); double cDarcy(); }; class EclipseToStimPlanCellTransmissibilityCalculator { public: - explicit EclipseToStimPlanCellTransmissibilityCalculator(const RimEclipseCase* caseToApply, cvf::Mat4f fractureTransform, double skinFactor, const RigStimPlanFracTemplateCell& stimPlanCell); + explicit EclipseToStimPlanCellTransmissibilityCalculator(const RimEclipseCase* caseToApply, + cvf::Mat4f fractureTransform, + double skinFactor, + double cDarcy, + const RigStimPlanFracTemplateCell& stimPlanCell); const std::vector& globalIndeciesToContributingEclipseCells(); const std::vector& contributingEclipseCellTransmissibilities(); private: - void calculateStimPlanCellsMatrixTransmissibility(); - - std::vector m_globalIndeciesToContributingEclipseCells; - std::vector m_contributingEclipseCellTransmissibilities; + void calculateStimPlanCellsMatrixTransmissibility(); + static std::vector getPotentiallyFracturedCellsForPolygon(std::vector polygon); + bool planeCellIntersectionPolygons(size_t cellindex, std::vector > & polygons, + cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ); + double calculateMatrixTransmissibility(double permX, double NTG, double Ay, double dx, double skinfactor, double fractureAreaWeightedlength); + + + const RimEclipseCase* m_case; + double m_cDarcy; + double m_fractureSkinFactor; + cvf::Mat4f m_fractureTransform; + const RigStimPlanFracTemplateCell& m_stimPlanCell; + + std::vector m_globalIndeciesToContributingEclipseCells; + std::vector m_contributingEclipseCellTransmissibilities; }; \ No newline at end of file diff --git a/ApplicationCode/ReservoirDataModel/RigStimPlanFracTemplateCell.h b/ApplicationCode/ReservoirDataModel/RigStimPlanFracTemplateCell.h index 1c4c422639..a9f0ea4f7e 100644 --- a/ApplicationCode/ReservoirDataModel/RigStimPlanFracTemplateCell.h +++ b/ApplicationCode/ReservoirDataModel/RigStimPlanFracTemplateCell.h @@ -37,7 +37,7 @@ public: virtual ~RigStimPlanFracTemplateCell(); - std::vector getPolygon() { return m_polygon; } + std::vector getPolygon() const { return m_polygon; } double getConductivtyValue() const { return m_concutivityValue; } double getDisplayValue() { return m_displayValue; } size_t getI() { return m_i; }