mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
#1487 Move the EclToStPlCellTransCalc to a new file
This commit is contained in:
Jacob Støren
parent
ab121fe80a
commit
cac9a9ed03
@ -23,8 +23,10 @@
|
||||
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigFracture.h"
|
||||
#include "RigFracture.h"
|
||||
#include "RigFractureTransCalc.h"
|
||||
#include "RigStimPlanFracTemplateCell.h"
|
||||
#include "RigMainGrid.h"
|
||||
#include "RigEclipseToStimPlanCellTransmissibilityCalculator.h"
|
||||
|
||||
#include "RimEclipseCase.h"
|
||||
#include "RimEclipseResultDefinition.h"
|
||||
@ -32,6 +34,9 @@
|
||||
#include "RimEclipseWell.h"
|
||||
#include "RimEllipseFractureTemplate.h"
|
||||
#include "RimFracture.h"
|
||||
#include "RimFractureTemplate.h"
|
||||
#include "RimSimWellFracture.h"
|
||||
#include "RimStimPlanFractureTemplate.h"
|
||||
#include "RimWellPath.h"
|
||||
|
||||
#include "cafProgressInfo.h"
|
||||
@ -39,11 +44,6 @@
|
||||
#include <QFile>
|
||||
#include <QString>
|
||||
#include <QTextStream>
|
||||
#include "RigFractureTransCalc.h"
|
||||
#include "RimFractureTemplate.h"
|
||||
#include "RimStimPlanFractureTemplate.h"
|
||||
#include "RigStimPlanFracTemplateCell.h"
|
||||
#include "RimSimWellFracture.h"
|
||||
|
||||
|
||||
|
||||
|
||||
@ -46,6 +46,7 @@ ${CEE_CURRENT_LIST_DIR}RigTransmissibilityCondenser.h
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseNativeStatCalc.h
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseNativeVisibleCellsStatCalc.h
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseMultiPropertyStatCalc.h
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseToStimPlanCellTransmissibilityCalculator.h
|
||||
${CEE_CURRENT_LIST_DIR}RigWellLogCurveData.h
|
||||
${CEE_CURRENT_LIST_DIR}RigTimeHistoryResultAccessor.h
|
||||
${CEE_CURRENT_LIST_DIR}RigCurveDataTools.h
|
||||
@ -60,7 +61,6 @@ ${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.h
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFracTemplateCell.h
|
||||
|
||||
|
||||
|
||||
)
|
||||
|
||||
set (SOURCE_GROUP_SOURCE_FILES
|
||||
@ -101,6 +101,7 @@ ${CEE_CURRENT_LIST_DIR}RigTransmissibilityCondenser.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseNativeStatCalc.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseNativeVisibleCellsStatCalc.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseMultiPropertyStatCalc.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigEclipseToStimPlanCellTransmissibilityCalculator.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigWellLogCurveData.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigTimeHistoryResultAccessor.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigCurveDataTools.cpp
|
||||
@ -114,7 +115,6 @@ ${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFractureDefinition.cpp
|
||||
${CEE_CURRENT_LIST_DIR}RigStimPlanFracTemplateCell.cpp
|
||||
|
||||
|
||||
)
|
||||
|
||||
list(APPEND CODE_HEADER_FILES
|
||||
|
||||
@ -0,0 +1,316 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2017 Statoil ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#include "RigEclipseToStimPlanCellTransmissibilityCalculator.h"
|
||||
|
||||
#include "RigStimPlanFracTemplateCell.h"
|
||||
#include "RigResultAccessorFactory.h"
|
||||
#include "RigEclipseCaseData.h"
|
||||
#include "RigActiveCellInfo.h"
|
||||
#include "RigCellGeometryTools.h"
|
||||
#include "RigMainGrid.h"
|
||||
|
||||
#include "RimEclipseCase.h"
|
||||
#include "RimEclipseView.h"
|
||||
#include "RimReservoirCellResultsStorage.h"
|
||||
|
||||
#include "RiaApplication.h"
|
||||
|
||||
#include "cvfGeometryTools.h"
|
||||
|
||||
|
||||
|
||||
|
||||
RigEclipseToStimPlanCellTransmissibilityCalculator::RigEclipseToStimPlanCellTransmissibilityCalculator(RimEclipseCase* caseToApply,
|
||||
cvf::Mat4f fractureTransform,
|
||||
double skinFactor,
|
||||
double cDarcy,
|
||||
const RigStimPlanFracTemplateCell& stimPlanCell)
|
||||
|
||||
: m_stimPlanCell(stimPlanCell)
|
||||
{
|
||||
m_case = caseToApply;
|
||||
m_fractureSkinFactor = skinFactor;
|
||||
m_cDarcy = cDarcy;
|
||||
m_fractureTransform = fractureTransform;
|
||||
}
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
const std::vector<size_t>& RigEclipseToStimPlanCellTransmissibilityCalculator::globalIndeciesToContributingEclipseCells()
|
||||
{
|
||||
if (m_globalIndeciesToContributingEclipseCells.size() < 1)
|
||||
{
|
||||
calculateStimPlanCellsMatrixTransmissibility();
|
||||
}
|
||||
|
||||
return m_globalIndeciesToContributingEclipseCells;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
const std::vector<double>& RigEclipseToStimPlanCellTransmissibilityCalculator::contributingEclipseCellTransmissibilities()
|
||||
{
|
||||
if (m_globalIndeciesToContributingEclipseCells.size() < 1)
|
||||
{
|
||||
calculateStimPlanCellsMatrixTransmissibility();
|
||||
}
|
||||
|
||||
return m_contributingEclipseCellTransmissibilities;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsMatrixTransmissibility()
|
||||
{
|
||||
|
||||
//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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<cvf::Vec3d> stimPlanPolygonTransformed;
|
||||
for (cvf::Vec3d v : m_stimPlanCell.getPolygon())
|
||||
{
|
||||
cvf::Vec3f stimPlanPolygonNode = static_cast<cvf::Vec3f>(v);
|
||||
stimPlanPolygonNode.transformPoint(m_fractureTransform);
|
||||
stimPlanPolygonTransformed.push_back(static_cast<cvf::Vec3d>(stimPlanPolygonNode));
|
||||
}
|
||||
|
||||
|
||||
std::vector<size_t> fracCells = getPotentiallyFracturedCellsForPolygon(stimPlanPolygonTransformed);
|
||||
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<std::vector<cvf::Vec3d> > 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<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
for (cvf::Vec3d& v : planeCellPolygon)
|
||||
{
|
||||
v.transformPoint(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
}
|
||||
}
|
||||
|
||||
cvf::Vec3d localZinFracPlane;
|
||||
localZinFracPlane = localZ;
|
||||
localZinFracPlane.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
|
||||
directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
|
||||
directionOfLength.normalize();
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > polygonsForStimPlanCellInEclipseCell;
|
||||
cvf::Vec3d areaVector;
|
||||
std::vector<cvf::Vec3d> stimPlanPolygon = m_stimPlanCell.getPolygon();
|
||||
|
||||
for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, stimPlanPolygon);
|
||||
for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
|
||||
{
|
||||
polygonsForStimPlanCellInEclipseCell.push_back(clippedPolygon);
|
||||
}
|
||||
}
|
||||
|
||||
if (polygonsForStimPlanCellInEclipseCell.size() == 0) continue;
|
||||
|
||||
double area;
|
||||
std::vector<double> areaOfFractureParts;
|
||||
double length;
|
||||
std::vector<double> lengthXareaOfFractureParts;
|
||||
double Ax = 0.0, Ay = 0.0, Az = 0.0;
|
||||
|
||||
for (std::vector<cvf::Vec3d> 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<cvf::Vec3d>(m_fractureTransform.translation()),
|
||||
static_cast<cvf::Vec3d>(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<size_t> RigEclipseToStimPlanCellTransmissibilityCalculator::getPotentiallyFracturedCellsForPolygon(std::vector<cvf::Vec3d> polygon)
|
||||
{
|
||||
std::vector<size_t> cellindecies;
|
||||
|
||||
RiaApplication* app = RiaApplication::instance();
|
||||
RimView* activeView = RiaApplication::instance()->activeReservoirView();
|
||||
if (!activeView) return cellindecies;
|
||||
|
||||
RimEclipseView* activeRiv = dynamic_cast<RimEclipseView*>(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 RigEclipseToStimPlanCellTransmissibilityCalculator::planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons,
|
||||
cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ)
|
||||
{
|
||||
|
||||
cvf::Plane fracturePlane;
|
||||
bool isCellIntersected = false;
|
||||
|
||||
fracturePlane.setFromPointAndNormal(static_cast<cvf::Vec3d>(m_fractureTransform.translation()),
|
||||
static_cast<cvf::Vec3d>(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<cvf::Vec3d>& 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<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
|
||||
|
||||
isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
|
||||
|
||||
RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons);
|
||||
|
||||
RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ);
|
||||
|
||||
return isCellIntersected;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigEclipseToStimPlanCellTransmissibilityCalculator::calculateMatrixTransmissibility(double perm, double NTG, double A, double cellSizeLength, double skinfactor, double fractureAreaWeightedlength)
|
||||
{
|
||||
double transmissibility;
|
||||
|
||||
double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D;
|
||||
transmissibility = 8 * m_cDarcy * (perm * NTG) * A / (cellSizeLength + slDivPi);
|
||||
|
||||
return transmissibility;
|
||||
}
|
||||
@ -0,0 +1,63 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2017 Statoil ASA
|
||||
//
|
||||
// ResInsight is free software: you can redistribute it and/or modify
|
||||
// it under the terms of the GNU General Public License as published by
|
||||
// the Free Software Foundation, either version 3 of the License, or
|
||||
// (at your option) any later version.
|
||||
//
|
||||
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
// FITNESS FOR A PARTICULAR PURPOSE.
|
||||
//
|
||||
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
|
||||
// for more details.
|
||||
//
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <vector>
|
||||
|
||||
#include "cvfBase.h"
|
||||
#include "cvfMatrix4.h"
|
||||
|
||||
class RimEclipseCase;
|
||||
class RigStimPlanFracTemplateCell;
|
||||
|
||||
//==================================================================================================
|
||||
///
|
||||
//==================================================================================================
|
||||
|
||||
class RigEclipseToStimPlanCellTransmissibilityCalculator
|
||||
{
|
||||
public:
|
||||
explicit RigEclipseToStimPlanCellTransmissibilityCalculator(RimEclipseCase* caseToApply,
|
||||
cvf::Mat4f fractureTransform,
|
||||
double skinFactor,
|
||||
double cDarcy,
|
||||
const RigStimPlanFracTemplateCell& stimPlanCell);
|
||||
|
||||
const std::vector<size_t>& globalIndeciesToContributingEclipseCells();
|
||||
const std::vector<double>& contributingEclipseCellTransmissibilities();
|
||||
|
||||
private:
|
||||
void calculateStimPlanCellsMatrixTransmissibility();
|
||||
static std::vector<size_t> getPotentiallyFracturedCellsForPolygon(std::vector<cvf::Vec3d> polygon);
|
||||
bool planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons,
|
||||
cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ);
|
||||
double calculateMatrixTransmissibility(double permX, double NTG, double Ay, double dx, double skinfactor, double fractureAreaWeightedlength);
|
||||
|
||||
|
||||
RimEclipseCase* m_case;
|
||||
double m_cDarcy;
|
||||
double m_fractureSkinFactor;
|
||||
cvf::Mat4f m_fractureTransform;
|
||||
const RigStimPlanFracTemplateCell& m_stimPlanCell;
|
||||
|
||||
std::vector<size_t> m_globalIndeciesToContributingEclipseCells;
|
||||
std::vector<double> m_contributingEclipseCellTransmissibilities;
|
||||
|
||||
};
|
||||
|
||||
@ -869,283 +869,3 @@ double RigFractureTransCalc::convertConductivtyValue(double Kw, RimDefines::Unit
|
||||
return cvf::UNDEFINED_DOUBLE;
|
||||
}
|
||||
|
||||
|
||||
|
||||
RigEclipseToStimPlanCellTransmissibilityCalculator::RigEclipseToStimPlanCellTransmissibilityCalculator(RimEclipseCase* caseToApply,
|
||||
cvf::Mat4f fractureTransform,
|
||||
double skinFactor,
|
||||
double cDarcy,
|
||||
const RigStimPlanFracTemplateCell& stimPlanCell)
|
||||
|
||||
: m_stimPlanCell(stimPlanCell)
|
||||
{
|
||||
m_case = caseToApply;
|
||||
m_fractureSkinFactor = skinFactor;
|
||||
m_cDarcy = cDarcy;
|
||||
m_fractureTransform = fractureTransform;
|
||||
}
|
||||
|
||||
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
const std::vector<size_t>& RigEclipseToStimPlanCellTransmissibilityCalculator::globalIndeciesToContributingEclipseCells()
|
||||
{
|
||||
if (m_globalIndeciesToContributingEclipseCells.size() < 1)
|
||||
{
|
||||
calculateStimPlanCellsMatrixTransmissibility();
|
||||
}
|
||||
|
||||
return m_globalIndeciesToContributingEclipseCells;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
const std::vector<double>& RigEclipseToStimPlanCellTransmissibilityCalculator::contributingEclipseCellTransmissibilities()
|
||||
{
|
||||
if (m_globalIndeciesToContributingEclipseCells.size() < 1)
|
||||
{
|
||||
calculateStimPlanCellsMatrixTransmissibility();
|
||||
}
|
||||
|
||||
return m_contributingEclipseCellTransmissibilities;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsMatrixTransmissibility()
|
||||
{
|
||||
|
||||
//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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<RigResultAccessor> 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<cvf::Vec3d> stimPlanPolygonTransformed;
|
||||
for (cvf::Vec3d v : m_stimPlanCell.getPolygon())
|
||||
{
|
||||
cvf::Vec3f stimPlanPolygonNode = static_cast<cvf::Vec3f>(v);
|
||||
stimPlanPolygonNode.transformPoint(m_fractureTransform);
|
||||
stimPlanPolygonTransformed.push_back(static_cast<cvf::Vec3d>(stimPlanPolygonNode));
|
||||
}
|
||||
|
||||
|
||||
std::vector<size_t> fracCells = getPotentiallyFracturedCellsForPolygon(stimPlanPolygonTransformed);
|
||||
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<std::vector<cvf::Vec3d> > 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<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
for (cvf::Vec3d& v : planeCellPolygon)
|
||||
{
|
||||
v.transformPoint(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
}
|
||||
}
|
||||
|
||||
cvf::Vec3d localZinFracPlane;
|
||||
localZinFracPlane = localZ;
|
||||
localZinFracPlane.transformVector(static_cast<cvf::Mat4d>(invertedTransMatrix));
|
||||
cvf::Vec3d directionOfLength = cvf::Vec3d::ZERO;
|
||||
directionOfLength.cross(localZinFracPlane, cvf::Vec3d(0, 0, 1));
|
||||
directionOfLength.normalize();
|
||||
|
||||
std::vector<std::vector<cvf::Vec3d> > polygonsForStimPlanCellInEclipseCell;
|
||||
cvf::Vec3d areaVector;
|
||||
std::vector<cvf::Vec3d> stimPlanPolygon = m_stimPlanCell.getPolygon();
|
||||
|
||||
for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
|
||||
{
|
||||
std::vector<std::vector<cvf::Vec3d> >clippedPolygons = RigCellGeometryTools::clipPolygons(planeCellPolygon, stimPlanPolygon);
|
||||
for (std::vector<cvf::Vec3d> clippedPolygon : clippedPolygons)
|
||||
{
|
||||
polygonsForStimPlanCellInEclipseCell.push_back(clippedPolygon);
|
||||
}
|
||||
}
|
||||
|
||||
if (polygonsForStimPlanCellInEclipseCell.size() == 0) continue;
|
||||
|
||||
double area;
|
||||
std::vector<double> areaOfFractureParts;
|
||||
double length;
|
||||
std::vector<double> lengthXareaOfFractureParts;
|
||||
double Ax = 0.0, Ay = 0.0, Az = 0.0;
|
||||
|
||||
for (std::vector<cvf::Vec3d> 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<cvf::Vec3d>(m_fractureTransform.translation()),
|
||||
static_cast<cvf::Vec3d>(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<size_t> RigEclipseToStimPlanCellTransmissibilityCalculator::getPotentiallyFracturedCellsForPolygon(std::vector<cvf::Vec3d> polygon)
|
||||
{
|
||||
std::vector<size_t> cellindecies;
|
||||
|
||||
RiaApplication* app = RiaApplication::instance();
|
||||
RimView* activeView = RiaApplication::instance()->activeReservoirView();
|
||||
if (!activeView) return cellindecies;
|
||||
|
||||
RimEclipseView* activeRiv = dynamic_cast<RimEclipseView*>(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 RigEclipseToStimPlanCellTransmissibilityCalculator::planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons,
|
||||
cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ)
|
||||
{
|
||||
|
||||
cvf::Plane fracturePlane;
|
||||
bool isCellIntersected = false;
|
||||
|
||||
fracturePlane.setFromPointAndNormal(static_cast<cvf::Vec3d>(m_fractureTransform.translation()),
|
||||
static_cast<cvf::Vec3d>(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<cvf::Vec3d>& 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<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
|
||||
|
||||
isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
|
||||
|
||||
RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons);
|
||||
|
||||
RigCellGeometryTools::findCellLocalXYZ(hexCorners, localX, localY, localZ);
|
||||
|
||||
return isCellIntersected;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigEclipseToStimPlanCellTransmissibilityCalculator::calculateMatrixTransmissibility(double perm, double NTG, double A, double cellSizeLength, double skinfactor, double fractureAreaWeightedlength)
|
||||
{
|
||||
double transmissibility;
|
||||
|
||||
double slDivPi = (skinfactor * fractureAreaWeightedlength) / cvf::PI_D;
|
||||
transmissibility = 8 * m_cDarcy * (perm * NTG) * A / (cellSizeLength + slDivPi);
|
||||
|
||||
return transmissibility;
|
||||
}
|
||||
|
||||
@ -80,34 +80,3 @@ private:
|
||||
RimDefines::UnitSystem m_unitForCalculation;
|
||||
|
||||
};
|
||||
|
||||
class RigEclipseToStimPlanCellTransmissibilityCalculator
|
||||
{
|
||||
public:
|
||||
explicit RigEclipseToStimPlanCellTransmissibilityCalculator(RimEclipseCase* caseToApply,
|
||||
cvf::Mat4f fractureTransform,
|
||||
double skinFactor,
|
||||
double cDarcy,
|
||||
const RigStimPlanFracTemplateCell& stimPlanCell);
|
||||
|
||||
const std::vector<size_t>& globalIndeciesToContributingEclipseCells();
|
||||
const std::vector<double>& contributingEclipseCellTransmissibilities();
|
||||
|
||||
private:
|
||||
void calculateStimPlanCellsMatrixTransmissibility();
|
||||
static std::vector<size_t> getPotentiallyFracturedCellsForPolygon(std::vector<cvf::Vec3d> polygon);
|
||||
bool planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons,
|
||||
cvf::Vec3d & localX, cvf::Vec3d & localY, cvf::Vec3d & localZ);
|
||||
double calculateMatrixTransmissibility(double permX, double NTG, double Ay, double dx, double skinfactor, double fractureAreaWeightedlength);
|
||||
|
||||
|
||||
RimEclipseCase* m_case;
|
||||
double m_cDarcy;
|
||||
double m_fractureSkinFactor;
|
||||
cvf::Mat4f m_fractureTransform;
|
||||
const RigStimPlanFracTemplateCell& m_stimPlanCell;
|
||||
|
||||
std::vector<size_t> m_globalIndeciesToContributingEclipseCells;
|
||||
std::vector<double> m_contributingEclipseCellTransmissibilities;
|
||||
|
||||
};
|
||||
Loading…
Reference in New Issue
Block a user