! *****************************************************************************/ ! Core Loss Calculation based on the reference: ! V. Raulin, A. Radun and I. Husain, "Modeling of losses in switched reluctance ! machines," IEEE Trans. Ind. Appl., vol. 40, no. 6, pp. 1560-1569, Nov/Dec 2004. ! ! Programmed by Jenwel (jingwei.zhang.1989@ieee.org) ! Revised at 2017/Aug/31 ! *****************************************************************************/ !> \ingroup Solvers !> \{ !------------------------------------------------------------------------------ !> Calculate the time step core loss !> Based on Flux Density Field !------------------------------------------------------------------------------ SUBROUTINE CoreLossSolver( Model,Solver,dt,Transient ) !------------------------------------------------------------------------------ USE DefUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver !< Linear & nonlinear equation solver options TYPE(Model_t) :: Model !< All model information (mesh, materials, BCs, etc...) REAL(KIND=dp) :: dt !< Timestep size for time dependent simulations LOGICAL :: Transient !< Steady state or transient simulation !------------------------------------------------------------------------------ ! Local variables !------------------------------------------------------------------------------ INTEGER :: CurrentTimeStep INTEGER :: i ! Nsize is size of field INTEGER :: Nsize TYPE(Variable_t), POINTER :: Ph, Peddy, Pcoreloss, Bx, By REAL(KIND=dp), ALLOCATABLE :: Bfield(:), PrevBField(:), dBField(:) TYPE(Mesh_t), POINTER :: Mesh REAL(KIND=dp) :: Hc,Kh,Keddy REAL(KIND=dp) :: AxialL LOGICAL :: Found TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff TYPE(Element_t), POINTER :: Element INTEGER :: Nelem !------------------------------------------------------------------------------ SAVE Bfield,PrevBField,dBField ! Find Constants Hc = GetConstReal(Model % Constants, "Hc", Found) IF(.NOT.Found) CALL Fatal("CoreLossSolver", "Unable to find Hc") Kh = GetConstReal(Model % Constants, "Kh", Found) IF(.NOT.Found) CALL Fatal("CoreLossSolver", "Unable to find Kh") Keddy = GetConstReal(Model % Constants, "Keddy", Found) IF(.NOT.Found) CALL Fatal("CoreLossSolver", "Unable to find Keddy") AxialL = GetConstReal(Model % Constants, "AxialL", Found) IF(.NOT.Found) CALL Fatal("CoreLossSolver", "Unable to find AxialL") ! Fetch the target field, e.g. magnetic flux density !-------------------------------------------------------------------------- Bx => VariableGet( Solver % Mesh % Variables, 'Magnetic Flux Density 1') IF( .NOT. ASSOCIATED( Bx ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Magnetic Flux Density 1') ) END IF By => VariableGet( Solver % Mesh % Variables, 'Magnetic Flux Density 2') IF( .NOT. ASSOCIATED( Bx ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Magnetic Flux Density 2') ) END IF ! Please also define the Ph, Peddy, and Pcoreloss in the Solver section of "MagnetoDynamics2D" in sif file ! like: ! Exported Variable 1 = Ph ! Exported Variable 2 = Peddy ! Exported Variable 3 = Pcoreloss Mesh => GetMesh() Ph => VariableGet( Mesh % Variables,'Ph' ) IF( .NOT. ASSOCIATED( Ph ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Ph') ) END IF Peddy => VariableGet( Mesh % Variables,'Peddy' ) IF( .NOT. ASSOCIATED( Peddy ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Peddy') ) END IF Pcoreloss => VariableGet( Mesh % Variables,'Pcoreloss' ) IF( .NOT. ASSOCIATED( Pcoreloss ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Pcoreloss') ) END IF ! Nsize is the number of nodes Nsize = SIZE( Ph % Values) IF(.NOT. ALLOCATED(PrevBField) ) THEN ALLOCATE( BField(Nsize), PrevBField(Nsize), dBField(Nsize)) END IF !--------------------------------------------------------------------------------- DO i = 1, Nsize BField(i) = sqrt(Bx % Values(i) ** 2 + By % Values(i) ** 2) END DO CurrentTimeStep = GetTimestep() PRINT *, "Peddy % Values(i) =", Peddy % Values(i) PRINT *, "BField(1) =", BField(1) PRINT *, "PrevBField(1) =", PrevBField(1) PRINT *, "Nsize =", Nsize PRINT *, "Ph % Values(1) =", Ph % Values(1) IF (CurrentTimeStep<2) THEN Ph % Values = 0.0_dp Peddy % Values = 0.0_dp dBField = 0.0_dp PRINT *, "CurrentTimeStep =", CurrentTimeStep DO i = 1, Nsize PrevBField(i) = BField(i) END DO ELSE DO i = 1, Nsize dBField(i) = BField(i) - PrevBField(i) Ph % Values(i) = Hc * ABS(dBField(i)/dt) + Kh * ABS(BField(i)) * ABS(dBField(i)/dt) Peddy % Values(i) = Keddy * ((dBField(i)/dt) ** 2) Pcoreloss % Values(i) = Ph % Values(i) + Peddy % Values(i) PrevBField(i) = BField(i) END DO END IF PRINT *, "Peddy % Values(1) =", Peddy % Values(1) PRINT *, "Ph % Values(1) =", Ph % Values(1) PRINT *, "BField(1) =", BField(1) !----------------------------------------------------------------------------- ! Nelem is the number of active elements Nelem = GetNOFActive() PRINT *, "Nelem =", Nelem Element => GetActiveElement( 1 ) IntegStuff = GaussPoints( Element ) !------------------------------------------------------------------------------ END SUBROUTINE CoreLossSolver !------------------------------------------------------------------------------ !> \}