Skip to content

Commit

Permalink
Multiple changes in everything magnetic
Browse files Browse the repository at this point in the history
  • Loading branch information
Ondřej Faltus committed Oct 3, 2024
1 parent 5639f37 commit 8039c0d
Show file tree
Hide file tree
Showing 4 changed files with 184 additions and 5 deletions.
21 changes: 17 additions & 4 deletions src/sm/Loads/hardmagneticboundarycondition.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,13 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );
{
ActiveBoundaryCondition::initializeFrom( ir );

FloatArray b;
FloatArray b, m;
IR_GIVE_FIELD( ir, b, _IFT_HardMagneticBoundaryCondition_b_ext );
b_ext = Tensor1_3d( FloatArrayF< 3 >(b) );
b_ext = Tensor1_3d( FloatArrayF< 3 >( b ) );
IR_GIVE_FIELD( ir, m, _IFT_HardMagneticBoundaryCondition_mjump );
m_jump = Tensor1_3d( FloatArrayF<3>( m ) );


IR_GIVE_FIELD( ir, ltf_index, _IFT_HardMagneticBoundaryCondition_ltf );

mu0 = BASE_VACUUM_PERMEABILITY_MU_0;
Expand Down Expand Up @@ -96,7 +100,16 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );
double pert = 1.e-6;
this->computeNumericalTangentFromElement( Ke_num, e, boundary, tStep, pert );

if ( pos == 3 ) {
//symmetry check
FloatMatrix KeT, Ke_numT, errorKe, errorKe_num;
KeT.beTranspositionOf( Ke );
Ke_numT.beTranspositionOf( Ke_num );
errorKe = KeT;
errorKe_num = Ke_numT;
errorKe.add( -1. * Ke );
errorKe_num.add( -1. * Ke_num );

if ( pos < 0 ) {

OOFEM_LOG_INFO( "Debugging surface %i\n", boundary );
OOFEM_LOG_INFO( "-------Analytical stiffness-----------\n" );
Expand Down Expand Up @@ -312,7 +325,7 @@ REGISTER_BoundaryCondition( HardMagneticBoundaryCondition );

for ( int j = 1; j <= 12; j++ ) {

answer.at( j, i ) = ( perturbedForceFront.at( j ) - perturbedForceBack.at( j ) ) / (2 * perturb);
answer.at( i, j ) = ( perturbedForceFront.at( j ) - perturbedForceBack.at( j ) ) / (2 * perturb);
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion src/sm/Loads/hardmagneticboundarycondition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@
#define _IFT_HardMagneticBoundaryCondition_Name "hardmagneticboundarycondition"
#define _IFT_HardMagneticBoundaryCondition_mu_0 "mu_0"
#define _IFT_HardMagneticBoundaryCondition_b_ext "b_ext"
#define _IFT_HardMagneticBoundaryCondition_mjump "mjump"
#define _IFT_HardMagneticBoundaryCondition_ltf "loadtimefunction" //load time function for applied field
//@}

Expand All @@ -61,7 +62,7 @@ class HardMagneticBoundaryCondition : public ActiveBoundaryCondition {
protected:

double mu0; //vacuum permeability
Tensor1_3d b_ext; //external magnetic field in free space
Tensor1_3d b_ext, m_jump; //external magnetic field in free space, jump in magnetization
//FloatArrayF<9> sigma_star; //precomputed free space stress
Tensor2_3d maxwell_stress;
int ltf_index; //index of load time function for applied load
Expand Down
157 changes: 157 additions & 0 deletions src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,12 @@ MooneyRivlinHardMagnetic::giveFirstPKStressVector_3d( const FloatArrayF<9> &vF,
/*case 5:
vP_me = computeFirstPKStressVector_3d_multiphysics( vF, gp, tStep );
break;*/
case 6:
vP_me = computeFirstPKStressVector_3d_universal( vF, gp, tStep );
break;
case 7:
vP_me = computeFirstPKStressVector_3d_universal_nobb( vF, gp, tStep );
break;
default:
OOFEM_ERROR( "Unknown hard magnetic material mode %i", materialMode );
break;
Expand Down Expand Up @@ -114,6 +120,12 @@ MooneyRivlinHardMagnetic::give3dMaterialStiffnessMatrix_dPdF( MatResponseMode ma
/* case 5:
vD_me = compute3dMaterialStiffnessMatrix_dPdF_multiphysics( matResponseMode, gp, tStep );
break;*/
case 6:
vD_me = compute3dMaterialStiffnessMatrix_dPdF_universal( matResponseMode, gp, tStep );
break;
case 7:
vD_me = compute3dMaterialStiffnessMatrix_dPdF_universal_nobb( matResponseMode, gp, tStep );
break;
default:
OOFEM_ERROR( "Unknown hard magnetic material mode %i", materialMode );
break;
Expand Down Expand Up @@ -350,6 +362,151 @@ FloatMatrixF<9, 9> MooneyRivlinHardMagnetic::compute3dMaterialStiffnessMatrix_dP
return vD_me;
}

//--------------------------------------------------------------------------- Universal -------------------------------------------------------------

FloatArrayF<9> MooneyRivlinHardMagnetic::computeFirstPKStressVector_3d_universal( const FloatArrayF<9> &vF, GaussPoint *gp, TimeStep *tStep ) const
{

double load_level = this->giveDomain()->giveFunction( ltf_index )->evaluateAtTime( tStep->giveIntrinsicTime() );
FloatArrayF<3> B_app_at_time = load_level * B_app;

Tensor2_3d F( vF ), P_me, delta( 1., 0., 0., 0., 1., 0., 0., 0., 1. ), Finv;
Finv = F.compute_inverse();

Tensor1_3d Bapp( B_app_at_time ), M( (1./mu_0) * B_res ), Bappref;
auto [J, cofF] = F.compute_determinant_and_cofactor();

if ( !referenceB ) {
Bappref( j_3 ) = Bapp( k_3 ) * cofF( k_3, j_3 );
} else {
Bappref = Bapp;
}

P_me( i_3, j_3 ) =
+ ( 1. / J ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
- ( 1. / J ) * F( i_3, m_3 ) * Bappref( m_3 ) * M( j_3 )
- ( 1. / J ) * F( i_3, m_3 ) * M( m_3 ) * Bappref( j_3 )
- ( 1. / ( 2. * mu_0 * J ) ) * F( m_3, n_3 ) * Bappref( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
+ ( 1. / (mu_0 * J) ) * F( i_3, m_3 ) * Bappref( m_3 ) * Bappref( j_3 );

auto vP_me = P_me.to_voigt_form();

return vP_me;
}

FloatMatrixF<9, 9> MooneyRivlinHardMagnetic::compute3dMaterialStiffnessMatrix_dPdF_universal( MatResponseMode matResponseMode, GaussPoint *gp, TimeStep *tStep ) const
{
StructuralMaterialStatus *status = static_cast<StructuralMaterialStatus *>( this->giveStatus( gp ) );
FloatArrayF<9> vF( status->giveTempFVector() );


double load_level = this->giveDomain()->giveFunction( ltf_index )->evaluateAtTime( tStep->giveIntrinsicTime() );
FloatArrayF<3> B_app_at_time = load_level * B_app;

Tensor2_3d F( vF );
Tensor2_3d Finv = F.compute_inverse();
Tensor4_3d D_me;

Tensor2_3d delta(1., 0., 0., 0., 1., 0., 0., 0., 1. );

Tensor1_3d Bapp ( B_app_at_time ), M( (1./mu_0)*B_res ), Bappref;
auto [J, cofF] = F.compute_determinant_and_cofactor();

if ( !referenceB ) {
Bappref( j_3 ) = Bapp( k_3 ) * cofF( k_3, j_3 );
} else {
Bappref = Bapp;
}

D_me( i_3, j_3, k_3, l_3 ) = -( 1. / ( J * J ) ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 ) * cofF( k_3, l_3 )
+ ( 1. / J ) * M( l_3 ) * F( k_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
+ ( 1. / J ) * F( k_3, n_3 ) * M( n_3 ) * Bappref( l_3 ) * Finv( j_3, i_3 )
- ( 1. / J ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, k_3 ) * Finv( l_3, i_3 )
+ ( 1. / ( J * J ) ) * F( i_3, m_3 ) * Bappref( m_3 ) * M( j_3 ) * cofF( k_3, l_3 )
- ( 1. / J ) * delta( i_3, k_3 ) * ( Bappref( l_3 ) * M( j_3 ) )
+ ( 1. / ( J * J ) ) * F( i_3, m_3 ) * M( m_3 ) * Bappref( j_3 ) * cofF( k_3, l_3 )
- ( 1. / J ) * delta( i_3, k_3 ) * ( M( l_3 ) * Bappref( j_3 ) )
+ ( 1. / ( 2. * mu_0 * J * J ) ) * F( m_3, n_3 ) * Bappref( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 ) * cofF( k_3, l_3 )
- ( 1. / ( 2. * mu_0 * J ) ) * Bappref( l_3 ) * F( k_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
- ( 1. / ( 2. * mu_0 * J ) ) * F( k_3, n_3 ) * Bappref( n_3 ) * Bappref( l_3 ) * Finv( j_3, i_3 )
+ ( 1. / ( 2. * mu_0 * J ) ) * F( m_3, n_3 ) * Bappref( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, k_3 ) * Finv( l_3, i_3 )
- ( 1. / ( mu_0 * J * J ) ) * F( i_3, m_3 ) * Bappref( m_3 ) * Bappref( j_3 ) * cofF( k_3, l_3 )
+ ( 1. / ( mu_0 * J ) ) * delta( i_3, k_3 ) * ( Bappref( l_3 ) * Bappref( j_3 ) );


auto vD_me = D_me.to_voigt_form();

return vD_me;
}

//--------------------------------------------------------------------------- Universal- No BB -------------------------------------------------------------

FloatArrayF<9> MooneyRivlinHardMagnetic::computeFirstPKStressVector_3d_universal_nobb( const FloatArrayF<9> &vF, GaussPoint *gp, TimeStep *tStep ) const
{

double load_level = this->giveDomain()->giveFunction( ltf_index )->evaluateAtTime( tStep->giveIntrinsicTime() );
FloatArrayF<3> B_app_at_time = load_level * B_app;

Tensor2_3d F( vF ), P_me, delta( 1., 0., 0., 0., 1., 0., 0., 0., 1. ), Finv;
Finv = F.compute_inverse();

Tensor1_3d Bapp( B_app_at_time ), M( (1./mu_0) * B_res ), Bappref;
auto [J, cofF] = F.compute_determinant_and_cofactor();

if ( !referenceB ) {
Bappref( j_3 ) = Bapp( k_3 ) * cofF( k_3, j_3 );
} else {
Bappref = Bapp;
}

P_me( i_3, j_3 ) = +( 1. / J ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
- ( 1. / J ) * F( i_3, m_3 ) * Bappref( m_3 ) * M( j_3 )
- ( 1. / J ) * F( i_3, m_3 ) * M( m_3 ) * Bappref( j_3 );

auto vP_me = P_me.to_voigt_form();

return vP_me;
}

FloatMatrixF<9, 9> MooneyRivlinHardMagnetic::compute3dMaterialStiffnessMatrix_dPdF_universal_nobb( MatResponseMode matResponseMode, GaussPoint *gp, TimeStep *tStep ) const
{
StructuralMaterialStatus *status = static_cast<StructuralMaterialStatus *>( this->giveStatus( gp ) );
FloatArrayF<9> vF( status->giveTempFVector() );


double load_level = this->giveDomain()->giveFunction( ltf_index )->evaluateAtTime( tStep->giveIntrinsicTime() );
FloatArrayF<3> B_app_at_time = load_level * B_app;

Tensor2_3d F( vF );
Tensor2_3d Finv = F.compute_inverse();
Tensor4_3d D_me;

Tensor2_3d delta(1., 0., 0., 0., 1., 0., 0., 0., 1. );

Tensor1_3d Bapp ( B_app_at_time ), M( (1./mu_0)*B_res ), Bappref;
auto [J, cofF] = F.compute_determinant_and_cofactor();

if ( !referenceB ) {
Bappref( j_3 ) = Bapp( k_3 ) * cofF( k_3, j_3 );
} else {
Bappref = Bapp;
}

D_me( i_3, j_3, k_3, l_3 ) = -( 1. / ( J * J ) ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 ) * cofF( k_3, l_3 )
+ ( 1. / J ) * M( l_3 ) * F( k_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, i_3 )
+ ( 1. / J ) * F( k_3, n_3 ) * M( n_3 ) * Bappref( l_3 ) * Finv( j_3, i_3 )
- ( 1. / J ) * F( m_3, n_3 ) * M( n_3 ) * F( m_3, p_3 ) * Bappref( p_3 ) * Finv( j_3, k_3 ) * Finv( l_3, i_3 )
+ ( 1. / ( J * J ) ) * F( i_3, m_3 ) * Bappref( m_3 ) * M( j_3 ) * cofF( k_3, l_3 )
- ( 1. / J ) * delta( i_3, k_3 ) * ( Bappref( l_3 ) * M( j_3 ) )
+ ( 1. / ( J * J ) ) * F( i_3, m_3 ) * M( m_3 ) * Bappref( j_3 ) * cofF( k_3, l_3 )
- ( 1. / J ) * delta( i_3, k_3 ) * ( M( l_3 ) * Bappref( j_3 ) );


auto vD_me = D_me.to_voigt_form();

return vD_me;
}

//--------------------------------------------------------------------------- Multiphysics ----------------------------------------------------------

//FloatArrayF<9> MooneyRivlinHardMagnetic::computeFirstPKStressVector_3d_multiphysics( const FloatArrayF<9> &vF, GaussPoint *gp, TimeStep *tStep ) const
Expand Down
8 changes: 8 additions & 0 deletions src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,14 @@ class MooneyRivlinHardMagnetic : public MooneyRivlinCompressibleMaterial//, publ
FloatMatrixF< 9, 9 > compute3dMaterialStiffnessMatrix_dPdF_ogdenpullback(MatResponseMode matResponseMode,
GaussPoint *gp, TimeStep *tStep ) const;

FloatArrayF< 9 > computeFirstPKStressVector_3d_universal(const FloatArrayF< 9 > &vF, GaussPoint *gp, TimeStep *tStep) const;
FloatMatrixF< 9, 9 > compute3dMaterialStiffnessMatrix_dPdF_universal(MatResponseMode matResponseMode,
GaussPoint *gp, TimeStep *tStep ) const;

FloatArrayF< 9 > computeFirstPKStressVector_3d_universal_nobb(const FloatArrayF< 9 > &vF, GaussPoint *gp, TimeStep *tStep) const;
FloatMatrixF< 9, 9 > compute3dMaterialStiffnessMatrix_dPdF_universal_nobb(MatResponseMode matResponseMode,
GaussPoint *gp, TimeStep *tStep ) const;

//FloatArrayF< 9 > computeFirstPKStressVector_3d_multiphysics(const FloatArrayF< 9 > &vF, GaussPoint *gp, TimeStep *tStep) const;
//FloatMatrixF<9, 9> compute3dMaterialStiffnessMatrix_dPdF_multiphysics( MatResponseMode matResponseMode, GaussPoint *gp, TimeStep *tStep ) const;
//FloatArrayF<3> computeMagneticInduction_3d_multiphysics( const FloatArrayF<3> &vH, GaussPoint *gp, TimeStep *tStep ) const;
Expand Down

0 comments on commit 8039c0d

Please sign in to comment.