diff --git a/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml index fa6d45549b0..2d9209a8872 100644 --- a/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml +++ b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml @@ -54,7 +54,8 @@ solidInternalEnergyModelName="rockInternalEnergy" /> - + + defaultDrainedLinearTEC="3e-7"/> + - + + defaultDrainedLinearTEC="3e-7"/> + + + + + diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/AnalyticalSol.py b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/AnalyticalSol.py new file mode 100755 index 00000000000..1eaf9bcd592 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/AnalyticalSol.py @@ -0,0 +1,50 @@ +import math +import numpy as np + +class DPAnalyticalSolution: + def __init__(self): + self.bulkMod = 0.5e9 # Pa + self.shearMod = 0.3e9 # Pa + self.lambdaCoeff = self.bulkMod - 2*self.shearMod/3. # Pa + self.cohesion = 5e3 # Pa + self.frictionAngle = 15.27 # deg + self.thermalExpansionCoeff = 3e-7 + self.phi = 6 * np.sin(np.deg2rad(self.frictionAngle)) / ( 3 - np.sin(np.deg2rad(self.frictionAngle)) ) + self.C = 6 * self.cohesion * np.cos(np.deg2rad(self.frictionAngle)) / ( 3 - np.sin(np.deg2rad(self.frictionAngle)) ) + + def compute_stress(self, deltaTemp): + sigma_x = 0 + sigma_z = 0 + + epsMech_y = - self.thermalExpansionCoeff * deltaTemp + + # We first assume all elastic deformation + epsMech_x = -epsMech_y * self.lambdaCoeff / 2 / (self.lambdaCoeff + self.shearMod) + epsMech_z = epsMech_x + + sigma_y = (self.lambdaCoeff + 2*self.shearMod) * epsMech_y + 2 * self.lambdaCoeff * epsMech_x + + P = sigma_y / 3. + Q = sigma_y + + yieldFunc = Q + self.phi * P - self.C + + if yieldFunc <= 0.0: + eps_x = epsMech_x + self.thermalExpansionCoeff * deltaTemp + + return sigma_x, sigma_y, eps_x + else: + sigma_y = self.C / (1 + self.phi/3.0) + epsMech_x = (sigma_y - (3*self.lambdaCoeff + 2*self.shearMod)*epsMech_y )/ (6*self.lambdaCoeff + 4*self.shearMod) + multipler = -(sigma_y - (self.lambdaCoeff + 2*self.shearMod)*epsMech_y - 2*self.lambdaCoeff*epsMech_x)/ (4/3 * self.shearMod) + eps_x = epsMech_x + self.thermalExpansionCoeff * deltaTemp + + return sigma_x, sigma_y, eps_x + + def compute_disp(self, deltaTemp): + _, _, esp_x = self.compute_stress(deltaTemp) + + disp_y = 0 + disp_x = esp_x*1 + + return disp_x, disp_y \ No newline at end of file diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/Example.rst b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/Example.rst new file mode 100755 index 00000000000..74bbe3fc852 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/Example.rst @@ -0,0 +1,224 @@ +.. _ExampleThermoMech1DCooling: + + +#################################################### +Thermally Induced Failure under Confined Cooling +#################################################### + +**Context** + + +When a rock cools but cannot contract, it develops tensile stress. In a confined +rock mass, this thermal stress combines with the existing in-situ stresses and +may cause the rock to fail—even when pore pressure and external loading stay constant. +This mechanism is responsible for thermal fracturing around cold-fluid injectors +and caprock damage during CO\ :sub:`2` storage, among other phenomena. + +This example simulates thermal contraction stresses using a one-dimensional +thermo-mechanical problem. We compare two constitutive models subjected to the same cooling history: + +- **Thermo-elastic rock** (``ElasticIsotropic``): The induced stress increases indefinitely as the rock continues cooling. + +- **Thermo-plastic rock** (``DruckerPrager``): The induced stress is limited by the yield surface. Beyond a critical temperature drop, the rock fails and deforms plastically. + +Both cases have closed-form solutions, making this example useful for verifying +the thermo-mechanical coupling and the Drucker-Prager return mapping under thermal loading. + + +**InputFile** + +This example uses no external input files. Everything required is contained within three GEOS +input files located at: + +.. code-block:: console + + inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml + +.. code-block:: console + + inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_fim_smoke.xml + +.. code-block:: console + + inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml + +--------------------------------------------------- +Description of the case +--------------------------------------------------- + +We consider a seven-meter column discretized with 14 elements along the ``y`` direction, and a +single element in the two other directions. The column is initially at a uniform temperature +of 100 K, and is cooled down to 20 K following a linear ramp imposed over the whole domain. + +The mechanical boundary conditions are the essential ingredient of the problem: the two ends +of the column (``yneg`` and ``ypos``) are fixed along ``y``, so the axial strain is prevented +(ε_yy = 0), while the `x` and `z` directions are only restrained on one face each and are +therefore free to deform. The lateral faces being traction-free, σ_xx = σ_zz = 0, and the only +non-zero stress component is σ_yy. + +.. math:: + + \varepsilon_{yy} = 0, \qquad \sigma_{xx} = \sigma_{zz} = 0 + +.. _thermoMech1DCoolingSketchFig: +.. figure:: xz_cross_section_uniaxial_stress.png + :align: center + :width: 500 + :figclass: align-center + + Sketch of the confined column: both ends are fixed along ``y``, while ``x`` and ``z`` are + restrained on a single face each and remain free to deform. + +.. literalinclude:: ../../../../../../../inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml + :language: xml + :start-after: + :end-before: + +The cooling history is prescribed by a ``TableFunction`` applied to the temperature field. +A short initial temperature plateau lets the mechanical equilibrium settle before the thermal +loading starts. + +.. literalinclude:: ../../../../../../../inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml + :language: xml + :start-after: + :end-before: + +The pore pressure is fixed to zero and the permeability is set to a negligible value +(:math:`10^{-100}` m\ :sup:`2`), so that no fluid flow takes place: the stress evolution is +entirely thermo-mechanical, like in the analytical solution. + +------------------------------------------------------------------ +Constitutive models +------------------------------------------------------------------ + +The two cases differ **only** by the solid model. The thermo-elastic case uses an +``ElasticIsotropic`` solid, with a drained linear thermal expansion coefficient +:math:`\alpha = 3 \times 10^{-7}` K\ :sup:`-1`: + +.. literalinclude:: ../../../../../../../inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_fim_smoke.xml + :language: xml + :start-after: + :end-before: + +The thermo-plastic case uses the same elastic properties and thermal expansion +coefficient, and adds a Drucker-Prager yield surface: + +.. literalinclude:: ../../../../../../../inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml + :language: xml + :start-after: + :end-before: + +------------------------------------------------------------------ +Analytical solution +------------------------------------------------------------------ + +**Thermo-elastic response.** In this uniaxial stress state (σ_xx = σ_zz = 0), the axial strain is blocked (ε_yy = 0) and the thermo-elastic constitutive law reduces to + +.. math:: + + \sigma_{yy} = -E \, \alpha \, \Delta T + +with :math:`E = 9KG/(3K+G)` the Young modulus. A cooling :math:`\Delta T < 0` therefore +produces a **tensile** stress that grows linearly with the temperature drop, without any +bound. + +**Onset of failure.** The Drucker-Prager yield function implemented in GEOS reads + +.. math:: + + F = Q + b \, P - c + +where :math:`P = \mathrm{tr}(\sigma)/3` is the mean stress and :math:`Q` the von Mises stress. +The two coefficients are obtained from the friction angle :math:`\varphi` and the cohesion so +that the cone passes through the triaxial compression corners of the Mohr-Coulomb surface: + +.. math:: + + b = \frac{6 \sin \varphi}{3 - \sin \varphi}, \qquad + c = \frac{6\, \mathrm{cohesion} \cos \varphi}{3 - \sin \varphi} + +For the uniaxial stress state of this problem, :math:`Q = \sigma_{yy}` and +:math:`P = \sigma_{yy}/3`, so the yield condition :math:`F = 0` gives a closed-form cap on the +thermally induced stress, and the corresponding critical cooling: + +.. math:: + + \sigma_{f} = \frac{c}{1 + b/3}, \qquad + \Delta T_{f} = -\frac{\sigma_{f}}{E \, \alpha} + +With the properties of this example, :math:`\sigma_{f} = 8868` Pa is reached after a cooling +of only :math:`\Delta T_{f} = -39.4` K, that is, less than half of the imposed temperature +drop. Beyond that point the rock deforms plastically and the stress stays on the yield +surface. + +**Lateral displacement.** Because the ``y`` direction is blocked while ``x`` and ``z`` are +free, all the deformation shows up laterally, and this gives a second, kinematic check that is +independent from the stress. In the elastic regime, + +.. math:: + + \varepsilon_{xx} = \alpha \, \Delta T \left( 1 + \frac{\lambda}{2(\lambda + G)} \right) + +Once the yield surface is reached, plastic flow adds lateral strain while the stress stays +put: the elasto-plastic column keeps contracting **faster** than the elastic one. Both +branches are implemented in ``AnalyticalSol.py``, which is the reference solution used below. + +------------------------------------------------------------------ +Running the case and post-processing +------------------------------------------------------------------ + +Both cases are run independently, each in its own directory: + +.. code-block:: console + + geosx -i ThermoElastic_1DCooling_fim_smoke.xml + geosx -i ThermoDruckerPrager_1DCooling_fim_smoke.xml + +Each run writes ``stressHistory.hdf5`` and ``displacementHistory.hdf5`` through the +``TimeHistory`` outputs. Those files are **not** stored in the repository; the curves shown +below are extracted once, locally, into a small CSV file with: + +.. code-block:: console + + python3 postprocess1DCooling.py -e -d + +The figure of this page is then generated at documentation build time from that CSV only. + +------------------------------------------------------------------ +Results +------------------------------------------------------------------ + +.. plot:: docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/plot1DCooling.py + +The figure on the left shows the stress induced by the confined cooling. Up to +:math:`-\Delta T \approx 39` K the two models are indistinguishable and follow the elastic +line :math:`-E \alpha \Delta T` exactly. Past that threshold, the elastic rock keeps +accumulating tensile stress and reaches 17.8 kPa at the end of the cooling, whereas the +Drucker-Prager rock **yields** and its stress saturates at the analytical cap +:math:`\sigma_{f}`, matching to machine precision. + +The figure in the middle explains the mechanism in the invariant plane. Because +:math:`\sigma_{xx} = \sigma_{zz} = 0`, the loading path is the straight line :math:`Q = 3P`, +whatever the amount of cooling. The elastic path simply crosses the Drucker-Prager envelope +and keeps going, which is physically inadmissible; the elasto-plastic path stops on the +envelope and slides along it. + +The figure on the right shows the kinematic counterpart. Up to the failure threshold, the two columns +contract identically. Beyond it, the roles reverse with respect to the stress plot: the +elasto-plastic column, whose stress is now frozen, contracts **more** than the elastic one, +reaching :math:`-32.7` against :math:`-29.7` µm. Plastic flow converts what would have been +additional stress into additional strain. GEOS matches the analytical displacement of both +branches to machine precision. + +The practical consequence is that **the safe amount of cooling is set by the strength of the +rock, not by its stiffness alone**: an elastic-only analysis of a cold injection would +over-predict the stress by a factor of two here, under-predict the deformation, and miss the +failure entirely. + +------------------------------------------------------------------ +To go further +------------------------------------------------------------------ + +**Feedback on this example** + +For any feedback on this example, please submit a `GitHub issue on the project's GitHub page `_. diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/cooling1D.csv b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/cooling1D.csv new file mode 100644 index 00000000000..57295a934f9 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/cooling1D.csv @@ -0,0 +1,102 @@ +time,temperature,dT,sigma_elastic_geos,sigma_druckerPrager_geos,sigma_elastic_analytical,sigma_yield,dp_b,dp_c,ux_elastic_geos,ux_druckerPrager_geos,ux_elastic_analytical,ux_druckerPrager_analytical +0.00000000e+00,1.00000000e+02,0.00000000e+00,0.00000000e+00,0.00000000e+00,-0.00000000e+00,8.86840336e+03,5.77427955e-01,1.05753580e+04,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00 +1.00000000e-10,1.00000000e+02,0.00000000e+00,0.00000000e+00,0.00000000e+00,-0.00000000e+00,8.86840336e+03,5.77427955e-01,1.05753580e+04,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00 +1.00000000e+00,9.92000000e+01,-8.00000000e-01,1.80000000e+02,1.80000000e+02,1.80000000e+02,8.86840336e+03,5.77427955e-01,1.05753580e+04,-3.00000000e-07,-3.00000000e-07,-3.00000000e-07,-3.00000000e-07 +2.00000000e+00,9.84000000e+01,-1.60000000e+00,3.60000000e+02,3.60000000e+02,3.60000000e+02,8.86840336e+03,5.77427955e-01,1.05753580e+04,-6.00000000e-07,-6.00000000e-07,-6.00000000e-07,-6.00000000e-07 +3.00000000e+00,9.76000000e+01,-2.40000000e+00,5.40000000e+02,5.40000000e+02,5.40000000e+02,8.86840336e+03,5.77427955e-01,1.05753580e+04,-9.00000000e-07,-9.00000000e-07,-9.00000000e-07,-9.00000000e-07 +4.00000000e+00,9.68000000e+01,-3.20000000e+00,7.20000000e+02,7.20000000e+02,7.20000000e+02,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.20000000e-06,-1.20000000e-06,-1.20000000e-06,-1.20000000e-06 +5.00000000e+00,9.60000000e+01,-4.00000000e+00,9.00000000e+02,9.00000000e+02,9.00000000e+02,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.50000000e-06,-1.50000000e-06,-1.50000000e-06,-1.50000000e-06 +6.00000000e+00,9.52000000e+01,-4.80000000e+00,1.08000000e+03,1.08000000e+03,1.08000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.80000000e-06,-1.80000000e-06,-1.80000000e-06,-1.80000000e-06 +7.00000000e+00,9.44000000e+01,-5.60000000e+00,1.26000000e+03,1.26000000e+03,1.26000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.10000000e-06,-2.10000000e-06,-2.10000000e-06,-2.10000000e-06 +8.00000000e+00,9.36000000e+01,-6.40000000e+00,1.44000000e+03,1.44000000e+03,1.44000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.40000000e-06,-2.40000000e-06,-2.40000000e-06,-2.40000000e-06 +9.00000000e+00,9.28000000e+01,-7.20000000e+00,1.62000000e+03,1.62000000e+03,1.62000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.70000000e-06,-2.70000000e-06,-2.70000000e-06,-2.70000000e-06 +1.00000000e+01,9.20000000e+01,-8.00000000e+00,1.80000000e+03,1.80000000e+03,1.80000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-3.00000000e-06,-3.00000000e-06,-3.00000000e-06,-3.00000000e-06 +1.10000000e+01,9.12000000e+01,-8.80000000e+00,1.98000000e+03,1.98000000e+03,1.98000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-3.30000000e-06,-3.30000000e-06,-3.30000000e-06,-3.30000000e-06 +1.20000000e+01,9.04000000e+01,-9.60000000e+00,2.16000000e+03,2.16000000e+03,2.16000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-3.60000000e-06,-3.60000000e-06,-3.60000000e-06,-3.60000000e-06 +1.30000000e+01,8.96000000e+01,-1.04000000e+01,2.34000000e+03,2.34000000e+03,2.34000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-3.90000000e-06,-3.90000000e-06,-3.90000000e-06,-3.90000000e-06 +1.40000000e+01,8.88000000e+01,-1.12000000e+01,2.52000000e+03,2.52000000e+03,2.52000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-4.20000000e-06,-4.20000000e-06,-4.20000000e-06,-4.20000000e-06 +1.50000000e+01,8.80000000e+01,-1.20000000e+01,2.70000000e+03,2.70000000e+03,2.70000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-4.50000000e-06,-4.50000000e-06,-4.50000000e-06,-4.50000000e-06 +1.60000000e+01,8.72000000e+01,-1.28000000e+01,2.88000000e+03,2.88000000e+03,2.88000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-4.80000000e-06,-4.80000000e-06,-4.80000000e-06,-4.80000000e-06 +1.70000000e+01,8.64000000e+01,-1.36000000e+01,3.06000000e+03,3.06000000e+03,3.06000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-5.10000000e-06,-5.10000000e-06,-5.10000000e-06,-5.10000000e-06 +1.80000000e+01,8.56000000e+01,-1.44000000e+01,3.24000000e+03,3.24000000e+03,3.24000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-5.40000000e-06,-5.40000000e-06,-5.40000000e-06,-5.40000000e-06 +1.90000000e+01,8.48000000e+01,-1.52000000e+01,3.42000000e+03,3.42000000e+03,3.42000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-5.70000000e-06,-5.70000000e-06,-5.70000000e-06,-5.70000000e-06 +2.00000000e+01,8.40000000e+01,-1.60000000e+01,3.60000000e+03,3.60000000e+03,3.60000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-6.00000000e-06,-6.00000000e-06,-6.00000000e-06,-6.00000000e-06 +2.10000000e+01,8.32000000e+01,-1.68000000e+01,3.78000000e+03,3.78000000e+03,3.78000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-6.30000000e-06,-6.30000000e-06,-6.30000000e-06,-6.30000000e-06 +2.20000000e+01,8.24000000e+01,-1.76000000e+01,3.96000000e+03,3.96000000e+03,3.96000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-6.60000000e-06,-6.60000000e-06,-6.60000000e-06,-6.60000000e-06 +2.30000000e+01,8.16000000e+01,-1.84000000e+01,4.14000000e+03,4.14000000e+03,4.14000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-6.90000000e-06,-6.90000000e-06,-6.90000000e-06,-6.90000000e-06 +2.40000000e+01,8.08000000e+01,-1.92000000e+01,4.32000000e+03,4.32000000e+03,4.32000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-7.20000000e-06,-7.20000000e-06,-7.20000000e-06,-7.20000000e-06 +2.50000000e+01,8.00000000e+01,-2.00000000e+01,4.50000000e+03,4.50000000e+03,4.50000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-7.50000000e-06,-7.50000000e-06,-7.50000000e-06,-7.50000000e-06 +2.60000000e+01,7.92000000e+01,-2.08000000e+01,4.68000000e+03,4.68000000e+03,4.68000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-7.80000000e-06,-7.80000000e-06,-7.80000000e-06,-7.80000000e-06 +2.70000000e+01,7.84000000e+01,-2.16000000e+01,4.86000000e+03,4.86000000e+03,4.86000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-8.10000000e-06,-8.10000000e-06,-8.10000000e-06,-8.10000000e-06 +2.80000000e+01,7.76000000e+01,-2.24000000e+01,5.04000000e+03,5.04000000e+03,5.04000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-8.40000000e-06,-8.40000000e-06,-8.40000000e-06,-8.40000000e-06 +2.90000000e+01,7.68000000e+01,-2.32000000e+01,5.22000000e+03,5.22000000e+03,5.22000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-8.70000000e-06,-8.70000000e-06,-8.70000000e-06,-8.70000000e-06 +3.00000000e+01,7.60000000e+01,-2.40000000e+01,5.40000000e+03,5.40000000e+03,5.40000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-9.00000000e-06,-9.00000000e-06,-9.00000000e-06,-9.00000000e-06 +3.10000000e+01,7.52000000e+01,-2.48000000e+01,5.58000000e+03,5.58000000e+03,5.58000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-9.30000000e-06,-9.30000000e-06,-9.30000000e-06,-9.30000000e-06 +3.20000000e+01,7.44000000e+01,-2.56000000e+01,5.76000000e+03,5.76000000e+03,5.76000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-9.60000000e-06,-9.60000000e-06,-9.60000000e-06,-9.60000000e-06 +3.30000000e+01,7.36000000e+01,-2.64000000e+01,5.94000000e+03,5.94000000e+03,5.94000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-9.90000000e-06,-9.90000000e-06,-9.90000000e-06,-9.90000000e-06 +3.40000000e+01,7.28000000e+01,-2.72000000e+01,6.12000000e+03,6.12000000e+03,6.12000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.02000000e-05,-1.02000000e-05,-1.02000000e-05,-1.02000000e-05 +3.50000000e+01,7.20000000e+01,-2.80000000e+01,6.30000000e+03,6.30000000e+03,6.30000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.05000000e-05,-1.05000000e-05,-1.05000000e-05,-1.05000000e-05 +3.60000000e+01,7.12000000e+01,-2.88000000e+01,6.48000000e+03,6.48000000e+03,6.48000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.08000000e-05,-1.08000000e-05,-1.08000000e-05,-1.08000000e-05 +3.70000000e+01,7.04000000e+01,-2.96000000e+01,6.66000000e+03,6.66000000e+03,6.66000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.11000000e-05,-1.11000000e-05,-1.11000000e-05,-1.11000000e-05 +3.80000000e+01,6.96000000e+01,-3.04000000e+01,6.84000000e+03,6.84000000e+03,6.84000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.14000000e-05,-1.14000000e-05,-1.14000000e-05,-1.14000000e-05 +3.90000000e+01,6.88000000e+01,-3.12000000e+01,7.02000000e+03,7.02000000e+03,7.02000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.17000000e-05,-1.17000000e-05,-1.17000000e-05,-1.17000000e-05 +4.00000000e+01,6.80000000e+01,-3.20000000e+01,7.20000000e+03,7.20000000e+03,7.20000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.20000000e-05,-1.20000000e-05,-1.20000000e-05,-1.20000000e-05 +4.10000000e+01,6.72000000e+01,-3.28000000e+01,7.38000000e+03,7.38000000e+03,7.38000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.23000000e-05,-1.23000000e-05,-1.23000000e-05,-1.23000000e-05 +4.20000000e+01,6.64000000e+01,-3.36000000e+01,7.56000000e+03,7.56000000e+03,7.56000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.26000000e-05,-1.26000000e-05,-1.26000000e-05,-1.26000000e-05 +4.30000000e+01,6.56000000e+01,-3.44000000e+01,7.74000000e+03,7.74000000e+03,7.74000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.29000000e-05,-1.29000000e-05,-1.29000000e-05,-1.29000000e-05 +4.40000000e+01,6.48000000e+01,-3.52000000e+01,7.92000000e+03,7.92000000e+03,7.92000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.32000000e-05,-1.32000000e-05,-1.32000000e-05,-1.32000000e-05 +4.50000000e+01,6.40000000e+01,-3.60000000e+01,8.10000000e+03,8.10000000e+03,8.10000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.35000000e-05,-1.35000000e-05,-1.35000000e-05,-1.35000000e-05 +4.60000000e+01,6.32000000e+01,-3.68000000e+01,8.28000000e+03,8.28000000e+03,8.28000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.38000000e-05,-1.38000000e-05,-1.38000000e-05,-1.38000000e-05 +4.70000000e+01,6.24000000e+01,-3.76000000e+01,8.46000000e+03,8.46000000e+03,8.46000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.41000000e-05,-1.41000000e-05,-1.41000000e-05,-1.41000000e-05 +4.80000000e+01,6.16000000e+01,-3.84000000e+01,8.64000000e+03,8.64000000e+03,8.64000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.44000000e-05,-1.44000000e-05,-1.44000000e-05,-1.44000000e-05 +4.90000000e+01,6.08000000e+01,-3.92000000e+01,8.82000000e+03,8.82000000e+03,8.82000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.47000000e-05,-1.47000000e-05,-1.47000000e-05,-1.47000000e-05 +5.00000000e+01,6.00000000e+01,-4.00000000e+01,9.00000000e+03,8.86840336e+03,9.00000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.50000000e-05,-1.50438655e-05,-1.50000000e-05,-1.50438655e-05 +5.10000000e+01,5.92000000e+01,-4.08000000e+01,9.18000000e+03,8.86840336e+03,9.18000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.53000000e-05,-1.54038655e-05,-1.53000000e-05,-1.54038655e-05 +5.20000000e+01,5.84000000e+01,-4.16000000e+01,9.36000000e+03,8.86840336e+03,9.36000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.56000000e-05,-1.57638655e-05,-1.56000000e-05,-1.57638655e-05 +5.30000000e+01,5.76000000e+01,-4.24000000e+01,9.54000000e+03,8.86840336e+03,9.54000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.59000000e-05,-1.61238655e-05,-1.59000000e-05,-1.61238655e-05 +5.40000000e+01,5.68000000e+01,-4.32000000e+01,9.72000000e+03,8.86840336e+03,9.72000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.62000000e-05,-1.64838655e-05,-1.62000000e-05,-1.64838655e-05 +5.50000000e+01,5.60000000e+01,-4.40000000e+01,9.90000000e+03,8.86840336e+03,9.90000000e+03,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.65000000e-05,-1.68438655e-05,-1.65000000e-05,-1.68438655e-05 +5.60000000e+01,5.52000000e+01,-4.48000000e+01,1.00800000e+04,8.86840336e+03,1.00800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.68000000e-05,-1.72038655e-05,-1.68000000e-05,-1.72038655e-05 +5.70000000e+01,5.44000000e+01,-4.56000000e+01,1.02600000e+04,8.86840336e+03,1.02600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.71000000e-05,-1.75638655e-05,-1.71000000e-05,-1.75638655e-05 +5.80000000e+01,5.36000000e+01,-4.64000000e+01,1.04400000e+04,8.86840336e+03,1.04400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.74000000e-05,-1.79238655e-05,-1.74000000e-05,-1.79238655e-05 +5.90000000e+01,5.28000000e+01,-4.72000000e+01,1.06200000e+04,8.86840336e+03,1.06200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.77000000e-05,-1.82838655e-05,-1.77000000e-05,-1.82838655e-05 +6.00000000e+01,5.20000000e+01,-4.80000000e+01,1.08000000e+04,8.86840336e+03,1.08000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.80000000e-05,-1.86438655e-05,-1.80000000e-05,-1.86438655e-05 +6.10000000e+01,5.12000000e+01,-4.88000000e+01,1.09800000e+04,8.86840336e+03,1.09800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.83000000e-05,-1.90038655e-05,-1.83000000e-05,-1.90038655e-05 +6.20000000e+01,5.04000000e+01,-4.96000000e+01,1.11600000e+04,8.86840336e+03,1.11600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.86000000e-05,-1.93638655e-05,-1.86000000e-05,-1.93638655e-05 +6.30000000e+01,4.96000000e+01,-5.04000000e+01,1.13400000e+04,8.86840336e+03,1.13400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.89000000e-05,-1.97238655e-05,-1.89000000e-05,-1.97238655e-05 +6.40000000e+01,4.88000000e+01,-5.12000000e+01,1.15200000e+04,8.86840336e+03,1.15200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.92000000e-05,-2.00838655e-05,-1.92000000e-05,-2.00838655e-05 +6.50000000e+01,4.80000000e+01,-5.20000000e+01,1.17000000e+04,8.86840336e+03,1.17000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.95000000e-05,-2.04438655e-05,-1.95000000e-05,-2.04438655e-05 +6.60000000e+01,4.72000000e+01,-5.28000000e+01,1.18800000e+04,8.86840336e+03,1.18800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-1.98000000e-05,-2.08038655e-05,-1.98000000e-05,-2.08038655e-05 +6.70000000e+01,4.64000000e+01,-5.36000000e+01,1.20600000e+04,8.86840336e+03,1.20600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.01000000e-05,-2.11638655e-05,-2.01000000e-05,-2.11638655e-05 +6.80000000e+01,4.56000000e+01,-5.44000000e+01,1.22400000e+04,8.86840336e+03,1.22400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.04000000e-05,-2.15238655e-05,-2.04000000e-05,-2.15238655e-05 +6.90000000e+01,4.48000000e+01,-5.52000000e+01,1.24200000e+04,8.86840336e+03,1.24200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.07000000e-05,-2.18838655e-05,-2.07000000e-05,-2.18838655e-05 +7.00000000e+01,4.40000000e+01,-5.60000000e+01,1.26000000e+04,8.86840336e+03,1.26000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.10000000e-05,-2.22438655e-05,-2.10000000e-05,-2.22438655e-05 +7.10000000e+01,4.32000000e+01,-5.68000000e+01,1.27800000e+04,8.86840336e+03,1.27800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.13000000e-05,-2.26038655e-05,-2.13000000e-05,-2.26038655e-05 +7.20000000e+01,4.24000000e+01,-5.76000000e+01,1.29600000e+04,8.86840336e+03,1.29600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.16000000e-05,-2.29638655e-05,-2.16000000e-05,-2.29638655e-05 +7.30000000e+01,4.16000000e+01,-5.84000000e+01,1.31400000e+04,8.86840336e+03,1.31400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.19000000e-05,-2.33238655e-05,-2.19000000e-05,-2.33238655e-05 +7.40000000e+01,4.08000000e+01,-5.92000000e+01,1.33200000e+04,8.86840336e+03,1.33200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.22000000e-05,-2.36838655e-05,-2.22000000e-05,-2.36838655e-05 +7.50000000e+01,4.00000000e+01,-6.00000000e+01,1.35000000e+04,8.86840336e+03,1.35000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.25000000e-05,-2.40438655e-05,-2.25000000e-05,-2.40438655e-05 +7.60000000e+01,3.92000000e+01,-6.08000000e+01,1.36800000e+04,8.86840336e+03,1.36800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.28000000e-05,-2.44038655e-05,-2.28000000e-05,-2.44038655e-05 +7.70000000e+01,3.84000000e+01,-6.16000000e+01,1.38600000e+04,8.86840336e+03,1.38600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.31000000e-05,-2.47638655e-05,-2.31000000e-05,-2.47638655e-05 +7.80000000e+01,3.76000000e+01,-6.24000000e+01,1.40400000e+04,8.86840336e+03,1.40400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.34000000e-05,-2.51238655e-05,-2.34000000e-05,-2.51238655e-05 +7.90000000e+01,3.68000000e+01,-6.32000000e+01,1.42200000e+04,8.86840336e+03,1.42200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.37000000e-05,-2.54838655e-05,-2.37000000e-05,-2.54838655e-05 +8.00000000e+01,3.60000000e+01,-6.40000000e+01,1.44000000e+04,8.86840336e+03,1.44000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.40000000e-05,-2.58438655e-05,-2.40000000e-05,-2.58438655e-05 +8.10000000e+01,3.52000000e+01,-6.48000000e+01,1.45800000e+04,8.86840336e+03,1.45800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.43000000e-05,-2.62038655e-05,-2.43000000e-05,-2.62038655e-05 +8.20000000e+01,3.44000000e+01,-6.56000000e+01,1.47600000e+04,8.86840336e+03,1.47600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.46000000e-05,-2.65638655e-05,-2.46000000e-05,-2.65638655e-05 +8.30000000e+01,3.36000000e+01,-6.64000000e+01,1.49400000e+04,8.86840336e+03,1.49400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.49000000e-05,-2.69238655e-05,-2.49000000e-05,-2.69238655e-05 +8.40000000e+01,3.28000000e+01,-6.72000000e+01,1.51200000e+04,8.86840336e+03,1.51200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.52000000e-05,-2.72838655e-05,-2.52000000e-05,-2.72838655e-05 +8.50000000e+01,3.20000000e+01,-6.80000000e+01,1.53000000e+04,8.86840336e+03,1.53000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.55000000e-05,-2.76438655e-05,-2.55000000e-05,-2.76438655e-05 +8.60000000e+01,3.12000000e+01,-6.88000000e+01,1.54800000e+04,8.86840336e+03,1.54800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.58000000e-05,-2.80038655e-05,-2.58000000e-05,-2.80038655e-05 +8.70000000e+01,3.04000000e+01,-6.96000000e+01,1.56600000e+04,8.86840336e+03,1.56600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.61000000e-05,-2.83638655e-05,-2.61000000e-05,-2.83638655e-05 +8.80000000e+01,2.96000000e+01,-7.04000000e+01,1.58400000e+04,8.86840336e+03,1.58400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.64000000e-05,-2.87238655e-05,-2.64000000e-05,-2.87238655e-05 +8.90000000e+01,2.88000000e+01,-7.12000000e+01,1.60200000e+04,8.86840336e+03,1.60200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.67000000e-05,-2.90838655e-05,-2.67000000e-05,-2.90838655e-05 +9.00000000e+01,2.80000000e+01,-7.20000000e+01,1.62000000e+04,8.86840336e+03,1.62000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.70000000e-05,-2.94438655e-05,-2.70000000e-05,-2.94438655e-05 +9.10000000e+01,2.72000000e+01,-7.28000000e+01,1.63800000e+04,8.86840336e+03,1.63800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.73000000e-05,-2.98038655e-05,-2.73000000e-05,-2.98038655e-05 +9.20000000e+01,2.64000000e+01,-7.36000000e+01,1.65600000e+04,8.86840336e+03,1.65600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.76000000e-05,-3.01638655e-05,-2.76000000e-05,-3.01638655e-05 +9.30000000e+01,2.56000000e+01,-7.44000000e+01,1.67400000e+04,8.86840336e+03,1.67400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.79000000e-05,-3.05238655e-05,-2.79000000e-05,-3.05238655e-05 +9.40000000e+01,2.48000000e+01,-7.52000000e+01,1.69200000e+04,8.86840336e+03,1.69200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.82000000e-05,-3.08838655e-05,-2.82000000e-05,-3.08838655e-05 +9.50000000e+01,2.40000000e+01,-7.60000000e+01,1.71000000e+04,8.86840336e+03,1.71000000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.85000000e-05,-3.12438655e-05,-2.85000000e-05,-3.12438655e-05 +9.60000000e+01,2.32000000e+01,-7.68000000e+01,1.72800000e+04,8.86840336e+03,1.72800000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.88000000e-05,-3.16038655e-05,-2.88000000e-05,-3.16038655e-05 +9.70000000e+01,2.24000000e+01,-7.76000000e+01,1.74600000e+04,8.86840336e+03,1.74600000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.91000000e-05,-3.19638655e-05,-2.91000000e-05,-3.19638655e-05 +9.80000000e+01,2.16000000e+01,-7.84000000e+01,1.76400000e+04,8.86840336e+03,1.76400000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.94000000e-05,-3.23238655e-05,-2.94000000e-05,-3.23238655e-05 +9.90000000e+01,2.08000000e+01,-7.92000000e+01,1.78200000e+04,8.86840336e+03,1.78200000e+04,8.86840336e+03,5.77427955e-01,1.05753580e+04,-2.97000000e-05,-3.26838655e-05,-2.97000000e-05,-3.26838655e-05 diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/fixed_1d_bar_cooling_3d_view.png b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/fixed_1d_bar_cooling_3d_view.png new file mode 100755 index 00000000000..bea20282775 Binary files /dev/null and b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/fixed_1d_bar_cooling_3d_view.png differ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/plot1DCooling.py b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/plot1DCooling.py new file mode 100644 index 00000000000..d9f95d0fd3e --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/plot1DCooling.py @@ -0,0 +1,127 @@ +"""Figure of the 1D tutorial: effect of cooling on failure. + +Executed by Sphinx (the `.. plot::` directive of Example.rst) on every documentation build. +It reads a single versioned CSV file sitting next to it: + + cooling1D.csv + +That CSV is produced offline by postprocess1DCooling.py from the TimeHistory (.hdf5) outputs +of both runs. The .hdf5 files are not versioned; this script therefore depends only on numpy +and matplotlib. + +To regenerate the data after a change in the cases: + geosx -i ThermoElastic_1DCooling_fim_smoke.xml # in elastic/ + geosx -i ThermoDruckerPrager_1DCooling_fim_smoke.xml # in druckerPrager/ + python3 postprocess1DCooling.py -e elastic -d druckerPrager -o +""" +import os + +import numpy as np +import matplotlib.pyplot as plt + +try: + HERE = os.path.dirname(os.path.abspath(__file__)) +except NameError: + # the Sphinx `.. plot::` directive may execute this file without defining __file__ + HERE = os.getcwd() + +fsize = 20 +msize = 7 +lw = 3 +mew = 1.8 +malpha = 0.6 +lalpha = 0.8 +N1 = 2 # one GEOS marker out of two, to keep the plot readable + + +def main(): + path = os.path.join(HERE, "cooling1D.csv") + if not os.path.exists(path): + raise FileNotFoundError( + f"cooling1D.csv not found in {HERE}. Regenerate it with postprocess1DCooling.py " + f"(see the header of this file).") + d = np.genfromtxt(path, delimiter=",", names=True) + + dT = d["dT"] + sy = float(d["sigma_yield"][0]) + b = float(d["dp_b"][0]) + c = float(d["dp_c"][0]) + # cooling at which the yield surface is reached + dT_yield = float(np.interp(sy, d["sigma_elastic_analytical"], dT)) + + cmap = plt.get_cmap("tab10") + fig, ax = plt.subplots(1, 3, figsize=(25, 7)) + + # ---------------------------------------------------------------- panel 1 + # stress induced by the confined cooling, as a function of -dT + a = ax[0] + x = -dT + a.plot(x, d["sigma_elastic_analytical"] / 1e3, '-', color=cmap(0), lw=lw, alpha=lalpha, + label=r'Analytical elastic: $\sigma_{yy} = -E\,\alpha\,\Delta T$') + a.plot(x[::N1], d["sigma_elastic_geos"][::N1] / 1e3, 'o', color=cmap(0), markersize=msize, + fillstyle='none', mew=mew, alpha=malpha, label='GEOS: thermo-elastic') + a.plot(x[::N1], d["sigma_druckerPrager_geos"][::N1] / 1e3, 's', color=cmap(1), + markersize=msize, fillstyle='none', mew=mew, alpha=malpha, + label='GEOS: thermo-plastic (Drucker-Prager)') + a.axhline(sy / 1e3, color=cmap(1), ls='--', lw=lw * 0.7, alpha=0.9, + label=r'Failure cap: $\sigma_{f} = c\,/\,(1 + b/3)$') + a.axvline(-dT_yield, color='0.4', ls=':', lw=lw * 0.7) + a.annotate(f'failure at $\\Delta T$ = {dT_yield:.1f} K', + xy=(-dT_yield, sy / 1e3), xytext=(-dT_yield + 4, sy / 1e3 * 0.45), + fontsize=fsize * 0.7, color='0.25', + arrowprops=dict(arrowstyle='->', color='0.45', lw=1.4)) + a.set_xlabel(r'Cooling $-\Delta T$ [K]', size=fsize, weight='bold') + a.set_ylabel(r'Total stress $\sigma_{yy}$ [kPa]', size=fsize, weight='bold') + a.set_title('Thermally induced stress under confined cooling', size=fsize * 0.85) + a.grid(True, alpha=0.35) + a.legend(loc='upper left', fontsize=fsize * 0.62) + a.tick_params(labelsize=fsize * 0.8) + + # ---------------------------------------------------------------- panel 2 + # loading path in the (P, Q) plane, together with the Drucker-Prager envelope + a = ax[1] + Pe, Qe = d["sigma_elastic_geos"] / 3.0, d["sigma_elastic_geos"] + Pd, Qd = d["sigma_druckerPrager_geos"] / 3.0, d["sigma_druckerPrager_geos"] + Pmax = max(Pe.max(), c / b) * 1.05 + Penv = np.linspace(0.0, c / b, 200) + a.plot(Penv / 1e3, (c - b * Penv) / 1e3, 'k--', lw=lw * 0.8, alpha=lalpha, + label=r'Drucker-Prager envelope: $Q = c - b\,P$') + a.plot(Pe / 1e3, Qe / 1e3, '-', color=cmap(0), lw=lw, alpha=lalpha, + label='thermo-elastic path') + a.plot(Pd[::N1] / 1e3, Qd[::N1] / 1e3, 's', color=cmap(1), markersize=msize, + fillstyle='none', mew=mew, alpha=malpha, label='thermo-plastic path') + a.plot([sy / 3.0 / 1e3], [sy / 1e3], '*', color='k', markersize=msize * 2.4, + label='failure point') + a.set_xlim(0.0, Pmax / 1e3) + a.set_xlabel(r'Mean stress $P = \mathrm{tr}(\sigma)/3$ [kPa]', size=fsize, weight='bold') + a.set_ylabel(r'Von Mises stress $Q$ [kPa]', size=fsize, weight='bold') + a.set_title('Stress path in the invariant plane', size=fsize * 0.85) + a.grid(True, alpha=0.35) + a.legend(loc='upper right', fontsize=fsize * 0.62) + a.tick_params(labelsize=fsize * 0.8) + + # ---------------------------------------------------------------- panel 3 + # lateral displacement of the free face: yielding adds to it + a = ax[2] + a.plot(x, d["ux_elastic_analytical"] * 1e6, '-', color=cmap(0), lw=lw, alpha=lalpha, + label='Analytical: elastic') + a.plot(x[::N1], d["ux_elastic_geos"][::N1] * 1e6, 'o', color=cmap(0), markersize=msize, + fillstyle='none', mew=mew, alpha=malpha, label='GEOS: thermo-elastic') + a.plot(x, d["ux_druckerPrager_analytical"] * 1e6, '-', color=cmap(1), lw=lw, alpha=lalpha, + label='Analytical: Drucker-Prager') + a.plot(x[::N1], d["ux_druckerPrager_geos"][::N1] * 1e6, 's', color=cmap(1), markersize=msize, + fillstyle='none', mew=mew, alpha=malpha, label='GEOS: thermo-plastic') + a.axvline(-dT_yield, color='0.4', ls=':', lw=lw * 0.7) + a.set_xlabel(r'Cooling $-\Delta T$ [K]', size=fsize, weight='bold') + a.set_ylabel(r'Lateral displacement $u_x$ [$\mu$m]', size=fsize, weight='bold') + a.set_title('Lateral displacement of the free face', size=fsize * 0.85) + a.grid(True, alpha=0.35) + a.legend(loc='lower left', fontsize=fsize * 0.62) + a.tick_params(labelsize=fsize * 0.8) + + fig.tight_layout() + plt.show() + + +if __name__ == "__main__": + main() diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/postprocess1DCooling.py b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/postprocess1DCooling.py new file mode 100644 index 00000000000..c5cf4c6b963 --- /dev/null +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/postprocess1DCooling.py @@ -0,0 +1,192 @@ +"""LOCAL post-processing of the 1D confined cooling case (thermo-elastic vs Drucker-Prager). + +Run by hand, after running both cases: + + geosx -i ThermoElastic_1DCooling_fim_smoke.xml # in an elastic/ directory + geosx -i ThermoDruckerPrager_1DCooling_fim_smoke.xml # in a druckerPrager/ directory + python3 postprocess1DCooling.py -e elastic -d druckerPrager -o + +It reads the GEOS TimeHistory outputs (stressHistory.hdf5, displacementHistory.hdf5) and +writes a single lightweight CSV: + + cooling1D.csv temperature, GEOS stress and displacement of both models, analytical solution + +That CSV is the file kept under version control. The .hdf5 files are not: the documentation +figure is produced by plot1DCooling.py, which reads this CSV only and depends solely on numpy +and matplotlib. + +Analytical solution +------------------- +The column is fixed along y (both yneg and ypos are constrained) and free along x and z, hence +a state of uniaxial strain: + eps_yy = 0, sigma_xx = sigma_zz = 0 +The thermo-elastic law then gives a stress that is purely induced by the cooling: + sigma_yy = -E * alpha * dT (tensile > 0 when dT < 0) + +The Drucker-Prager criterion used by GEOS reads F = Q + b*P - c (DruckerPrager.hpp), with +P = tr(sigma)/3, Q the von Mises stress, and (DruckerPrager.cpp): + b = 6 sin(phi) / (3 - sin(phi)) + c = 6 * cohesion * cos(phi) / (3 - sin(phi)) +For this uniaxial state, Q = sigma_yy and P = sigma_yy/3, which yields the failure cap: + sigma_failure = c / (1 + b/3) +and the cooling required to reach it: dT_failure = -sigma_failure / (E * alpha). +""" +import argparse +import math +import os +import sys + +import numpy as np +import h5py + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from AnalyticalSol import DPAnalyticalSolution # noqa: E402 + + +# --- case properties (ThermoElastic/ThermoDruckerPrager_1DCooling_fim_smoke.xml) --- +BULK_MODULUS = 0.5e9 # defaultBulkModulus +SHEAR_MODULUS = 0.3e9 # defaultShearModulus +ALPHA = 3e-7 # defaultDrainedLinearTEC +COHESION = 5e3 # defaultCohesion +FRICTION_ANGLE = 15.27 # defaultFrictionAngle (degrees) +T_INITIAL = 100.0 # initialTemperature +T_FINAL = 20.0 # final value of timeFunction +T_RAMP_END = 100.0 # maxTime: end of the cooling ramp +COLUMN_WIDTH = 1.0 # extent along x (xCoords = { 0, 1 }) + + +def youngModulus(): + K, G = BULK_MODULUS, SHEAR_MODULUS + return 9.0 * K * G / (3.0 * K + G) + + +def druckerPragerCap(): + """(sigma_failure, dT_failure, b, c) for the uniaxial stress state of this case.""" + phi = math.radians(FRICTION_ANGLE) + b = 6.0 * math.sin(phi) / (3.0 - math.sin(phi)) + c = 6.0 * COHESION * math.cos(phi) / (3.0 - math.sin(phi)) + sigma = c / (1.0 + b / 3.0) + dT = -sigma / (youngModulus() * ALPHA) + return sigma, dT, b, c + + +def imposedTemperature(t): + """Ramp imposed by the TableFunction of the XML: 100 up to 1e-10, then linear down to 20.""" + return np.interp(t, [0.0, 1e-10, T_RAMP_END], [T_INITIAL, T_INITIAL, T_FINAL]) + + +def elasticLateralDisplacement(dT): + """Lateral displacement u_x of the free face, purely elastic branch. + + Same formulation as AnalyticalSol.DPAnalyticalSolution, without the plasticity check: + eps_y^mech = -alpha*dT, eps_x^mech = -eps_y^mech * lambda / (2 (lambda + G)) + eps_x = eps_x^mech + alpha*dT + The column is 1 m wide along x, hence u_x = eps_x * 1. + """ + lam = BULK_MODULUS - 2.0 * SHEAR_MODULUS / 3.0 + epsMech_y = -ALPHA * dT + epsMech_x = -epsMech_y * lam / (2.0 * (lam + SHEAR_MODULUS)) + return (epsMech_x + ALPHA * dT) * COLUMN_WIDTH + + +def _validLength(t): + """The HDF5 buffer may be pre-allocated: cut at the first time value falling back to zero.""" + for j in range(1, len(t)): + if t[j] < 1e-12: + return j + return len(t) + + +def readCase(case_dir): + """(time, sigma_yy, u_x): stress at the centre of the column, lateral displacement of the + free face x = 1. Both fields are uniform in this 1D setup, which is checked here.""" + ps = os.path.join(case_dir, "stressHistory.hdf5") + pd = os.path.join(case_dir, "displacementHistory.hdf5") + for p in (ps, pd): + if not os.path.exists(p): + raise FileNotFoundError( + f"{p} not found. Run the corresponding case, then run this script again.") + + hs = h5py.File(ps, 'r') + t = np.array(hs.get('averageStress Time')).ravel() + S = np.array(hs.get('averageStress')) + n = _validLength(t) + t, S = t[:n], S[:n] + spread = float(np.abs(S[-1, :, 1] - S[-1, :, 1].mean()).max()) + if spread > 1e-6 * max(1.0, abs(S[-1, :, 1].mean())): + print(f" [warning] {case_dir}: sigma_yy is not uniform (spread {spread:.3e} Pa)") + sigma = S[:, S.shape[1] // 2, 1] + + hd = h5py.File(pd, 'r') + td = np.array(hd.get('totalDisplacement Time')).ravel()[:n] + U = np.array(hd.get('totalDisplacement'))[:n] + X = np.array(hd.get('totalDisplacement ReferencePosition'))[0] + if not np.allclose(t, td): + raise ValueError(f"{case_dir}: stress and displacement time bases differ") + free = np.where(np.abs(X[:, 0] - X[:, 0].max()) < 1e-9)[0] # free face x = 1 + ux = U[:, free, 0].mean(axis=1) + uy = float(np.abs(U[:, :, 1]).max()) + if uy > 1e-12: + print(f" [warning] {case_dir}: u_y is not zero ({uy:.3e} m), the column should be fixed") + return t, sigma, ux + + +HEADER = ("time,temperature,dT,sigma_elastic_geos,sigma_druckerPrager_geos," + "sigma_elastic_analytical,sigma_yield,dp_b,dp_c," + "ux_elastic_geos,ux_druckerPrager_geos," + "ux_elastic_analytical,ux_druckerPrager_analytical") + + +def main(): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("-e", "--elastic", default="elastic", help="directory of the thermo-elastic run") + ap.add_argument("-d", "--druckerPrager", default="druckerPrager", help="directory of the Drucker-Prager run") + ap.add_argument("-o", "--outdir", default=".", help="directory where the CSV is written") + args = ap.parse_args() + os.makedirs(args.outdir, exist_ok=True) + + te, se, uxe = readCase(args.elastic) + td, sd, uxd = readCase(args.druckerPrager) + if not np.allclose(te, td): + raise ValueError("both runs do not share the same time base: " + "check that the Events of the two XML files are identical") + + T = imposedTemperature(te) + dT = T - T_INITIAL + E = youngModulus() + sigma_ana = -E * ALPHA * dT # thermo-elastic, confined deformation + sigma_yield, dT_yield, b, c = druckerPragerCap() + + # lateral displacement: elastic branch in closed form above, elasto-plastic branch from + # AnalyticalSol.DPAnalyticalSolution (reference solution shipped with the PR) + ux_ana_el = elasticLateralDisplacement(dT) + dp = DPAnalyticalSolution() + ux_ana_dp = np.array([dp.compute_disp(x)[0] * COLUMN_WIDTH for x in dT]) + + # b and c are constant: stored as columns so that the plotting script can draw the failure + # envelope without having to re-import the material properties + data = np.column_stack([te, T, dT, se, sd, sigma_ana, + np.full_like(te, sigma_yield), + np.full_like(te, b), np.full_like(te, c), + uxe, uxd, ux_ana_el, ux_ana_dp]) + out = os.path.join(args.outdir, "cooling1D.csv") + np.savetxt(out, data, delimiter=",", header=HEADER, comments="", fmt="%.8e") + + print(f"{out} ({data.shape[0]} time steps)") + print(f" E = {E:.4g} Pa, alpha = {ALPHA:.1e} 1/K") + print(f" Drucker-Prager: b = {b:.6f}, c = {c:.2f} Pa") + print(f" sigma_failure = {sigma_yield:.2f} Pa reached at dT = {dT_yield:.2f} K " + f"(T = {T_INITIAL + dT_yield:.2f})") + err_el = np.abs(se - sigma_ana).max() + print(f" max deviation GEOS elastic / analytical: {err_el:.3e} Pa " + f"({100 * err_el / max(abs(sigma_ana).max(), 1e-30):.3f} %)") + print(f" Drucker-Prager plateau, GEOS: {sd[-1]:.2f} Pa " + f"(analytical {sigma_yield:.2f} Pa)") + print(f" final u_x, elastic : GEOS {uxe[-1]:+.6e} m / analytical {ux_ana_el[-1]:+.6e} m " + f"(deviation {abs(uxe[-1] - ux_ana_el[-1]):.2e})") + print(f" final u_x, Drucker-Prager: GEOS {uxd[-1]:+.6e} m / analytical {ux_ana_dp[-1]:+.6e} m " + f"(deviation {abs(uxd[-1] - ux_ana_dp[-1]):.2e})") + + +if __name__ == "__main__": + main() diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/xz_cross_section_uniaxial_stress.png b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/xz_cross_section_uniaxial_stress.png new file mode 100755 index 00000000000..d74f9da5e14 Binary files /dev/null and b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/1DCooling/xz_cross_section_uniaxial_stress.png differ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/Index.rst b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/Index.rst index 03e6d1dda41..12927df634d 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/Index.rst +++ b/src/docs/sphinx/advancedExamples/validationStudies/thermoPoromechanics/Index.rst @@ -9,4 +9,6 @@ Thermoporomechanics thermalConsolidation/Example + 1DCooling/Example + diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedContactThermoElasticWellbore/thermoElastic_casedContactWellbore_stress.py b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedContactThermoElasticWellbore/thermoElastic_casedContactWellbore_stress.py index bcac4ebd57d..c44c02f0772 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedContactThermoElasticWellbore/thermoElastic_casedContactWellbore_stress.py +++ b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedContactThermoElasticWellbore/thermoElastic_casedContactWellbore_stress.py @@ -39,9 +39,9 @@ def main(): temperature = np.asarray( hf_temperature.get('temperature') ) # Compute total stress - stress_xx_total = stress[:,:,0] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) - stress_yy_total = stress[:,:,1] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) - stress_zz_total = stress[:,:,2] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) + stress_xx_total = stress[:,:,0] + stress_yy_total = stress[:,:,1] + stress_zz_total = stress[:,:,2] stress_yz_total = stress[:,:,3] stress_xz_total = stress[:,:,4] stress_xy_total = stress[:,:,5] diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedThermoElasticWellbore/thermoElastic_casedWellbore_stress.py b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedThermoElasticWellbore/thermoElastic_casedWellbore_stress.py index 0ece31db68c..74874cfbc86 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedThermoElasticWellbore/thermoElastic_casedWellbore_stress.py +++ b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/casedThermoElasticWellbore/thermoElastic_casedWellbore_stress.py @@ -41,9 +41,9 @@ def main(): # Compute total stress: # With the actual version of GEOS, the output stress need to be combined with the temperature contribution to obtain the total tress as follows: - stress_xx_total = stress[:,:,0] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) - stress_yy_total = stress[:,:,1] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) - stress_zz_total = stress[:,:,2] - 3 * bulkModuli[idx] * thermalExpansionCoefficients[idx] * (temperature-100) + stress_xx_total = stress[:,:,0] + stress_yy_total = stress[:,:,1] + stress_zz_total = stress[:,:,2] stress_yz_total = stress[:,:,3] stress_xz_total = stress[:,:,4] stress_xy_total = stress[:,:,5] diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/THM_wellbore_stress.py b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/THM_wellbore_stress.py index 4448fa04bc8..e349cf99dc9 100644 --- a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/THM_wellbore_stress.py +++ b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/THM_wellbore_stress.py @@ -53,9 +53,9 @@ def main(): pressure = np.asarray( hf_pressure.get('pressure') ) # Compute total stress - stress_xx_total = stress[:,:,0] - BiotCoefficient * pressure - 3 * bulkModulus * thermalExpansionCoefficients * temperature - stress_yy_total = stress[:,:,1] - BiotCoefficient * pressure - 3 * bulkModulus * thermalExpansionCoefficients * temperature - stress_zz_total = stress[:,:,2] - BiotCoefficient * pressure - 3 * bulkModulus * thermalExpansionCoefficients * temperature + stress_xx_total = stress[:,:,0] - BiotCoefficient * pressure + stress_yy_total = stress[:,:,1] - BiotCoefficient * pressure + stress_zz_total = stress[:,:,2] - BiotCoefficient * pressure stress_yz_total = stress[:,:,3] stress_xz_total = stress[:,:,4] stress_xy_total = stress[:,:,5] diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/displacement.png b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/displacement.png index 854f7babfcb..6bb6f5f7361 100644 Binary files a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/displacement.png and b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/displacement.png differ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/stress.png b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/stress.png index c8659a07528..7f1487e9112 100644 Binary files a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/stress.png and b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/stress.png differ diff --git a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/temperature_pressure.png b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/temperature_pressure.png index cb086e3e0b7..317642fbf12 100644 Binary files a/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/temperature_pressure.png and b/src/docs/sphinx/advancedExamples/validationStudies/wellboreProblems/thermoPoroElasticWellbore/temperature_pressure.png differ