From 034ed07547a294280baebe4b60b151817cc5b288 Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Thu, 2 Apr 2026 16:33:11 +0200 Subject: [PATCH 1/4] Documentation is added for cap zones and their return mappings --- .../custom_constitutive/README_cap.md | 210 +++++++++++++++++- 1 file changed, 208 insertions(+), 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md index e316633f397e..db903e48678b 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md +++ b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md @@ -66,7 +66,7 @@ The figure below shows a typical Mohr–Coulomb yield surface extended with tens Here, we need to convert the compression cap yield surface from $(p, q)$ coordinates to $(\sigma, \tau)$ coordinates. The conversion is to be followed ... -### Plastic Potential for the compression cap +### Plastic potential for the compression cap For the cap branch, plastic deformation is primarily volumetric (compaction), and the plastic potential is usually taken to be associated. The flow function is then: @@ -77,4 +77,210 @@ The derivative of the flow function is the: ```math \frac{\partial G_{cap}}{\partial \sigma} = \frac{2 q}{X^2} \frac{\partial q}{\partial \sigma} + 2 p \frac{\partial p}{\partial \sigma} -``` \ No newline at end of file +``` + +or + +```math + \frac{\partial G_{cap}}{\partial \sigma_i} = \frac{1}{3} G_{,p} + \frac{3 G_{,q}}{2 q} \left(\sigma_i - p \right) +``` + +### Cap corner point + +The point where the compression cap yield surface intersects the Mohr-Coulomb yield surface is called the cap corner point. This point can be calculated by extracting $q$ of Coulomb yield surface: + +```math + q = - \frac{6 \sin{\phi}}{3 - \sin{\phi}} p + \frac{6 c \cos⁡{\phi}}{3 - \sin{\phi}} +``` + +then subsituting in the cap yield surface, it leads to the following equation: + +```math + A p_{corner}^2 + B p_{corner} + C = 0 +``` + +Where, + +```math + A = 1 + b_1 a_2^2 +``` +```math + B = -2 b_1 a_2 c_2 +``` +```math + C = b_1 c_2^2 - c_1 +``` + +```math + b_1 = 1 / X^2 +``` +```math + c_1 = p_c^2 +``` +```math + a_2 = \frac{6 \sin{\phi}}{3 - \sin{\phi}} +``` +```math + c_2 = \frac{6 c \cos{\phi}}{3 - \sin{\phi}} +``` + +A second order polynomial equation needs to be solved, and the minimum root needs to be selected, because $p \lt 0$ + +```math + p_{corner} = \frac{ -B - \sqrt{B^2 - 4 A C}}{2A} +``` +then +```math + q_{corner} = -a_2 p + c_2 +``` + +### Return mapping from cap compression zone +The cap compression zone is the rigion where the trial principal stresses are, +1. Outside the compression cap yield surface +```math + \frac{q^2}{X^2} + p^2 - p_c^2 > 0 +``` +2. Under the line which passes from the cap corner point and in the direction normal to the flow function of the cap yield surface. +```math + q - q_{corner} - \left( G_{cap,p}/G_{cap,q} \right) (p - p_{corner}) < 0 +``` + +Then the trial principal stresses need to be mapped to the cap yield surface by: +```math + \sigma = \sigma^{trial} + \lambda C \frac{\partial G_{cap}}{\partial \sigma} +``` + +### Return mapping from cap corner zone +The cap compression zone is the rigion where the trial principal stresses are, + +1. above the line which passes from the cap corner point and in the direction normal to the flow function of the cap yield surface. +```math + q - q_{corner} - \left( G_{cap,p}/G_{cap,q} \right) (p - p_{corner}) > 0 +``` + +2. Under the line which passes from the cap corner point and in the direction normal to the flow function of the Coulomb yield surface. +```math + q - q_{corner} - \left( G_{MC,p}/G_{MC,q} \right) (p - p_{corner}) < 0 +``` + +Then the trial principal stresses need to be mapped to the cap yield surface by: +```math + \sigma = \sigma^{trial} + \lambda_{MC} C \frac{\partial G_{MC}}{\partial \sigma} + + \lambda_{cap} C \frac{\partial G_{cap}}{\partial \sigma} +``` + +Subsituting this traisl stresses in compression cap and Coulomb yield surfaces, it leads to two equations and two unknowns. +```math + c_1 \lambda_{MC} + c_2 \lambda_{cap} = c_3 +``` +```math + c_4 \lambda_{MC}^2 + c_5 \lambda_{cap}^2 + c_6 \lambda_{MC} + c_7 \lambda_{cap} + + c_8 \lambda_{MC} \lambda_{cap} = c_9 +``` + +Combining those two equations, it leads to a second order polynomial equation for $\lambda_{cap}$. + +```math + A \lambda_{cap}^2 + B \lambda_{cap} + C = 0 +``` + +where, + + +```math + c_0 = \frac{6 \sin{\phi}}{3 - \sin{\phi}} +``` +```math + \left[ p_{MC}^{cor}, q_{MC}^{cor} \right]^T = C G_{MC} +``` +```math + \left[ p_{cap}^{cor}, q_{cap}^{cor} \right]^T = C G_{cap} +``` +```math + c_1 = q_{MC}^{cor} + c_0 p_{MC}^{cor} +``` +```math + c_2 = q_{cap}^{cor} + c_0 p_{cap}^{cor} +``` +```math + c_3 = -F_{MC} \left( \sigma^{trial} \right) +``` +```math + c_4 = q_{MC}^{cor^2} / X^2 + p_{MC}^{cor^2} +``` +```math + c_5 = q_{cap}^{cor^2} / X^2 + p_{cap}^{cor^2} +``` +```math + c_6 = 2 \left( q^{trial} p_{MC}^{cor} / X^2 + p^{trial} p_{MC}^{cor} \right) +``` +```math + c_7 = 2 \left( q^{trial} q_{cap}^{cor} / X^2 + p^{trial} p_{cap}^{cor} \right) +``` +```math + c_8 = 2 \left( q_{MC}^{cor} q_{cap}^{cor} / X^2 + p_{MC}^{cor} p_{cap}^{cor} \right) +``` +```math + c_9 = -F_{cap} \left( \sigma^{trial} \right) +``` + +```math + A = \frac{c_2}{c_1} \left( \frac{c_2 c_4}{c_1} - c_8 \right) + c_5 +``` +```math + B = \frac{1}{c_1} \left( -2 \frac{c_2 c_3 c_4}{c_1} - c_2 c_6 + c_3 c_8 \right) + c_7 +``` +```math + C = \frac{c_3}{c_1} \left( \frac{c_3 c_4}{c_1} + c_6 \right) - c9 +``` + +Then, solving the recond order polynomial, it gives +```math + \lambda_{cap} = \frac{-B + \sqrt{B^2 - 4 A C}}{2A} +``` +```math + \lambda_{MC} = \frac{c_3 - c_2 \lambda_{cap}}{c_1} +``` + +### Plastic multiplier for compression cap +The plastic multiplier above, $\lambda_{cap}$ needs to be calculated. The mapping is, +```math + \sigma = \sigma^{trial} + \lambda_{cap} C \frac{\partial G_{cap}}{\partial \sigma} +``` + +Now, by considering the cap yield surface function, and extracting $p$ and $q$ in the form of $\sigma$'s, +```math + \frac{1}{2 X^2} \left[ \left( \sigma_1 - \sigma_2 \right)^2 + \left( \sigma_2 - \sigma_3 \right)^2 + \left( \sigma_3 - \sigma_1 \right)^2 \right] + \frac{1}{9} \left( \sigma_1 + \sigma_2 + \sigma_3 \right)^2 - p_c^2 = 0 +``` +We define the following vectors, +```math + \Delta \sigma = \left[\sigma_1 - \sigma_2 \; , \; \sigma_2 - \sigma_3 \; , \; \sigma_3 - \sigma_1 \right]^T +``` +```math + \Delta \sigma^{cor} = \left[\sigma_1^{cor} - \sigma_2^{cor} \; , \; \sigma_2^{cor} - \sigma_3^{cor} \; , \; \sigma_3^{cor} - \sigma_1^{cor} \right]^T +``` + +where +```math + \sigma^{cor} = \lambda_{cap} C \frac{\partial G_{cap}}{\partial \sigma} +``` + +Subsituting the stresses with the mapped stresses, we get the following relations. + +```math + A = \frac{\Delta \sigma^{cor} \cdot \Delta \sigma^{cor}}{2 X^2} + \frac{1}{9} \left( \sum_{i=1}^3{\sigma_i^{cor}} \right)^2 +``` +```math + B = \frac{ \Delta \sigma \cdot \Delta \sigma^{cor}}{X^2} + \frac{2}{9} + \sum_{i=1}^3{\sigma_i} \sum_{i=1}^3{\Delta \sigma_i^{cor}} +``` +```math + C = F_{cap} \left( \sigma^{trial} \right) +``` + +the solution of this equation, gives: + + ```math + \lambda_{cap} = \frac{-B + \sqrt{B^2 - 4 A C}}{2A} +``` + From df90deebc3ac97ed4614c60797e6909c5455380b Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Fri, 3 Apr 2026 11:51:44 +0200 Subject: [PATCH 2/4] fixes --- .../custom_constitutive/README_cap.md | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md index db903e48678b..b5fb52acf1b8 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md +++ b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md @@ -21,7 +21,7 @@ where: In stress-invariant form, the MC yield function is typically written as: ```math - F_{MC}(p, q) = q - \frac{6 \sin{\phi}}{3 - \sin{\phi}} p - \frac{6 c \cos⁡{\phi}}{3 - \sin{\phi}} + F_{MC}(p, q) = q + \frac{6 \sin{\phi}}{3 - \sin{\phi}} p - \frac{6 c \cos⁡{\phi}}{3 - \sin{\phi}} ``` where: @@ -73,7 +73,7 @@ For the cap branch, plastic deformation is primarily volumetric (compaction), an ```math G_{cap} \left(p, q \right) = F_{cap} \left(p, q \right) ``` -The derivative of the flow function is the: +The derivative of the flow function is: ```math \frac{\partial G_{cap}}{\partial \sigma} = \frac{2 q}{X^2} \frac{\partial q}{\partial \sigma} + 2 p \frac{\partial p}{\partial \sigma} @@ -87,7 +87,7 @@ or ### Cap corner point -The point where the compression cap yield surface intersects the Mohr-Coulomb yield surface is called the cap corner point. This point can be calculated by extracting $q$ of Coulomb yield surface: +The point where the compression cap yield surface intersects the Mohr-Coulomb yield surface is called the cap corner point. This point can be calculated by solving the Coulomb yield surface for $q$: ```math q = - \frac{6 \sin{\phi}}{3 - \sin{\phi}} p + \frac{6 c \cos⁡{\phi}}{3 - \sin{\phi}} @@ -151,7 +151,7 @@ Then the trial principal stresses need to be mapped to the cap yield surface by: ``` ### Return mapping from cap corner zone -The cap compression zone is the rigion where the trial principal stresses are, +The cap compression zone is the region where the trial principal stresses are, 1. above the line which passes from the cap corner point and in the direction normal to the flow function of the cap yield surface. ```math @@ -163,13 +163,13 @@ The cap compression zone is the rigion where the trial principal stresses are, q - q_{corner} - \left( G_{MC,p}/G_{MC,q} \right) (p - p_{corner}) < 0 ``` -Then the trial principal stresses need to be mapped to the cap yield surface by: +Then the trial principal stresses need to be mapped to the cap yield surface and Coulomb yield surface by: ```math \sigma = \sigma^{trial} + \lambda_{MC} C \frac{\partial G_{MC}}{\partial \sigma} + \lambda_{cap} C \frac{\partial G_{cap}}{\partial \sigma} ``` -Subsituting this traisl stresses in compression cap and Coulomb yield surfaces, it leads to two equations and two unknowns. +Substuting this trial stresses in compression cap and Coulomb yield surfaces, it leads to two equations and two unknowns. ```math c_1 \lambda_{MC} + c_2 \lambda_{cap} = c_3 ``` @@ -212,7 +212,7 @@ where, c_5 = q_{cap}^{cor^2} / X^2 + p_{cap}^{cor^2} ``` ```math - c_6 = 2 \left( q^{trial} p_{MC}^{cor} / X^2 + p^{trial} p_{MC}^{cor} \right) + c_6 = 2 \left( q^{trial} q_{MC}^{cor} / X^2 + p^{trial} p_{MC}^{cor} \right) ``` ```math c_7 = 2 \left( q^{trial} q_{cap}^{cor} / X^2 + p^{trial} p_{cap}^{cor} \right) @@ -231,10 +231,10 @@ where, B = \frac{1}{c_1} \left( -2 \frac{c_2 c_3 c_4}{c_1} - c_2 c_6 + c_3 c_8 \right) + c_7 ``` ```math - C = \frac{c_3}{c_1} \left( \frac{c_3 c_4}{c_1} + c_6 \right) - c9 + C = \frac{c_3}{c_1} \left( \frac{c_3 c_4}{c_1} + c_6 \right) - c_9 ``` -Then, solving the recond order polynomial, it gives +Then, solving the second order polynomial, it gives ```math \lambda_{cap} = \frac{-B + \sqrt{B^2 - 4 A C}}{2A} ``` @@ -272,7 +272,7 @@ Subsituting the stresses with the mapped stresses, we get the following relation ``` ```math B = \frac{ \Delta \sigma \cdot \Delta \sigma^{cor}}{X^2} + \frac{2}{9} - \sum_{i=1}^3{\sigma_i} \sum_{i=1}^3{\Delta \sigma_i^{cor}} + \sum_{i=1}^3{\sigma_i} \sum_{i=1}^3{\sigma_i^{cor}} ``` ```math C = F_{cap} \left( \sigma^{trial} \right) From 6998df969d79e0a57ec8dc86d5271da5820693cb Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Tue, 7 Apr 2026 22:52:06 +0200 Subject: [PATCH 3/4] fixes based on review comments --- .../custom_constitutive/README_cap.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md index b5fb52acf1b8..7aedd8dd3b0d 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md +++ b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md @@ -82,7 +82,7 @@ The derivative of the flow function is: or ```math - \frac{\partial G_{cap}}{\partial \sigma_i} = \frac{1}{3} G_{,p} + \frac{3 G_{,q}}{2 q} \left(\sigma_i - p \right) + \frac{\partial G_{cap}}{\partial \sigma_i} = \frac{1}{3} \frac{\partial G_{cap}}{\partial p} + \frac{3}{2 q} \frac{\partial G_{cap}}{\partial q} \left(\sigma_i - p \right) ``` ### Cap corner point @@ -124,7 +124,7 @@ Where, c_2 = \frac{6 c \cos{\phi}}{3 - \sin{\phi}} ``` -A second order polynomial equation needs to be solved, and the minimum root needs to be selected, because $p \lt 0$ +A second order polynomial equation needs to be solved, and the minimum root needs to be selected, ```math p_{corner} = \frac{ -B - \sqrt{B^2 - 4 A C}}{2A} @@ -153,12 +153,12 @@ Then the trial principal stresses need to be mapped to the cap yield surface by: ### Return mapping from cap corner zone The cap compression zone is the region where the trial principal stresses are, -1. above the line which passes from the cap corner point and in the direction normal to the flow function of the cap yield surface. +1. above the line which passes from the cap corner point and in the direction normal to the cap flow function $G_{cap}$. ```math q - q_{corner} - \left( G_{cap,p}/G_{cap,q} \right) (p - p_{corner}) > 0 ``` -2. Under the line which passes from the cap corner point and in the direction normal to the flow function of the Coulomb yield surface. +2. Under the line which passes from the cap corner point and in the direction normal to the Coulomb flow function $G_{MC}$. ```math q - q_{corner} - \left( G_{MC,p}/G_{MC,q} \right) (p - p_{corner}) < 0 ``` From 9710c5d9ac9d9b5c47b2c2ba8903fd146ad4da7a Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Tue, 7 Apr 2026 22:56:57 +0200 Subject: [PATCH 4/4] fix --- .../GeoMechanicsApplication/custom_constitutive/README_cap.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md index 7aedd8dd3b0d..8fa8c81efe7d 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md +++ b/applications/GeoMechanicsApplication/custom_constitutive/README_cap.md @@ -151,7 +151,7 @@ Then the trial principal stresses need to be mapped to the cap yield surface by: ``` ### Return mapping from cap corner zone -The cap compression zone is the region where the trial principal stresses are, +The cap compression zone is the region where the trial $\left(p, q \right)$ stress invariants are demarcated by the following two lines, 1. above the line which passes from the cap corner point and in the direction normal to the cap flow function $G_{cap}$. ```math @@ -262,7 +262,7 @@ We define the following vectors, where ```math - \sigma^{cor} = \lambda_{cap} C \frac{\partial G_{cap}}{\partial \sigma} + \sigma^{cor} = C \frac{\partial G_{cap}}{\partial \sigma} ``` Subsituting the stresses with the mapped stresses, we get the following relations.