add OOP designed NeoHookean material

MatMechLab · Apr 11, 2021 · 67acda1 · 67acda1
1 parent 3227575
commit 67acda1
Show file tree

Hide file tree

Showing 9 changed files with 370 additions and 0 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -219,6 +219,9 @@ set(src ${src} src/MateSystem/ConstDiffusionMaterial.cpp)
 # for elastic material
 set(inc ${inc} include/MateSystem/LinearElasticMaterial.h)
 set(src ${src} src/MateSystem/LinearElasticMaterial.cpp)
+### For NeoHookean material
+set(inc ${inc} include/MateSystem/NeoHookeanMaterial.h)
+set(src ${src} src/MateSystem/NeoHookeanMaterial.cpp)
 ### For free energy based materials
 set(inc ${inc} include/MateSystem/FreeEnergyMaterialBase.h)
 #   for double well free energy materials

diff --git a/include/MateSystem/BulkMateSystem.h b/include/MateSystem/BulkMateSystem.h
@@ -42,12 +42,14 @@
 #include "MateSystem/ConstDiffusionMaterial.h"
 #include "MateSystem/DoubleWellFreeEnergyMaterial.h"
 #include "MateSystem/LinearElasticMaterial.h"
+#include "MateSystem/NeoHookeanMaterial.h"
 #include "MateSystem/MieheFractureMaterial.h"
 
 class BulkMateSystem: public ConstPoissonMaterial,
                       public ConstDiffusionMaterial,
                       public DoubleWellFreeEnergyMaterial,
                       public LinearElasticMaterial,
+                      public NeoHookeanMaterial,
                       public MieheFractureMaterial{
 public:
     BulkMateSystem();

diff --git a/include/MateSystem/NeoHookeanMaterial.h b/include/MateSystem/NeoHookeanMaterial.h
@@ -0,0 +1,42 @@
+//****************************************************************
+//* This file is part of the AsFem framework
+//* A Simple Finite Element Method program (AsFem)
+//* All rights reserved, Yang Bai @ CopyRight 2021
+//* https://github.com/yangbai90/AsFem.git
+//* Licensed under GNU GPLv3, please see LICENSE for details
+//* https://www.gnu.org/licenses/gpl-3.0.en.html
+//****************************************************************
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+//+++ Author : Yang Bai
+//+++ Date   : 2021.04.11
+//+++ Purpose: Implement the calculation of neo-hookean type
+//+++          hyperlelastic material
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#pragma once
+
+#include "MateSystem/MechanicsMaterialBase.h"
+
+class NeoHookeanMaterial: public MechanicsMaterialBase{
+public:
+    virtual void InitMaterialProperties(const int &nDim,const Vector3d &gpCoord,const vector<double> &InputParams,
+                                        const vector<double> &gpU,const vector<double> &gpUdot,
+                                        const vector<Vector3d> &gpGradU,const vector<Vector3d> &gpGradUdot,
+                                        Materials &Mate) override;
+
+    virtual void ComputeMaterialProperties(const double &t, const double &dt,const int &nDim,
+                                           const Vector3d &gpCoord,const vector<double> &InputParams,
+                                           const vector<double> &gpU,const vector<double> &gpUOld,
+                                           const vector<double> &gpUdot,const vector<double> &gpUdotOld,
+                                           const vector<Vector3d> &gpGradU,const vector<Vector3d> &gpGradUOld,
+                                           const vector<Vector3d> &gpGradUdot,const vector<Vector3d> &gpGradUdotOld,
+                                           const Materials &MateOld, Materials &Mate) override;
+
+private:
+    virtual void ComputeStrain(const int &nDim,const vector<Vector3d> &GradDisp, RankTwoTensor &Strain) override;
+    virtual void ComputeStressAndJacobian(const vector<double> &InputParams,const RankTwoTensor &Strain,RankTwoTensor &Stress,RankFourTensor &Jacobian) override;
+
+private:
+    RankTwoTensor _GradU,_Strain,_Stress,_I,_devStress,_F;
+    RankTwoTensor _C,_Cinv;
+    RankFourTensor _Jac;
+};
diff --git a/src/MateSystem/InitBulkMateLibs.cpp b/src/MateSystem/InitBulkMateLibs.cpp
@@ -37,6 +37,10 @@ void BulkMateSystem::InitBulkMateLibs(const MateType &imate,const int &mateindex
             LinearElasticMaterial::InitMaterialProperties(nDim,gpCoord,_BulkMateBlockList[mateindex-1]._Parameters,
                                                           gpU,gpUdot,gpGradU,gpGradUdot,_Materials);
             break;
+        case MateType::NEOHOOKEANMATE:
+            NeoHookeanMaterial::InitMaterialProperties(nDim,gpCoord,_BulkMateBlockList[mateindex-1]._Parameters,
+                                                       gpU,gpUdot,gpGradU,gpGradUdot,_Materials);
+            break;
         case MateType::MIEHEFRACTUREMATE:
             MieheFractureMaterial::InitMaterialProperties(nDim,gpCoord,_BulkMateBlockList[mateindex-1]._Parameters,
                                                           gpU,gpUdot,gpGradU,gpGradUdot,_Materials);

diff --git a/src/MateSystem/NeoHookeanMaterial.cpp b/src/MateSystem/NeoHookeanMaterial.cpp
@@ -0,0 +1,86 @@
+//****************************************************************
+//* This file is part of the AsFem framework
+//* A Simple Finite Element Method program (AsFem)
+//* All rights reserved, Yang Bai @ CopyRight 2021
+//* https://github.com/yangbai90/AsFem.git
+//* Licensed under GNU GPLv3, please see LICENSE for details
+//* https://www.gnu.org/licenses/gpl-3.0.en.html
+//****************************************************************
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+//+++ Author : Yang Bai
+//+++ Date   : 2021.04.11
+//+++ Purpose: Implement the calculation of neo-hookean type
+//+++          hyperlelastic material
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#include "MateSystem/NeoHookeanMaterial.h"
+
+void NeoHookeanMaterial::InitMaterialProperties(const int &nDim, const Vector3d &gpCoord,
+                                                const vector<double> &InputParams, const vector<double> &gpU,
+                                                const vector<double> &gpUdot, const vector<Vector3d> &gpGradU,
+                                                const vector<Vector3d> &gpGradUdot, Materials &Mate) {
+    // Here we do not consider any initial internal strains, stress
+    if(nDim||gpCoord(1)||InputParams.size()||gpU[0]||gpUdot[0]||
+       gpGradU[0](1)||gpGradUdot[0](1)||Mate.ScalarMaterials.size()){}
+}
+//*****************************************************
+void NeoHookeanMaterial::ComputeStrain(const int &nDim, const vector<Vector3d> &GradDisp, RankTwoTensor &Strain) {
+    // here we calculate the deformation gradient tensor as well as euler-lagrange strain tensor
+    if(nDim==2){
+        _GradU.SetFromGradU(GradDisp[1],GradDisp[2]);
+    }
+    else if(nDim==3){
+        _GradU.SetFromGradU(GradDisp[1],GradDisp[2],GradDisp[3]);
+    }
+    _I.SetToIdentity();
+    _F=_GradU+_I;// F=I+U_{i,j}
+    _C=_F.Transpose()*_F; //C=F^T*F
+    _Cinv=_C.Inverse();
+    Strain=(_C-_I)*0.5;
+}
+//*****************************************************
+void NeoHookeanMaterial::ComputeStressAndJacobian(const vector<double> &InputParams, const RankTwoTensor &Strain,
+                                                  RankTwoTensor &Stress, RankFourTensor &Jacobian) {
+    // here Strain is E
+    if(Strain(1,1)){}// get rid of unused warning
+
+
+    double EE=InputParams[0];
+    double nu=InputParams[1];
+
+    double lambda=EE*nu/((1+nu)*(1-2*nu));
+    double mu=EE/(2*(1+nu));
+    double J=_F.Det();
+
+    Stress=(_I-_Cinv)*mu+_Cinv*lambda*(J-1)*J;
+    Jacobian=_Cinv.ODot(_Cinv)*mu*2.0
+            +_Cinv.CrossDot(_Cinv)*lambda*(2*J-1)*J
+            -_Cinv.ODot(_Cinv)*lambda*(J-1)*J*2;
+}
+//*****************************************************
+void NeoHookeanMaterial::ComputeMaterialProperties(const double &t, const double &dt, const int &nDim,
+                                                   const Vector3d &gpCoord, const vector<double> &InputParams,
+                                                   const vector<double> &gpU, const vector<double> &gpUOld,
+                                                   const vector<double> &gpUdot, const vector<double> &gpUdotOld,
+                                                   const vector<Vector3d> &gpGradU, const vector<Vector3d> &gpGradUOld,
+                                                   const vector<Vector3d> &gpGradUdot,const vector<Vector3d> &gpGradUdotOld,
+                                                   const Materials &MateOld,Materials &Mate) {
+    //***************************************************************
+    //*** get rid of unused warning
+    //***************************************************************
+    if(t||dt||gpCoord(1)||gpU[0]||gpUOld[0]||
+       gpUdot[0]||gpUdotOld[0]||gpGradU[0](1)||gpGradUOld[0](1)||
+       gpGradUdot[0](1)||gpGradUdotOld[0](1)||MateOld.ScalarMaterials.size()){}// get rid of unused warning
+
+    ComputeStrain(nDim,gpGradU,_Strain);
+    ComputeStressAndJacobian(InputParams,_Strain,_Stress,_Jac);
+
+    _devStress=_Stress-_I*(_Stress.Trace()/3.0);
+
+    Mate.ScalarMaterials["vonMises"]=sqrt(1.5*_devStress.DoubleDot(_devStress));
+
+    Mate.Rank2Materials["strain"]=_Strain;
+    Mate.Rank2Materials["stress"]=_Stress;
+    Mate.Rank4Materials["jacobian"]=_Jac;
+
+}
diff --git a/src/MateSystem/RunBulkMateLibs.cpp b/src/MateSystem/RunBulkMateLibs.cpp
@@ -48,6 +48,12 @@ void BulkMateSystem::RunBulkMateLibs(const MateType &imate,const int &mateindex,
                                                          gpGradU,gpGradUOld,gpGradUdot,gpGradUdotOld,
                                                          _MaterialsOld,_Materials);
         break;
+    case MateType::NEOHOOKEANMATE:
+        NeoHookeanMaterial::ComputeMaterialProperties(t,dt,nDim,gpCoord,_BulkMateBlockList[mateindex-1]._Parameters,
+                                                      gpU,gpUOld,gpUdot,gpUdotOld,
+                                                      gpGradU,gpGradUOld,gpGradUdot,gpGradUdotOld,
+                                                      _MaterialsOld,_Materials);
+        break;
     case MateType::MIEHEFRACTUREMATE:
         MieheFractureMaterial::ComputeMaterialProperties(t,dt,nDim,gpCoord,_BulkMateBlockList[mateindex-1]._Parameters,
                                                          gpU,gpUOld,gpUdot,gpUdotOld,

diff --git a/tests/mechanics/neohookean2d.i b/tests/mechanics/neohookean2d.i
@@ -0,0 +1,73 @@
+*** This is an input file for the compressive neohookean model
+
+[mesh]
+  type=asfem
+  dim=2
+  xmax=5
+  ymax=5
+  nx=50
+  ny=50
+  meshtype=quad4
+[end]
+
+[dofs]
+name=ux uy
+[end]
+
+[projection]
+scalarmate=vonMises
+rank2mate=stress strain
+[end]
+
+[elmts]
+  [mechanics]
+    type=mechanics
+    dofs=ux uy
+    mate=neohookean
+    domain=alldomain
+  [end]
+[end]
+
+[mates]
+  [neohookean]
+    type=neohookean
+    params=100.0 0.3
+  [end]
+[end]
+
+
+
+[bcs]
+  [FixUx]
+    type=dirichlet
+    dof=ux
+    boundary=bottom
+    value=0.0
+  [end]
+  [FixUy]
+    type=dirichlet
+    dof=uy
+    boundary=bottom top
+    value=0.0
+  [end]
+  [loadUx]
+    type=dirichlet
+    dof=ux
+    value=1.0*t
+    boundary=top
+  [end]
+[end]
+
+[timestepping]
+  type=be
+  dt=2.0e-3
+  endtime=1.0e-1
+  adaptive=false
+  optiters=3
+  dtmax=1.0e-1
+[end]
+
+[job]
+  type=transient
+  debug=dep
+[end]
diff --git a/tests/mechanics/neohookean2dtensile.i b/tests/mechanics/neohookean2dtensile.i
@@ -0,0 +1,75 @@
+*** This is an input file for the compressive neohookean model
+
+[mesh]
+  type=asfem
+  dim=2
+  xmax=2
+  ymax=2
+  nx=80
+  ny=80
+  meshtype=quad4
+[end]
+
+
+
+[dofs]
+name=ux uy
+[end]
+
+[projection]
+scalarmate=vonMises
+rank2mate=stress strain
+[end]
+
+[elmts]
+  [mechanics]
+    type=mechanics
+    dofs=ux uy
+    mate=neohookean
+    domain=alldomain
+  [end]
+[end]
+
+[mates]
+  [neohookean]
+    type=neohookean
+    params=100.0 0.3
+  [end]
+[end]
+
+
+
+[bcs]
+  [FixUx]
+    type=dirichlet
+    dof=ux
+    boundary=left
+    value=0.0
+  [end]
+  [FixUy]
+    type=dirichlet
+    dof=uy
+    boundary=bottom
+    value=0.0
+  [end]
+  [loadUx]
+    type=dirichlet
+    dof=uy
+    value=1.0*t
+    boundary=top
+  [end]
+[end]
+
+[timestepping]
+  type=be
+  dt=5.0e-3
+  endtime=8.0e-1
+  adaptive=false
+  optiters=3
+  dtmax=1.0e-1
+[end]
+
+[job]
+  type=transient
+  debug=dep
+[end]