Merge pull request #2179 from wenqing/benchmark

Another test of BGRa creep model

ProcessLib/ThermoMechanics/Tests.cmake

@@ -102,3 +102,17 @@ AddTest(
     DIFF_DATA
     SimpleAxisymmetricCreepWithAnalyticSolution.vtu SimpleAxisymmetricCreepWithAnalyticalSolution_pcs_0_ts_1000_t_100.000000.vtu   analytic_strain_yy  epsilon_yy 1e-10 1e-4
 )
+
+AddTest(
+    NAME ThermoMechanics_CreepAfterExcavation
+    PATH ThermoMechanics/CreepBGRa/CreepAfterExcavation
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS CreepAfterExcavation.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu sigma_xx sigma_xx 1e-16 1e-11
+    ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu sigma_yy sigma_yy 1e-16 1e-11
+    ExpectedCreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu CreepAfterExcavation_pcs_0_ts_61_t_4320000.000000.vtu displacement displacement 1e-16 1e-9
+)

Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>mesh.vtu</mesh>
+    <geometry>geometry.gml</geometry>
+    <search_length_algorithm>
+        <type>fixed</type>
+        <value>1e-1</value>
+    </search_length_algorithm>
+    <processes>
+        <process>
+            <name>ThermoMechanics</name>
+            <type>THERMO_MECHANICS</type>
+            <integration_order>2</integration_order>
+            <constitutive_relation>
+                <type>CreepBGRa</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+                <a>2.0833333333333334e-06</a>
+                <n>4.9</n>
+                <sigma0>1e6</sigma0>
+                <q>54000</q>
+                <nonlinear_solver>
+                    <maximum_iterations>1000</maximum_iterations>
+                    <error_tolerance>2e-8</error_tolerance>
+                </nonlinear_solver>
+            </constitutive_relation>
+            <reference_solid_density>rho_sr</reference_solid_density>
+            <linear_thermal_expansion_coefficient>alpha</linear_thermal_expansion_coefficient>
+            <specific_heat_capacity>cs</specific_heat_capacity>
+            <thermal_conductivity>lambda</thermal_conductivity>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <temperature>temperature</temperature>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma_xx" output_name="sigma_xx"/>
+                <secondary_variable type="static" internal_name="sigma_yy" output_name="sigma_yy"/>
+                <secondary_variable type="static" internal_name="sigma_zz" output_name="sigma_zz"/>
+                <secondary_variable type="static" internal_name="sigma_xy" output_name="sigma_xy"/>
+                <secondary_variable type="static" internal_name="epsilon_xx" output_name="epsilon_xx"/>
+                <secondary_variable type="static" internal_name="epsilon_yy" output_name="epsilon_yy"/>
+                <secondary_variable type="static" internal_name="epsilon_zz" output_name="epsilon_zz"/>
+                <secondary_variable type="static" internal_name="epsilon_xy" output_name="epsilon_xy"/>
+            </secondary_variables>
+            <specific_body_force>0 -9.81</specific_body_force>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="ThermoMechanics">
+                <nonlinear_solver>basic_newton</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1e-11</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <output>
+                    <variables>
+                        <variable>displacement</variable>
+                        <variable>temperature</variable>
+                        <variable>sigma_xx</variable>
+                        <variable>sigma_yy</variable>
+                        <variable>sigma_zz</variable>
+                        <variable>sigma_xy</variable>
+                        <variable>epsilon_xx</variable>
+                        <variable>epsilon_yy</variable>
+                        <variable>epsilon_zz</variable>
+                        <variable>epsilon_xy</variable>
+                    </variables>
+                </output>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>43200e+2</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>86400.0e-3</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>10</repeat>
+                            <delta_t>86400.0e-2</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>500</repeat>
+                            <delta_t>86400.0</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>CreepAfterExcavation</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>10</repeat>
+                    <each_steps>2</each_steps>
+                </pair>
+                <pair>
+                    <repeat>1</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+        </output>
+    </time_loop>
+
+    <parameters>
+        <parameter>
+            <name>E</name>
+            <type>MeshElement</type>
+            <field_name>E</field_name>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>MeshElement</type>
+            <field_name>nu</field_name>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>MeshElement</type>
+            <field_name>rho</field_name>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>4e-5</value>
+        </parameter>
+        <parameter>
+            <name>cs</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>lambda</name>
+            <type>Constant</type>
+            <value>5</value>
+        </parameter>
+        <parameter>
+            <name>temperature_ic</name>
+            <type>MeshNode</type>
+            <field_name>T0</field_name>
+        </parameter>
+        <parameter>
+            <name>temperature_t</name>
+            <type>Constant</type>
+            <value>310</value>
+        </parameter>
+        <parameter>
+            <name>temperature_circ</name>
+            <type>Constant</type>
+            <value>300</value>
+        </parameter>
+        <parameter>
+            <name>temperature_b</name>
+            <type>Constant</type>
+            <value>320</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>u_bc</name>
+            <type>Constant</type>
+            <value>0.</value>
+        </parameter>
+        <parameter>
+            <name>s_top</name>
+            <type>Constant</type>
+            <value>-1.0e+7</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>1</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>left</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>u_bc</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>right</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>u_bc</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>u_bc</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Neumann</type>
+                    <component>1</component>
+                    <parameter>s_top</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>temperature_ic</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>top</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>temperature_t</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>bottom</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>temperature_b</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>geo</geometrical_set>
+                    <geometry>circle</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>temperature_circ</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_newton</name>
+            <type>Newton</type>
+            <max_iter>100</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i bicgstab -p jacobi -tol 1e-11 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>BiCGSTAB</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>

web/content/docs/benchmarks/creep-after-excavation-bgra/CreepAfterExcavation.md

++++
+date = "2018-08-10T11:46:45+01:00"
+title = "Creep in a rectangle domain with hole under pressure on its top surface"
+weight = 171
+project = "ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj"
+author = "Wenqing Wang"
+
+[menu]
+  [menu.benchmarks]
+    parent = "thermo-mechanics"
+
++++
+
+{{< data-link >}}
+
+-------------------------------------------------------------------------
+Following up the benchmark about the BGRa creep model described on the 
+[this page](https://www.opengeosys.org/docs/benchmarks/creepbgra/creepbrga/), 
+this example represents the creep in the near field of
+drift in the deep rock salt after excavation. The domain and the
+geometry are shown in the following figure:
+<figure>
+    <img src="../mesh.png" alt="Mesh and Geometry" id="fig_6"
+ style="height:400px;width:490px;">
+    <figcaption>Mesh and Geometry.</figcaption>
+</figure>
+
+The domain has two material groups, which are highlighted by different
+colors. The material group that is in the
+top part of the domain represents a cap rock type, while the other
+material group is for rock salt. The material properties of the rocks are given in the following table:
+
+---------------------- ---------- ----------- ------------
+|         |Cap rock |  Rock salt|   Unit|
+| --- |:---------:| -----:|----:|
+|Density            |    2000   |    2170 |       kg/m$^{3}$|
+|Young’s Modulus    |    7.0    |    7.65  |      GPa       |
+|Poisson ratio      |    0.3    |    2.7   |      -          |
+|Thermal conductivity |   5     |     5    |       W/(mK)    |
+---------------------- ---------- ----------- ------------
+
+The parameters of the BGRa creep model are $A=0.18\, \mbox{d}^{-1}$,
+$m=5$, $Q=54 \mbox{ kJ/mol}$.
+
+
+The width
+and the height of of the domain are 300 m and 340 m, respectively. The
+height of the cap rock portion is 40 m. The drift to be excavated has a
+radius of 50 m.
+
+Here we consider the creep of the rock after excavation. Therefore, we
+assume a steady state of stress and temperature of excavation at the
+beginning of the current simulation. For this assumption, the boundary
+conditions are given as:
+
+-   top boundary: $\tau_x = \sigma_x=0$, $\tau_y=\sigma_y=70$ MPa,
+    $T=310$ K.
+
+-   two lateral boundaries: normal displacement is fixed, and no heat
+    flux.
+
+-   bottom boundary: normal displacement is fixed, and $T=320$ K.
+
+-   circle of drift surface: traction free
+    (${ \mathbf\sigma}\cdot \mathbf n = 0$), and $T=300$ K as excavation
+    conditions.
+
+The boundary conditions of the mechanical process lead to a distinct
+stress field after excavation. Besides, heat capacity is neglected in
+order to pose a steady state temperature field after excavation. The
+initial stresses are zero and the initial temperature is linearly
+distributed from top to bottom with the boundary values of temperature.
+
+The time step sizes of the simulation are: One step of 0.001 day, 10 steps of 0.1 day, and the
+remaining steps of 1 day.
+
+The following three figures are plotted by using the results of the simulation of 1000 days creep. 
+The three figures display the distribution of horizontal and vertical stresses at times of
+108 days, 409 days and 1000 days, respectively. 
+
+In theses three figures,
+the left sub-figure show the time variations of horizontal and vertical stresses
+ at a position just close to the top of the drift, and the
+green vertical line in it marks the time of the displayed stress field.
+
+<figure>
+    <img src="../stress_xx_yy_20.png" alt="Stress distribution at the time of 109 days." id="fig_2">
+    <figcaption>Stress distribution at the time of 109 days.</figcaption>
+</figure>
+
+<figure>
+    <img src="../stress_xx_yy_50.png" alt="Stress distribution at the time of 409 days." id="fig_3">
+    <figcaption>Stress distribution at the time of 409 days.</figcaption>
+</figure>
+
+<figure>
+    <img src="../stress_xx_yy_110.png" alt="Stress distribution at the time of 1000 days." id="fig_4">
+    <figcaption>Stress distribution at the time of 1000 days.</figcaption>
+</figure>
+
+
+The above three figures show that the absolute value of the horizontal stress increase
+signification in the areas above and beneath the drift due to the creep,
+and the change of the vertical stress is slow compared to that of the
+horizontal stress.
+
+The following figure shows the strain distribution at the
+end of the simulation time at 1000 days.
+<figure>
+    <img src="../strain.png" alt="Strain distribution at the time of 1000 days." id="fig_5">
+    <figcaption>Strain distribution at the time of 1000 days.</figcaption>
+</figure>
+
+
+The steady-state temperature distribution is displayed in the following figure
+<figure>
+    <img src="../T.png" alt="Temperature distribution at the time of 1000 days." id="fig_6">
+    <figcaption>Temperature distribution at the time of 1000 days.</figcaption>
+</figure>
+
+## Note:
+For the automatic benchmarking, the time duration of creep is reduced to 50 days in order to reduce the run time.
+ If one wants to test this benchmark for 1000 days' creep, please change the end time in the tag of `<t_end> `
+in the project file as
+ ```
+<t_end>1000</t_end>
+```
\ No newline at end of file