Skip to content

Commit

Permalink
Implementation of magnetoelasticity
Browse files Browse the repository at this point in the history
  • Loading branch information
@nitramkaroh
nitramkaroh committed Sep 8, 2024
1 parent 6d66ceb commit 0f92203
Show file tree
Hide file tree
Showing 9 changed files with 347 additions and 121 deletions.
94 changes: 72 additions & 22 deletions src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,8 @@
#include "classfactory.h"
#include "mathfem.h"
#include "tensor/tensor3.h"

#include "domain.h"
#include "function.h"


namespace oofem {
Expand All @@ -56,21 +57,14 @@ HardMagneticMooneyRivlinCompressibleMaterial :: give_FirstPKStressVector_Magneti
Tensor1_3d H(vH), Q, B;
Tensor2_3d F(vF), P;
//Q = (H+M)
Q(i_3) = H(i_3) + M(i_3);
double m_level = this->giveDomain()->giveFunction( m_ltf )->evaluateAtTime( tStep->giveIntrinsicTime() );
Q(i_3) = H(i_3) + m_level * M(i_3);
// compute the first Piola-Kirchhoff
auto [J, cofF] = F.compute_determinant_and_cofactor();
//
P( i_3, j_3 ) = C1 * this->compute_dI1_Cdev_dF(F)(i_3, j_3) + C2 * this->compute_dI2_Cdev_dF(F)(i_3, j_3) + this->compute_dVolumetricEnergy_dF(F)(i_3, j_3)
- ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 ));
//+ (mu_0 / (2*J*J)) * cofF(m_3, n_3) * Q(n_3) * cofF(m_3, k_3) * Q(k_3) * cofF(i_3, j_3);

/*
Tensor2_3d TT, TL;
TT( i_3, j_3 ) = - ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 ));
Tensor3_3d eps;
eps.be_Levi_Civita();
TL( i_3, j_3 ) = - ( mu_0 / J ) * eps(i_3, k_3, m_3) * F( k_3, l_3) * eps(j_3, l_3, n_3) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 ));
*/
- ( mu_0 / J ) * F.compute_tensor_cross_product()( i_3, j_3, m_3, n_3 ) * (cofF( m_3, q_3 ) * Q( q_3 ) * Q( n_3 ))
+ (mu_0 / (2*J*J)) * cofF(m_3, n_3) * Q(n_3) * cofF(m_3, k_3) * Q(k_3) * cofF(i_3, j_3);

B( i_3 ) = mu_0 / J * cofF( j_3, i_3 ) * cofF( j_3, k_3 ) * Q( k_3 );
auto vP = P.to_voigt_form();
Expand Down Expand Up @@ -110,24 +104,51 @@ HardMagneticMooneyRivlinCompressibleMaterial :: giveConstitutiveMatrices_dPdF_dB

//
//Q = (H+M)
Q(i_3) = H(i_3) + M(i_3);
double m_level = this->giveDomain()->giveFunction( m_ltf )->evaluateAtTime( tStep->giveIntrinsicTime() );
Q(i_3) = H(i_3) + m_level * M(i_3);
GQQ( i_3, j_3 ) = G( i_3, k_3 ) * Q( k_3 ) * Q( j_3 );
auto GQQcross = GQQ.compute_tensor_cross_product();
//

//
Tensor4_3d dPdF;
Tensor3_3d dPdH;
Tensor3_3d dBdF;
Tensor2_3d dBdH;
dPdF(i_3, j_3, k_3, l_3) = C1 * this->compute_d2I1_Cdev_dF2(F)(i_3, j_3, k_3, l_3)
+ C2 * this->compute_d2I2_Cdev_dF2(F)(i_3, j_3, k_3, l_3)
+ this->compute_d2VolumetricEnergy_dF2(F)(i_3, j_3, k_3, l_3)
+ C2 * this->compute_d2I2_Cdev_dF2(F)(i_3, j_3, k_3, l_3)
+ this->compute_d2VolumetricEnergy_dF2(F)(i_3, j_3, k_3, l_3)
- mu_0 / J * (GQQcross(i_3, j_3, k_3, l_3) + ( Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3) ) * ( Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 ) ))
+ mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3) ;
// + ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) )
// - ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3)
// + ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3)));
+ mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3)
+ ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) )
- ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3))
+ ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3));
//

Tensor4_3d A1, A2, A2t, A3;
//A1(i_3, j_3, k_3, l_3) = + ( mu_0 / ( 2 * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q( k_3 ) * Fcross( i_3, j_3, k_3, l_3 ) );
//A1(i_3, j_3, k_3, l_3) = - ( mu_0 / ( J * J * J ) ) * ( G( m_3, n_3 ) * Q( n_3 ) * G( m_3, k_3 ) * Q(k_3) * G(i_3, j_3) * G(k_3, l_3);
// A1(i_3, j_3, k_3, l_3) = + ( mu_0 / ( J * J ) ) * G(i_3, j_3) * (Fcross(k_3, l_3, m_3, n_3) * G(m_3,q_3) * Q(q_3) * Q(n_3)));

/*
A1(i_3, j_3, k_3, l_3) = - mu_0 / J * GQQcross(i_3, j_3, k_3, l_3);
Tensor3_3d FcQ1, FcQ2, FcQ3, FcQ4;
FcQ1(i_3, j_3, m_3) = Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3);
FcQ2(m_3, k_3, l_3) = Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 );
FcQ3(m_3, q_3, l_3) = Q( k_3 ) * Fcross( m_3, q_3, k_3, l_3 );
FcQ4(m_3, q_3, k_3) = Q( l_3 ) * Fcross( m_3, q_3, k_3, l_3 );
auto Fc = FloatMatrix(Fcross.to_voigt_form());
auto FcQ1m = FloatMatrix(FcQ1.to_voigt_form_3x9());
auto FcQ2m = FloatMatrix(FcQ2.to_voigt_form_3x9());
auto FcQ3m = FloatMatrix(FcQ3.to_voigt_form_3x9());
auto FcQ4m = FloatMatrix(FcQ4.to_voigt_form_3x9());
A2t(i_3, j_3, k_3, l_3) = FcQ1(i_3, j_3, m_3) * FcQ2(m_3, k_3, l_3);
A2(i_3, j_3, k_3, l_3) =- mu_0 / J * ( Fcross( i_3, j_3, m_3, n_3 ) * Q(n_3) ) * ( Fcross( m_3, q_3, k_3, l_3 ) * Q( q_3 ) );
A3(i_3, j_3, k_3, l_3) = + mu_0/J/J * Fcross( i_3, j_3, m_3, n_3 ) * GQQ(m_3, n_3) * G(k_3, l_3) ;
//
auto A1v = FloatMatrix(A1.to_voigt_form());
auto A2v = FloatMatrix(A2.to_voigt_form());
auto A3v = FloatMatrix(A3.to_voigt_form());
*/
dBdH( i_3, k_3 ) = -mu_0 / J * G( j_3, i_3 ) * G( j_3, k_3 );
//
dBdF( i_3, k_3, l_3 ) = mu_0 / J * (G( j_3, m_3 ) * Q( m_3 )* Fcross( j_3, i_3, k_3, l_3 )) + G( j_3, i_3 ) * (Q( m_3 ) * Fcross( j_3, m_3, k_3, l_3 )) - mu_0 / ( J * J ) * ( G( j_3, i_3 ) * G( j_3, m_3 ) * Q( m_3 ) * G( k_3, l_3 ) );
Expand All @@ -139,7 +160,7 @@ HardMagneticMooneyRivlinCompressibleMaterial :: giveConstitutiveMatrices_dPdF_dB

auto [K1, K2, K3, K4] = this->compute_stiff_num(vF, vH, gp, tStep);


auto D1 = FloatMatrix(dPdF.to_voigt_form());
auto D2 = FloatMatrix(dPdH.to_voigt_form_9x3());
auto D3 = FloatMatrix(dBdF.to_voigt_form_3x9());
Expand Down Expand Up @@ -223,11 +244,40 @@ HardMagneticMooneyRivlinCompressibleMaterial::initializeFrom(InputRecord &ir)
IR_GIVE_FIELD(ir, C2, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c2);
FloatArray m;
IR_GIVE_FIELD(ir, m, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_magnetization);
IR_GIVE_FIELD(ir, m_ltf, _IFT_HardMagneticMooneyRivlinCompressibleMaterial_m_ltf);
if(m.giveSize() != 3){
OOFEM_ERROR("Magnetization has to be vector of size 3");
} else {
M = Tensor1_3d(FloatArrayF<3>(m));
}
}




int
HardMagneticMooneyRivlinCompressibleMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{

MagnetoElasticMaterialStatus *status = static_cast< MagnetoElasticMaterialStatus * >( this->giveStatus(gp) );
if ( type == IST_DeformationGradientTensor ) {
answer = status->giveFVector();
return 1;
} else if ( type == IST_FirstPKStressTensor ) {
answer = status->givePVector();
return 1;
} else if ( type == IST_MagneticInductionVector ) {
answer = status->giveBVector();
return 1;
} else if ( type == IST_MagneticFieldVector ) {
answer = status->giveHVector();
return 1;
} else {
return Material::giveIPValue(answer, gp, type, tStep);
}
}




} // end namespace oofem
6 changes: 6 additions & 0 deletions src/mpm/MagnetoElasticity/Materials/hardmagneticmooneyrivlincompressiblematerial.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
//@{
#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_Name "hardmagneticmooneyrivlincompressiblemat"
#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_magnetization "magnetization"
#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_m_ltf "m_ltf"
#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c1 "c1"
#define _IFT_HardMagneticMooneyRivlinCompressibleMaterial_c2 "c2"

Expand All @@ -65,6 +66,7 @@ class HardMagneticMooneyRivlinCompressibleMaterial : public MagnetoElasticMateri
double C1;
double C2;
double mu_0 = 1.25663706143e-6;
int m_ltf = 0;



Expand All @@ -86,6 +88,10 @@ class HardMagneticMooneyRivlinCompressibleMaterial : public MagnetoElasticMateri

const char *giveInputRecordName() const override { return _IFT_HardMagneticMooneyRivlinCompressibleMaterial_Name; }
const char *giveClassName() const override { return "HardMagneticMooneyRivlinCompressibleMaterial"; }

int giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep);


};
} // end namespace oofem
#endif
82 changes: 81 additions & 1 deletion src/mpm/magnetoelasticelement.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include "intarray.h"
#include "classfactory.h"
#include "gaussintegrationrule.h"
#include "crosssection.h"

#include "fei2dquadquad.h"
#include "fei2dquadlin.h"
Expand Down Expand Up @@ -95,7 +96,8 @@ class MagnetoElasticElement : public MPElement {
} else {
OOFEM_ERROR("Unknown characteristic vector type");
}
}
}

};


Expand Down Expand Up @@ -185,5 +187,83 @@ REGISTER_Element(MagnetoElasticQuad_ql)



/**
* @brief 2D Magneto elatic element with linear interpolation of displacements and magnetic potential
*
*/
class MagnetoElasticQuad_ll : public MagnetoElasticElement {
protected:
const static FEInterpolation & phiInterpol;
const static FEInterpolation & uInterpol;
const static Variable& phi;
const static Variable& u;

public:
MagnetoElasticQuad_ll(int n, Domain* d): MagnetoElasticElement(n,d)
{
numberOfDofMans = 4;
numberOfGaussPoints = 4;
this->computeGaussPoints();
}

int getNumberOfSurfaceDOFs() const override {return 0;}
int getNumberOfEdgeDOFs() const override {return 0;}
void getSurfaceLocalCodeNumbers(oofem::IntArray&, oofem::Variable::VariableQuantity) const override {;}

void getDofManLocalCodeNumbers (IntArray& answer, const Variable::VariableQuantity q, int num ) const override {
/* dof ordering: u1 v1 phi1 u2 v2 phi2 u3 v3 phi3 u4 v4 phi4 u5 v5 u6 v6 */
if (q == Variable::VariableQuantity::Displacement) {
//answer={1,2, 4,5, 7,8, 10,11, 13,14 15,16, 17,18 19,20};
int o = (num-1)*3+1;
answer={o, o+1};
} else if (q == Variable::VariableQuantity::MagneticPotential) {
answer={num*3};
}
}
void getInternalDofManLocalCodeNumbers (IntArray& answer, const Variable::VariableQuantity q, int num ) const override {
answer={};
}

void giveDofManDofIDMask(int inode, IntArray &answer) const override {
answer = {1,2,55};
}
int giveNumberOfDofs() override { return 12; }
const char *giveInputRecordName() const override {return "MagnetoElasticQuad_ll";}
const char *giveClassName() const override { return "MagnetoElasticQuad_ll"; }


const FEInterpolation& getGeometryInterpolation() const override {return this->uInterpol;}

Element_Geometry_Type giveGeometryType() const override {
return EGT_quad_1;
}
void getEdgeLocalCodeNumbers(IntArray& answer, const Variable::VariableQuantity q) const override {}


private:
virtual int giveNumberOfUDofs() const override {return 8;}
virtual int giveNumberOfPhiDofs() const override {return 4;}
virtual const Variable& getU() const override {return u;}
virtual const Variable& getPhi() const override {return phi;}
void computeGaussPoints() override {
if ( integrationRulesArray.size() == 0 ) {
integrationRulesArray.resize( 1 );
integrationRulesArray [ 0 ] = std::make_unique<GaussIntegrationRule>(1, this);
integrationRulesArray [ 0 ]->SetUpPointsOnSquare(numberOfGaussPoints, _PlaneStrain);
}
}
};

const FEInterpolation & MagnetoElasticQuad_ll :: uInterpol = FEI2dQuadLin(1,2);
const FEInterpolation & MagnetoElasticQuad_ll :: phiInterpol = FEI2dQuadLin(1,2);
const Variable& MagnetoElasticQuad_ll::phi = Variable( MagnetoElasticQuad_ll::phiInterpol, Variable::VariableQuantity::MagneticPotential, Variable::VariableType::scalar, 1, NULL, {55});
const Variable& MagnetoElasticQuad_ll::u = Variable(MagnetoElasticQuad_ll::uInterpol, Variable::VariableQuantity::Displacement, Variable::VariableType::vector, 2, NULL, {1,2});

#define _IFT_MagnetoElasticQuad_ll_Name "magnetoelasticquad_ll"
REGISTER_Element(MagnetoElasticQuad_ll)





} // end namespace oofem
2 changes: 1 addition & 1 deletion src/mpm/termlibrary4.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,7 +229,7 @@ MagnetoElasticity_GradGrad_Term :: evaluate_lin (FloatMatrix& answer, MPElement&
// Grad contains deformation gradient and magnetic field
auto size = this->computeGradientFields(vGrad, vB, cell, gp->giveNaturalCoordinates(), gp->giveMaterialMode(), tstep);
FloatMatrix D1,D2,D3,D4;
this->compute_lin_num(D1,D2,D3,D4, cell, gp, tstep);
// this->compute_lin_num(D1,D2,D3,D4, cell, gp, tstep);
//
auto cs = cell.giveCrossSection();
auto mcs = dynamic_cast<MagnetoElasticCrossSection *> (cs);
Expand Down
5 changes: 4 additions & 1 deletion src/oofemlib/internalstatetype.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,10 @@ namespace oofem {
ENUM_ITEM_WITH_VALUE(IST_VolumeFraction, 146) \
ENUM_ITEM_WITH_VALUE(IST_X_LCS, 147) /*Unit vector in local coordinate system in the x direction (usable for diagrams of internal forces for VTK export)*/ \
ENUM_ITEM_WITH_VALUE(IST_Y_LCS, 148) \
ENUM_ITEM_WITH_VALUE(IST_Z_LCS, 149)
ENUM_ITEM_WITH_VALUE(IST_Z_LCS, 149) \
ENUM_ITEM_WITH_VALUE(IST_MagneticPotential, 150) \
ENUM_ITEM_WITH_VALUE(IST_MagneticFieldVector, 151) \
ENUM_ITEM_WITH_VALUE(IST_MagneticInductionVector, 152)


/**
Expand Down
Loading

0 comments on commit 0f92203

Please sign in to comment.