From f1cdcabb14b6f427883be9dd2708b6e2cad1c3fc Mon Sep 17 00:00:00 2001
From: Patrick Jaap <patrick.jaap@tu-dresden.de>
Date: Wed, 29 Jun 2022 11:05:19 +0200
Subject: [PATCH] Fix in Mooney-Rivlin density: Use proper 4th term for Ciarlet
In Ciarlet the last factor is left unspecified. I could not follow the
choice used in this implementation. I guess it was intended to minimize
the energy under a hydrostatic stress F = t*I at the unit matrix, i.e.,
for t=1. A simple calculation shows that (independent on the dimension)
d = 2a + 4b + 2c
has to hold in this case. Values in the test were adjusted.
---
dune/elasticity/materials/mooneyrivlindensity.hh | 10 ++++++----
test/mooneyrivlintest.cc | 8 ++++----
2 files changed, 10 insertions(+), 8 deletions(-)
diff --git a/dune/elasticity/materials/mooneyrivlindensity.hh b/dune/elasticity/materials/mooneyrivlindensity.hh
index 2b4c35d..37ba272 100644
--- a/dune/elasticity/materials/mooneyrivlindensity.hh
+++ b/dune/elasticity/materials/mooneyrivlindensity.hh
@@ -64,13 +64,14 @@ public:
/* The Mooney-Rivlin-Density is given as a function of the eigenvalues of the right Cauchy-Green-Deformation tensor C = F^TF
or the right Cauchy-Green-Deformation tensor B = FF^T.
C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
-
+
There are three Mooney-Rivlin-Variants:
ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF^2 + mooneyrivlin_c*(detF)^2 -
- ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF)
+ (2*mooneyrivlin_a + 4*mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF),
+ where the last term is chosen s.t. W( t*I ) is minimal for t=1
log: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2
square: W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2
-
+
For log and square: I1 and I2 are the first two invariants of C (or B), multiplied with a factor depending on det(F):
I1 = (det(F)^(-2/dim)) * [ first invariant of C ]
= (det(F)^(-2/dim)) * (sum of all eigenvalues of C)
@@ -99,7 +100,8 @@ public:
gradientInverse.invert();
field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
using std::log;
- return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
+ return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF
+ - (2.0*mooneyrivlin_a + 4.0*mooneyrivlin_b + 2.0*mooneyrivlin_c)*log(detF);
}
else
{
diff --git a/test/mooneyrivlintest.cc b/test/mooneyrivlintest.cc
index 8af6cba..330aedf 100644
--- a/test/mooneyrivlintest.cc
+++ b/test/mooneyrivlintest.cc
@@ -167,10 +167,10 @@ int main (int argc, char *argv[])
parameters["mooneyrivlin_a"] = "2.42e+6";
parameters["mooneyrivlin_b"] = "6.52e+6";
parameters["mooneyrivlin_c"] = "-7.34e+6";
- expectedEnergy = 68753239.6;
- expectedGradientTwoNorm = 31244643.1;
- expectedGradientInfinityNorm = 9673572.39;
- expectedMatrixFrobeniusNorm = 1.67660965e+09;
+ expectedEnergy = 76302830.4;
+ expectedGradientTwoNorm = 40670527.3;
+ expectedGradientInfinityNorm = 11116511.8;
+ expectedMatrixFrobeniusNorm = 2.1978108e+09;
int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
return testCiarlet + testLog + testSquare;
--
GitLab