OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6095 from OPM/pore-compressibility-fix-6031

Browse Source 
Pore compressibility fix 6031
This commit is contained in:
Kristian Bendiksen
2020-06-19 11:16:31 +02:00
committed by GitHub
parent 034f75f64e 5276001787
commit 6a54458551
6 changed files with 181 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ApplicationCode/Application/Tools/RiaEclipseUnitTools.h
View File

@@ -45,8 +45,10 @@ public:
static double mmToMeter( double mm ) { return mm / 1000.0; }
static double meterToMm( double meter ) { return 1000.0 * meter; }
static double barToPascal( double bar ) { return bar * 100000.0; }
static double barToPsi( double bar ) { return bar * 14.5038; }
static double barPerMeterToPsiPerFeet( double barPerMeter ) { return barPerMeter * 4.42075; }
static double gigaPascalToPascal( double gigaPascal ) { return gigaPascal * 1.0e9; }
static double darcysConstant( UnitSystem unitSystem );

60
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorPoreCompressibility.cpp
View File

@@ -18,6 +18,7 @@
#include "RigFemPartResultCalculatorPoreCompressibility.h"
#include "RiaEclipseUnitTools.h"
#include "RiaLogging.h"
#include "RigFemPart.h"
@@ -140,16 +141,20 @@ RigFemScalarResultFrames*
frameCountProgress.setNextProgressIncrement( 1u );
int referenceFrameIdx = m_resultCollection->referenceTimeStep();
int frameCount = srcEVDataFrames->frameCount();
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
const std::vector<float>& evData = srcEVDataFrames->frameData( fIdx );
const std::vector<float>& verticalStrainData = verticalStrainDataFrames->frameData( fIdx );
const std::vector<float>& youngsModuliData = youngsModuliFrames->frameData( fIdx );
const std::vector<float>& poissonRatioData = poissonRatioFrames->frameData( fIdx );
const std::vector<float>& voidRatioData = voidRatioFrames->frameData( fIdx );
const std::vector<float>& initialPorFrameData = srcPORDataFrames->frameData( 0 );
const std::vector<float>& porFrameData = srcPORDataFrames->frameData( fIdx );
const std::vector<float>& evData = srcEVDataFrames->frameData( fIdx );
const std::vector<float>& referenceEvData = srcEVDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& verticalStrainData = verticalStrainDataFrames->frameData( fIdx );
const std::vector<float>& referenceVerticalStrainData = verticalStrainDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& youngsModuliData = youngsModuliFrames->frameData( fIdx );
const std::vector<float>& poissonRatioData = poissonRatioFrames->frameData( fIdx );
const std::vector<float>& voidRatioData = voidRatioFrames->frameData( referenceFrameIdx );
const std::vector<float>& referencePorFrameData = srcPORDataFrames->frameData( referenceFrameIdx );
const std::vector<float>& porFrameData = srcPORDataFrames->frameData( fIdx );
std::vector<float>& poreCompressibilityFrameData = poreCompressibilityFrames->frameData( fIdx );
std::vector<float>& verticalCompressibilityFrameData = verticalCompressibilityFrames->frameData( fIdx );
@@ -187,16 +192,16 @@ RigFemScalarResultFrames*
size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
if ( elmNodResIdx < evData.size() )
{
if ( fIdx == 0 )
if ( fIdx == referenceFrameIdx )
{
// Geostatic step: result not defined
// The time step and the reference time step are the same: results undefined
poreCompressibilityFrameData[elmNodResIdx] = inf;
verticalCompressibilityFrameData[elmNodResIdx] = inf;
verticalCompressibilityRatioFrameData[elmNodResIdx] = inf;
}
else
{
// Use biot coefficient for all other (not Geostatic) timesteps
// Use biot coefficient for all timesteps
double biotCoefficient = 1.0;
if ( biotData.empty() )
{
@@ -210,9 +215,10 @@ RigFemScalarResultFrames*
int nodeIdx = femPart->nodeIdxFromElementNodeResultIdx( elmNodResIdx );
// Calculate bulk modulus for solids (grains)
double poissonRatio = poissonRatioData[elmIdx];
double youngsModuli = youngsModuliData[elmIdx];
// Calculate bulk modulus for solids (grains).
// Incoming unit for Young's Modulus is GPa: convert to Pa.
double poissonRatio = poissonRatioData[elmIdx];
double youngsModuli = RiaEclipseUnitTools::gigaPascalToPascal( youngsModuliData[elmIdx] );
double bulkModulusFrame = youngsModuli / ( 3.0 * ( 1.0 - 2.0 * poissonRatio ) );
double bulkModulus = bulkModulusFrame / ( 1.0 - biotCoefficient );
@@ -220,17 +226,19 @@ RigFemScalarResultFrames*
double voidr = voidRatioData[elmNodResIdx];
double porosity = voidr / ( 1.0 + voidr );
// Calculate difference in pore pressure between reference state and this state
double initialPorePressure = initialPorFrameData[nodeIdx];
double framePorePressure = porFrameData[nodeIdx];
double deltaPorePressure = framePorePressure - initialPorePressure;
// Calculate difference in pore pressure between reference state and this state,
// and convert unit from Bar to Pascal.
double referencePorePressure = referencePorFrameData[nodeIdx];
double framePorePressure = porFrameData[nodeIdx];
double deltaPorePressure =
RiaEclipseUnitTools::barToPascal( framePorePressure - referencePorePressure );
// Calculate pore compressibility
double poreCompressibility = inf;
if ( deltaPorePressure != 0.0 && porosity != 0.0 )
{
poreCompressibility =
-( biotCoefficient * evData[elmNodResIdx] ) / ( deltaPorePressure * porosity );
double deltaEv = evData[elmNodResIdx] - referenceEvData[elmNodResIdx];
poreCompressibility = -( biotCoefficient * deltaEv ) / ( deltaPorePressure * porosity );
// Guard against divide by zero: second term can be ignored when bulk modulus is zero,
// which can happens when biot coefficient is 1.0
if ( biotCoefficient != 1.0 && porosity != 1.0 )
@@ -238,16 +246,19 @@ RigFemScalarResultFrames*
poreCompressibility += ( 1.0 / bulkModulus ) * ( biotCoefficient / porosity - 1.0 );
}
}
poreCompressibilityFrameData[elmNodResIdx] = poreCompressibility;
// Convert from 1/Pa to 1/GPa
poreCompressibilityFrameData[elmNodResIdx] = poreCompressibility * 1.0e9;
double verticalCompressibility = inf;
double verticalCompressibilityRatio = inf;
if ( biotCoefficient != 0.0 && deltaPorePressure != 0.0 )
{
// Calculate vertical compressibility
verticalCompressibility =
-verticalStrainData[elmNodResIdx] / ( biotCoefficient * deltaPorePressure );
double deltaStrain = verticalStrainData[elmNodResIdx] -
referenceVerticalStrainData[elmNodResIdx];
// Calculate vertical compressibility (unit: 1/Pa)
verticalCompressibility = -deltaStrain / ( biotCoefficient * deltaPorePressure );
// Calculate vertical compressibility ratio
verticalCompressibilityRatio =
@@ -255,7 +266,8 @@ RigFemScalarResultFrames*
( ( 1.0 + poissonRatio ) * ( 1.0 - 2.0 * poissonRatio ) );
}
verticalCompressibilityFrameData[elmNodResIdx] = verticalCompressibility;
// Convert from 1/Pa to 1/GPa
verticalCompressibilityFrameData[elmNodResIdx] = verticalCompressibility * 1.0e9;
verticalCompressibilityRatioFrameData[elmNodResIdx] = verticalCompressibilityRatio;
}
}

76
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp
View File

@@ -111,6 +111,8 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf
m_biotFixedFactor = 1.0;
m_biotResultAddress = "";
m_referenceTimeStep = 0;
m_resultCalculators.push_back(
std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorTimeLapse( *this ) ) );
m_resultCalculators.push_back(
@@ -395,6 +397,28 @@ void RigFemPartResultsCollection::setBiotCoefficientParameters( double biotFixed
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartResultsCollection::setReferenceTimeStep( int referenceTimeStep )
{
m_referenceTimeStep = referenceTimeStep;
std::set<RigFemResultAddress> results = referenceCaseDependentResults();
for ( auto result : results )
{
deleteResult( result );
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigFemPartResultsCollection::referenceTimeStep() const
{
return m_referenceTimeStep;
}
//--------------------------------------------------------------------------------------------------
/// Will always return a valid object, but it can be empty
//--------------------------------------------------------------------------------------------------
@@ -1150,6 +1174,23 @@ std::vector<RigFemResultAddress>
return addresses;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isResultInSet( const RigFemResultAddress& result,
const std::set<RigFemResultAddress>& results )
{
for ( auto res : results )
{
if ( res.resultPosType == result.resultPosType && res.fieldName == result.fieldName &&
res.componentName == result.componentName )
{
return true;
}
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@@ -1182,15 +1223,36 @@ std::set<RigFemResultAddress> RigFemPartResultsCollection::normalizedResults()
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isNormalizableResult( const RigFemResultAddress& result )
{
for ( auto normRes : normalizedResults() )
return isResultInSet( result, normalizedResults() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::set<RigFemResultAddress> RigFemPartResultsCollection::referenceCaseDependentResults()
{
std::set<RigFemResultAddress> results;
for ( auto elementType : {RIG_ELEMENT_NODAL, RIG_INTEGRATION_POINT} )
{
if ( normRes.resultPosType == result.resultPosType && normRes.fieldName == result.fieldName &&
normRes.componentName == result.componentName )
{
return true;
}
results.insert(
RigFemResultAddress( elementType, "COMPRESSIBILITY", "PORE", RigFemResultAddress::allTimeLapsesValue() ) );
results.insert(
RigFemResultAddress( elementType, "COMPRESSIBILITY", "VERTICAL", RigFemResultAddress::allTimeLapsesValue() ) );
results.insert( RigFemResultAddress( elementType,
"COMPRESSIBILITY",
"VERTICAL-RATIO",
RigFemResultAddress::allTimeLapsesValue() ) );
}
return false;
return results;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultsCollection::isReferenceCaseDependentResult( const RigFemResultAddress& result )
{
return isResultInSet( result, referenceCaseDependentResults() );
}
//--------------------------------------------------------------------------------------------------

10
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h
View File

@@ -120,6 +120,8 @@ public:
double* globalPosClosestToZero,
double* globalNegClosestToZero );
static bool isResultInSet( const RigFemResultAddress& result, const std::set<RigFemResultAddress>& results );
static std::vector<RigFemResultAddress> tensorComponentAddresses( const RigFemResultAddress& resVarAddr );
static std::vector<RigFemResultAddress> tensorPrincipalComponentAdresses( const RigFemResultAddress& resVarAddr );
static std::set<RigFemResultAddress> normalizedResults();
@@ -128,6 +130,12 @@ public:
void setNormalizationAirGap( double normalizationAirGap );
double normalizationAirGap() const;
void setReferenceTimeStep( int referenceTimeStep );
int referenceTimeStep() const;
static std::set<RigFemResultAddress> referenceCaseDependentResults();
static bool isReferenceCaseDependentResult( const RigFemResultAddress& result );
RigFemScalarResultFrames* findOrLoadScalarResult( int partIndex, const RigFemResultAddress& resVarAddr );
RigFemScalarResultFrames* createScalarResult( int partIndex, const RigFemResultAddress& resVarAddr );
@@ -154,6 +162,8 @@ private:
double m_biotFixedFactor;
QString m_biotResultAddress;
int m_referenceTimeStep;
std::vector<std::unique_ptr<RigFemPartResultCalculator>> m_resultCalculators;
RigStatisticsDataCache* statistics( const RigFemResultAddress& resVarAddr );

70
ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp
View File

@@ -93,6 +93,7 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition( void )
"",
"",
"" );
CAF_PDM_InitField( &m_referenceTimeStep, "ReferenceTimeStep", 0, "Reference Time Step", "", "", "" );
CAF_PDM_InitField( &m_compactionRefLayer, "CompactionRefLayer", 0, "Compaction Ref Layer", "", "", "" );
m_compactionRefLayer.uiCapability()->setUiHidden( true );
@@ -170,8 +171,17 @@ void RimGeoMechResultDefinition::defineUiOrdering( QString uiConfigName, caf::Pd
if ( m_resultPositionTypeUiField() != RIG_FORMATION_NAMES )
{
caf::PdmUiGroup* timeLapseGr = uiOrdering.addNewGroup( "Difference Options" );
timeLapseGr->add( &m_timeLapseBaseTimestep );
bool isReferenceCaseDependent = referenceCaseDependentResultSelected();
if ( isReferenceCaseDependent )
{
caf::PdmUiGroup* referenceCaseGr = uiOrdering.addNewGroup( "Reference Case" );
referenceCaseGr->add( &m_referenceTimeStep );
}
else
{
caf::PdmUiGroup* timeLapseGr = uiOrdering.addNewGroup( "Difference Options" );
timeLapseGr->add( &m_timeLapseBaseTimestep );
}
}
if ( m_resultPositionTypeUiField() == RIG_NODAL )
@@ -264,6 +274,21 @@ QList<caf::PdmOptionItemInfo>
static_cast<int>( stepIdx ) ) );
}
}
else if ( &m_referenceTimeStep == fieldNeedingOptions )
{
std::vector<std::string> stepNames;
if ( m_geomCase->geoMechData() )
{
stepNames = m_geomCase->geoMechData()->femPartResults()->filteredStepNames();
}
for ( size_t stepIdx = 0; stepIdx < stepNames.size(); ++stepIdx )
{
options.push_back(
caf::PdmOptionItemInfo( QString( "%1 (#%2)" ).arg( QString::fromStdString( stepNames[stepIdx] ) ).arg( stepIdx ),
static_cast<int>( stepIdx ) ) );
}
}
else if ( &m_compactionRefLayerUiField == fieldNeedingOptions )
{
if ( m_geomCase->geoMechData() )
@@ -319,7 +344,8 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
}
}
if ( &m_resultPositionTypeUiField == changedField || &m_timeLapseBaseTimestep == changedField )
if ( &m_resultPositionTypeUiField == changedField || &m_timeLapseBaseTimestep == changedField ||
&m_referenceTimeStep == changedField )
{
std::map<std::string, std::vector<std::string>> fieldCompNames = getResultMetaDataForUIFieldSetting();
QStringList uiVarNames;
@@ -336,6 +362,11 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
{
m_resultVariableUiField = "";
}
if ( &m_referenceTimeStep == changedField )
{
m_geomCase->geoMechData()->femPartResults()->setReferenceTimeStep( m_referenceTimeStep() );
}
}
if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationAirGap == 0.0 )
{
@@ -357,7 +388,7 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch
if ( &m_resultVariableUiField == changedField || &m_compactionRefLayerUiField == changedField ||
&m_timeLapseBaseTimestep == changedField || &m_normalizeByHydrostaticPressure == changedField ||
&m_normalizationAirGap == changedField )
&m_normalizationAirGap == changedField || &m_referenceTimeStep == changedField )
{
QStringList fieldComponentNames = m_resultVariableUiField().split( QRegExp( "\\s+" ) );
if ( fieldComponentNames.size() > 0 )
@@ -667,9 +698,15 @@ QString RimGeoMechResultDefinition::resultComponentName() const
QString RimGeoMechResultDefinition::diffResultUiName() const
{
QString diffResultString;
if ( m_timeLapseBaseTimestep != RigFemResultAddress::noTimeLapseValue() )
if ( m_geomCase->geoMechData() )
{
if ( m_geomCase->geoMechData() )
if ( referenceCaseDependentResultSelected() )
{
std::vector<std::string> stepNames = m_geomCase->geoMechData()->femPartResults()->filteredStepNames();
QString timeStepString = QString::fromStdString( stepNames[m_referenceTimeStep()] );
diffResultString += QString( "<b>Reference Time Step</b>: %1" ).arg( timeStepString );
}
else if ( m_timeLapseBaseTimestep != RigFemResultAddress::noTimeLapseValue() )
{
std::vector<std::string> stepNames = m_geomCase->geoMechData()->femPartResults()->filteredStepNames();
QString timeStepString = QString::fromStdString( stepNames[m_timeLapseBaseTimestep()] );
@@ -685,9 +722,13 @@ QString RimGeoMechResultDefinition::diffResultUiName() const
QString RimGeoMechResultDefinition::diffResultUiShortName() const
{
QString diffResultString;
if ( m_timeLapseBaseTimestep != RigFemResultAddress::noTimeLapseValue() )
if ( m_geomCase->geoMechData() )
{
if ( m_geomCase->geoMechData() )
if ( referenceCaseDependentResultSelected() )
{
diffResultString += QString( "Ref. Time: #%1" ).arg( m_referenceTimeStep() );
}
else if ( m_timeLapseBaseTimestep != RigFemResultAddress::noTimeLapseValue() )
{
diffResultString += QString( "Base Time: #%1" ).arg( m_timeLapseBaseTimestep() );
}
@@ -782,6 +823,11 @@ QString RimGeoMechResultDefinition::currentResultUnits() const
{
return "GPa";
}
else if ( this->resultFieldName() == "COMPRESSIBILITY" &&
( this->resultComponentName() == "PORE" || this->resultComponentName() == "VERTICAL" ) )
{
return "1/GPa";
}
else
{
for ( auto resultName : RiaDefines::wbsDerivedResultNames() )
@@ -854,6 +900,14 @@ bool RimGeoMechResultDefinition::normalizableResultSelected() const
return RigFemPartResultsCollection::isNormalizableResult( this->resultAddress() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimGeoMechResultDefinition::referenceCaseDependentResultSelected() const
{
return RigFemPartResultsCollection::isReferenceCaseDependentResult( this->resultAddress() );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

2
ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h
View File

@@ -112,6 +112,7 @@ private:
static QString convertToUiResultFieldName( QString resultFieldName );
bool normalizableResultSelected() const;
bool referenceCaseDependentResultSelected() const;
// Data Fields
@@ -122,6 +123,7 @@ private:
caf::PdmField<int> m_compactionRefLayer;
caf::PdmField<bool> m_normalizeByHydrostaticPressure;
caf::PdmField<double> m_normalizationAirGap;
caf::PdmField<int> m_referenceTimeStep;
// UI Fields only
friend class RimGeoMechPropertyFilter; // Property filter needs the ui fields