#6064 Compute stress and stress gradients for fracture model plot.

ApplicationCode/ProjectDataModel/RimFractureModelPlot.cpp

@@ -18,6 +18,8 @@
 #include "RimFractureModelPlot.h"
 
 #include "RiaDefines.h"
+#include "RiaLogging.h"
+
 #include "RicfCommandObject.h"
 
 #include "RimEclipseCase.h"
@@ -40,6 +42,18 @@ CAF_PDM_SOURCE_INIT( RimFractureModelPlot, "FractureModelPlot" );
 RimFractureModelPlot::RimFractureModelPlot()
 {
     CAF_PDM_InitScriptableObject( "Fracture Model Plot", "", "", "A fracture model plot" );
+
+    CAF_PDM_InitFieldNoDefault( &m_fractureModel, "FractureModel", "Fracture Model", "", "", "" );
+    m_fractureModel.uiCapability()->setUiTreeChildrenHidden( true );
+    m_fractureModel.uiCapability()->setUiHidden( true );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+void RimFractureModelPlot::setFractureModel( RimFractureModel* fractureModel )
+{
+    m_fractureModel = fractureModel;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -75,6 +89,30 @@ void RimFractureModelPlot::applyDataSource()
     this->updateConnectedEditors();
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFractureModelPlot::convertToPsiPerFeetFromBarPerMeter( double value )
+{
+    return value * 4.42075;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFractureModelPlot::convertToFeetFromMeter( double value )
+{
+    return value * 3.28084;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFractureModelPlot::convertToPsiFromBar( double value )
+{
+    return value * 14.5038;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -122,11 +160,30 @@ void RimFractureModelPlot::calculateLayers( std::vector<std::pair<double, double
         {
             layerBoundaryDepths.push_back( std::make_pair( depths[startIndex], depths[i] ) );
             layerBoundaryIndexes.push_back( std::make_pair( startIndex, i ) );
-            startIndex = i + 1;
+            startIndex = i;
         }
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+double RimFractureModelPlot::computeValueAtDepth( const std::vector<double>&              values,
+                                                  std::vector<std::pair<double, double>>& layerBoundaryDepths,
+                                                  double                                  depth )
+{
+    for ( size_t i = 0; i < layerBoundaryDepths.size(); i++ )
+    {
+        if ( layerBoundaryDepths[i].first <= depth && layerBoundaryDepths[i].second >= depth )
+        {
+            return values[i];
+        }
+    }
+
+    RiaLogging::error( QString( "Failed to compute value at depth: %1" ).arg( depth ) );
+    return std::numeric_limits<double>::infinity();
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -180,11 +237,10 @@ std::vector<double> RimFractureModelPlot::calculateTrueVerticalDepth() const
     std::vector<double> tvdTopZone;
     for ( auto p : layerBoundaryDepths )
     {
-        std::cout << "Layer boundaries: " << p.first << " - " << p.second << std::endl;
-        tvdTopZone.push_back( p.first );
+        double depthInFeet = convertToFeetFromMeter( p.first );
+        tvdTopZone.push_back( depthInFeet );
     }
 
-    // TODO: convert to feet!!!!
     return tvdTopZone;
 }
 
@@ -196,7 +252,7 @@ std::vector<double> RimFractureModelPlot::findCurveAndComputeLayeredAverage( con
     RimWellLogExtractionCurve* curve = findCurveByName( curveName );
     if ( !curve )
     {
-        std::cerr << "NO " << curveName.toStdString() << " FOUND!!!" << std::endl;
+        RiaLogging::error( QString( "No curve named '%1' found" ).arg( curveName ) );
         return std::vector<double>();
     }
 
@@ -248,7 +304,139 @@ std::vector<double> RimFractureModelPlot::calculateVerticalPermeability() const
 //--------------------------------------------------------------------------------------------------
 std::vector<double> RimFractureModelPlot::calculateStress() const
 {
-    return std::vector<double>();
+    std::vector<double> stress;
+    std::vector<double> stressGradients;
+    calculateStressWithGradients( stress, stressGradients );
+    return stress;
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RimFractureModelPlot::calculateStressWithGradients( std::vector<double>& stress,
+                                                         std::vector<double>& stressGradients ) const
+{
+    // Reference stress
+    const double verticalStressRef         = m_fractureModel->verticalStress();
+    const double verticalStressGradientRef = m_fractureModel->verticalStressGradient();
+    const double stressDepthRef            = m_fractureModel->stressDepth();
+
+    std::vector<std::pair<double, double>> layerBoundaryDepths;
+    std::vector<std::pair<size_t, size_t>> layerBoundaryIndexes;
+    calculateLayers( layerBoundaryDepths, layerBoundaryIndexes );
+
+    // Biot coefficient
+    RimWellLogExtractionCurve* biotCurve = findCurveByName( "Biot Coefficient" );
+    if ( !biotCurve )
+    {
+        RiaLogging::error( "Biot coefficient data not found." );
+        return false;
+    }
+    std::vector<double> biotData = biotCurve->curveData()->xValues();
+
+    // Biot coefficient
+    RimWellLogExtractionCurve* k0Curve = findCurveByName( "k0" );
+    if ( !k0Curve )
+    {
+        RiaLogging::error( "k0 data not found." );
+        return false;
+    }
+    std::vector<double> k0Data = k0Curve->curveData()->xValues();
+
+    // Pressure at the give time step
+    RimWellLogExtractionCurve* timeStepPressureCurve = findCurveByName( "PRESSURE" );
+    if ( !timeStepPressureCurve )
+    {
+        RiaLogging::error( "Pressure data for time step not found." );
+        return false;
+    }
+    std::vector<double> timeStepPressureData = timeStepPressureCurve->curveData()->xValues();
+
+    // Initial pressure
+    RimWellLogExtractionCurve* initialPressureCurve = findCurveByName( "INITIAL PRESSURE" );
+    if ( !initialPressureCurve )
+    {
+        RiaLogging::error( "Initial pressure data not found." );
+        return false;
+    }
+    std::vector<double> initialPressureData = initialPressureCurve->curveData()->xValues();
+
+    // Poissons ratio
+    RimWellLogExtractionCurve* poissonsRatioCurve = findCurveByName( "Poisson's Ratio" );
+    if ( !poissonsRatioCurve )
+    {
+        RiaLogging::error( "Poisson's ratio data not found." );
+        return false;
+    }
+    std::vector<double> poissonsRatioData = poissonsRatioCurve->curveData()->xValues();
+
+    std::vector<double> stressForGradients;
+    std::vector<double> pressureForGradients;
+    std::vector<double> depthForGradients;
+
+    // Calculate the stress
+    for ( size_t i = 0; i < layerBoundaryDepths.size(); i++ )
+    {
+        double depthTopOfZone    = layerBoundaryDepths[i].first;
+        double depthBottomOfZone = layerBoundaryDepths[i].second;
+
+        // Data from curves at the top zone depth
+        double k0               = computeValueAtDepth( k0Data, layerBoundaryDepths, depthTopOfZone );
+        double biot             = computeValueAtDepth( biotData, layerBoundaryDepths, depthTopOfZone );
+        double poissonsRatio    = computeValueAtDepth( poissonsRatioData, layerBoundaryDepths, depthTopOfZone );
+        double initialPressure  = computeValueAtDepth( initialPressureData, layerBoundaryDepths, depthTopOfZone );
+        double timeStepPressure = computeValueAtDepth( timeStepPressureData, layerBoundaryDepths, depthTopOfZone );
+
+        // Vertical stress
+        // Use difference between reference depth and depth of top of zone
+        double depthDiff = depthTopOfZone - stressDepthRef;
+        double Sv        = verticalStressRef + verticalStressGradientRef * depthDiff;
+
+        double Sh_init      = k0 * Sv + initialPressure * ( 1.0 - k0 );
+        double pressureDiff = timeStepPressure - initialPressure;
+
+        // Vertical stress diff assumed to be zero
+        double Sv_diff               = 0.0;
+        double deltaHorizontalStress = poissonsRatio / ( 1.0 - poissonsRatio ) * ( Sv_diff - biot * pressureDiff ) +
+                                       ( biot * pressureDiff );
+
+        double depletionStress = Sh_init + deltaHorizontalStress;
+        stress.push_back( convertToPsiFromBar( depletionStress ) );
+
+        // Cache some results for the gradients calculation
+        stressForGradients.push_back( Sv );
+        pressureForGradients.push_back( initialPressure );
+        depthForGradients.push_back( depthTopOfZone );
+
+        if ( i == layerBoundaryDepths.size() - 1 )
+        {
+            // Use the bottom of the last layer to compute gradient for last layer
+            double bottomInitialPressure =
+                computeValueAtDepth( initialPressureData, layerBoundaryDepths, depthBottomOfZone );
+            double bottomDepthDiff = depthBottomOfZone - stressDepthRef;
+            double bottomSv        = verticalStressRef + verticalStressGradientRef * bottomDepthDiff;
+            stressForGradients.push_back( bottomSv );
+            pressureForGradients.push_back( bottomInitialPressure );
+            depthForGradients.push_back( depthBottomOfZone );
+        }
+    }
+
+    assert( stressForGradients.size() == layerBoundaryDepths.size() + 1 );
+    assert( pressureForGradients.size() == layerBoundaryDepths.size() + 1 );
+    assert( depthForGradients.size() == layerBoundaryDepths.size() + 1 );
+
+    // Second pass to calculate the stress gradients
+    for ( size_t i = 0; i < layerBoundaryDepths.size(); i++ )
+    {
+        double diffStress   = stressForGradients[i + 1] - stressForGradients[i];
+        double diffPressure = pressureForGradients[i + 1] - pressureForGradients[i];
+        double diffDepth    = depthForGradients[i + 1] - depthForGradients[i];
+        double k0           = computeValueAtDepth( k0Data, layerBoundaryDepths, depthForGradients[i] );
+        double gradient     = ( diffStress * k0 + diffPressure * ( 1.0 - k0 ) ) / diffDepth;
+        stressGradients.push_back( convertToPsiPerFeetFromBarPerMeter( gradient ) );
+    }
+
+    return true;
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -256,7 +444,10 @@ std::vector<double> RimFractureModelPlot::calculateStress() const
 //--------------------------------------------------------------------------------------------------
 std::vector<double> RimFractureModelPlot::calculateStressGradient() const
 {
-    return std::vector<double>();
+    std::vector<double> stress;
+    std::vector<double> stressGradients;
+    calculateStressWithGradients( stress, stressGradients );
+    return stressGradients;
 }
 
 //--------------------------------------------------------------------------------------------------