diff --git a/src/sm/Loads/hardmagneticboundarycondition.C b/src/sm/Loads/hardmagneticboundarycondition.C index aa6d21e75..5c945deb2 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.C +++ b/src/sm/Loads/hardmagneticboundarycondition.C @@ -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; @@ -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" ); @@ -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); } } } diff --git a/src/sm/Loads/hardmagneticboundarycondition.h b/src/sm/Loads/hardmagneticboundarycondition.h index 1ce8bc9e9..b1e09666b 100644 --- a/src/sm/Loads/hardmagneticboundarycondition.h +++ b/src/sm/Loads/hardmagneticboundarycondition.h @@ -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 //@} @@ -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 diff --git a/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.C b/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.C index bd1719eca..4908f18da 100644 --- a/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.C +++ b/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.C @@ -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; @@ -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; @@ -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( 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( 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 diff --git a/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.h b/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.h index e26b8400e..21bda0efb 100644 --- a/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.h +++ b/src/sm/Materials/MagnetoelasticMaterials/mooneyrivlinhardmagnetic.h @@ -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;