From 9ebaaf2375fc7eef88e098b6324f6fc41cdb9afe Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 27 Oct 2014 15:52:26 +0000
Subject: [PATCH] Implement and use the long form of the nonquadratic membrane
energy
This form allows to set the shear correction factor kappa. It is still
not ruled out that we need to twiddle with this factor in order to
weaken the material a little bit.
[[Imported from SVN: r9937]]
---
dune/gfe/cosseratenergystiffness.hh | 39 +++++++++++++++++++++++++++--
1 file changed, 37 insertions(+), 2 deletions(-)
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index a7b31d83..6d7f7414 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -129,7 +129,7 @@ public:
mu_ = parameters.template get<double>("mu");
lambda_ = parameters.template get<double>("lambda");
- // Cosserat couple modulus, preferably 0
+ // Cosserat couple modulus
mu_c_ = parameters.template get<double>("mu_c");
// Length scale parameter
@@ -208,6 +208,40 @@ public:
+ (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
}
+ /** \brief The energy \f$ W_{mp}(\overline{U}) \f$, as written in
+ * the second equation of (4.4) in Neff's paper
+ */
+ RT longNonquadraticMembraneEnergy(const Dune::GFE::CosseratStrain<field_type,3,gridDim>& U) const
+ {
+ RT result = 0;
+
+ // shear-stretch energy
+ Dune::FieldMatrix<field_type,dim-1,dim-1> sym2x2;
+ for (int i=0; i<dim-1; i++)
+ for (int j=0; j<dim-1; j++)
+ sym2x2[i][j] = 0.5 * (U.matrix()[i][j] + U.matrix()[j][i]) - (i==j);
+
+ result += mu_ * sym2x2.frobenius_norm2();
+
+ // first order drill energy
+ Dune::FieldMatrix<field_type,dim-1,dim-1> skew2x2;
+ for (int i=0; i<dim-1; i++)
+ for (int j=0; j<dim-1; j++)
+ skew2x2[i][j] = 0.5 * (U.matrix()[i][j] - U.matrix()[j][i]);
+
+ result += mu_c_ * skew2x2.frobenius_norm2();
+
+
+ // classical transverse shear energy
+ result += kappa_ * (mu_ + mu_c_)/2 * (U.matrix()[2][0]*U.matrix()[2][0] + U.matrix()[2][1]*U.matrix()[2][1]);
+
+ // elongational stretch energy
+ RT detU = U.determinant();
+ result += (mu_*lambda_)/(2*mu_ + lambda_) * 0.5 * ((detU-1)*(detU-1) + (1.0/detU -1)*(1.0/detU -1));
+
+ return result;
+ }
+
RT curvatureEnergy(const Tensor3<field_type,3,3,3>& DR) const
{
#ifdef DONT_USE_CURL
@@ -329,7 +363,8 @@ energy(const Entity& element,
//energy += weight * thickness_ * quadraticMembraneEnergy(U.matrix());
energy += weight * thickness_ * longQuadraticMembraneEnergy(U);
#else
- energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
+ //energy += weight * thickness_ * nonquadraticMembraneEnergy(U);
+ energy += weight * thickness_ * longNonquadraticMembraneEnergy(U);
#endif
energy += weight * thickness_ * curvatureEnergy(DR);
energy += weight * std::pow(thickness_,3) / 12.0 * bendingEnergy(R,DR);
--
GitLab