From b7d348ce943f7c6636e3b99bf27048380fd5fdfe Mon Sep 17 00:00:00 2001
From: tnagel <nagelt@tcd.ie>
Date: Fri, 9 Aug 2019 14:08:10 +0200
Subject: [PATCH] Changed Lubby2 from rate form to incremental form.
---
MaterialLib/SolidModels/Lubby2.cpp | 48 ++++++++++++++----------------
1 file changed, 23 insertions(+), 25 deletions(-)
diff --git a/MaterialLib/SolidModels/Lubby2.cpp b/MaterialLib/SolidModels/Lubby2.cpp
index 7d844ec1da9..c824ad0d7be 100644
--- a/MaterialLib/SolidModels/Lubby2.cpp
+++ b/MaterialLib/SolidModels/Lubby2.cpp
@@ -25,7 +25,7 @@ namespace Lubby2
template <int DisplacementDim>
Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
Lubby2<DisplacementDim>::KelvinVectorSize>
-calculatedGdEBurgers(double const dt)
+calculatedGdEBurgers()
{
Eigen::Matrix<double, Lubby2<DisplacementDim>::JacobianResidualSize,
Lubby2<DisplacementDim>::KelvinVectorSize>
@@ -35,17 +35,16 @@ calculatedGdEBurgers(double const dt)
dGdE.template topLeftCorner<Lubby2<DisplacementDim>::KelvinVectorSize,
Lubby2<DisplacementDim>::KelvinVectorSize>()
.diagonal()
- .setConstant(-2. / dt);
+ .setConstant(-2.);
return dGdE;
}
template <int DisplacementDim, typename LinearSolver>
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> tangentStiffnessA(
- double const GM0, double const KM0, double const dt,
- LinearSolver const& linear_solver)
+ double const GM0, double const KM0, LinearSolver const& linear_solver)
{
// Calculate dGdE for time step
- auto const dGdE = calculatedGdEBurgers<DisplacementDim>(dt);
+ auto const dGdE = calculatedGdEBurgers<DisplacementDim>();
// Consistent tangent from local Newton iteration of material
// functionals.
@@ -185,7 +184,6 @@ Lubby2<DisplacementDim>::integrateStress(
KelvinMatrix C =
tangentStiffnessA<DisplacementDim>(local_lubby2_properties.GM0,
local_lubby2_properties.KM0,
- dt,
linear_solver);
// Hydrostatic part for the stress and the tangent.
@@ -219,21 +217,21 @@ void Lubby2<DisplacementDim>::calculateResidualBurgers(
{
// calculate stress residual
res.template segment<KelvinVectorSize>(0).noalias() =
- (stress_curr - stress_t) / dt -
- 2. / dt *
- ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
- (strain_Max_curr - strain_Max_t));
+ (stress_curr - stress_t) -
+ 2. * ((strain_curr - strain_t) - (strain_Kel_curr - strain_Kel_t) -
+ (strain_Max_curr - strain_Max_t));
// calculate Kelvin strain residual
res.template segment<KelvinVectorSize>(KelvinVectorSize).noalias() =
- 1. / dt * (strain_Kel_curr - strain_Kel_t) -
- 1. / (2. * properties.etaK) * (properties.GM0 * stress_curr -
- 2. * properties.GK * strain_Kel_curr);
+ (strain_Kel_curr - strain_Kel_t) -
+ dt / (2. * properties.etaK) *
+ (properties.GM0 * stress_curr -
+ 2. * properties.GK * strain_Kel_curr);
// calculate Maxwell strain residual
res.template segment<KelvinVectorSize>(2 * KelvinVectorSize).noalias() =
- 1. / dt * (strain_Max_curr - strain_Max_t) -
- 0.5 * properties.GM0 / properties.etaM * stress_curr;
+ (strain_Max_curr - strain_Max_t) -
+ dt * 0.5 * properties.GM0 / properties.etaM * stress_curr;
}
template <int DisplacementDim>
@@ -252,23 +250,23 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
// build G_11
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
.diagonal()
- .setConstant(1. / dt);
+ .setConstant(1.);
// build G_12
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0, KelvinVectorSize)
.diagonal()
- .setConstant(2. / dt);
+ .setConstant(2.);
// build G_13
Jac.template block<KelvinVectorSize, KelvinVectorSize>(0,
2 * KelvinVectorSize)
.diagonal()
- .setConstant(2. / dt);
+ .setConstant(2.);
// build G_21
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize, 0)
.noalias() =
- -0.5 * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
+ -0.5 * dt * properties.GM0 / properties.etaK * KelvinMatrix::Identity();
if (s_eff > 0.)
{
KelvinVector const eps_K_aid =
@@ -281,15 +279,15 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
properties.etaK / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
0)
- .noalias() += 0.5 * eps_K_aid * dmu_vK.transpose() +
- 1. / properties.etaK * eps_K_i * dG_K.transpose();
+ .noalias() += 0.5 * dt * eps_K_aid * dmu_vK.transpose() +
+ dt / properties.etaK * eps_K_i * dG_K.transpose();
}
// build G_22
Jac.template block<KelvinVectorSize, KelvinVectorSize>(KelvinVectorSize,
KelvinVectorSize)
.diagonal()
- .setConstant(1. / dt + properties.GK / properties.etaK);
+ .setConstant(1. + dt * properties.GK / properties.etaK);
// nothing to do for G_23
@@ -297,14 +295,14 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
0)
.noalias() =
- -0.5 * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
+ -0.5 * dt * properties.GM0 / properties.etaM * KelvinMatrix::Identity();
if (s_eff > 0.)
{
KelvinVector const dmu_vM = 1.5 * _mp.mvM(t, x)[0] * properties.GM0 *
properties.etaM / s_eff * sig_i;
Jac.template block<KelvinVectorSize, KelvinVectorSize>(
2 * KelvinVectorSize, 0)
- .noalias() += 0.5 * properties.GM0 /
+ .noalias() += 0.5 * dt * properties.GM0 /
(properties.etaM * properties.etaM) * sig_i *
dmu_vM.transpose();
}
@@ -315,7 +313,7 @@ void Lubby2<DisplacementDim>::calculateJacobianBurgers(
Jac.template block<KelvinVectorSize, KelvinVectorSize>(2 * KelvinVectorSize,
2 * KelvinVectorSize)
.diagonal()
- .setConstant(1. / dt);
+ .setConstant(1.);
}
template class Lubby2<2>;
--
GitLab