! *****************************************************************************/ ! 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_init0( Model,Solver,dt,Transient ) !------------------------------------------------------------------------------ USE DefUtils IMPLICIT NONE !------------------------------------------------------------------------------ !TYPE(Solver_t) :: Solver !< Linear & nonlinear equation solver options TYPE(Solver_t), TARGET :: Solver 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 !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ TYPE(Solver_t), POINTER :: Solvers(:), PSolver INTEGER :: n, m, mysolver, i TYPE(ValueList_t),POINTER :: DGSolverParams INTEGER, POINTER :: Active(:) LOGICAL :: Found !------------------------------------------------------------------------------ ! Ok, Dirty way of adding DG field for postprocessing ! Add a new solver with exactly the same active locations as the current one. PSolver => Solver DO mysolver=1,Model % NumberOfSolvers IF ( ASSOCIATED(PSolver,Model % Solvers(mysolver)) ) EXIT END DO n = Model % NumberOfSolvers DO i=1,Model % NumberOFEquations Active => ListGetIntegerArray(Model % Equations(i) % Values, & 'Active Solvers', Found) m = SIZE(Active) IF ( ANY(Active==mysolver) ) & CALL ListAddIntegerArray( Model % Equations(i) % Values, & 'Active Solvers', m+1, [Active, n+1] ) END DO ! Create DG solver structures on-the-fly without actually solving the matrix equations. ALLOCATE(Solvers(n+1)) Solvers(1:n) = Model % Solvers Solvers(n+1) % Values => ListAllocate() DGSolverParams => Solvers(n+1) % Values CALL ListAddLogical( DGSolverParams, 'Discontinuous Galerkin', .TRUE. ) Solvers(n+1) % DG = .TRUE. Solvers(n+1) % PROCEDURE = 0 Solvers(n+1) % ActiveElements => NULL() CALL ListAddString( DGSolverParams, 'Exec Solver', 'never' ) CALL ListAddLogical( DGSolverParams, 'No Matrix',.TRUE.) CALL ListAddLogical( DGSolverParams, 'Optimize Bandwidth',.FALSE.) CALL ListAddString( DGSolverParams, 'Equation', 'CoreLoss_Dummy' ) CALL ListAddString( DGSolverParams, 'Procedure', & 'CoreLoss CoreLossSolver_Dummy',.FALSE. ) CALL ListAddString( DGSolverParams, "Exported Variable 1 =", "Core Loss Phys e" ) CALL ListAddString( DGSolverParams, "Exported Variable 2 =", "Core Loss Peddy e") !CALL ListAddString( DGSolverParams, 'Variable','Core Loss 1e' ) !CALL ListAddString( DGSolverParams, 'Exported Variable 1','Core Loss 2e' ) CALL ListAddString( DGSolverParams, "Exported Variable 3 =", "CoreLoss_e" ) !PRINT *, "Initial Done", i DEALLOCATE(Model % Solvers) Model % Solvers => Solvers Model % NumberOfSolvers = n+1 END SUBROUTINE CoreLossSolver_init0 !-------------------------------------------------------------------------------- !------------------------------------------------------------------------------ SUBROUTINE CoreLossSolver_Dummy(Model,Solver,dt,Transient) !------------------------------------------------------------------------------ USE MagnetoDynamicsUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver TYPE(Model_t) :: Model REAL(KIND=dp) :: dt LOGICAL :: Transient !------------------------------------------------------------------------------ END SUBROUTINE CoreLossSolver_Dummy !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ SUBROUTINE CoreLossSolver( Model,Solver,dt,Transient ) !------------------------------------------------------------------------------ USE DefUtils IMPLICIT NONE !------------------------------------------------------------------------------ TYPE(Solver_t) :: Solver !< Linear & nonlinear equation solver options !TYPE(Solver_t), TARGET :: Solver 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, n TYPE(Variable_t), POINTER :: Ph, Peddy, Pcoreloss, Bx, By REAL(KIND=dp), ALLOCATABLE :: Bfield(:), PrevBField(:), dBField(:) REAL(KIND=dp), ALLOCATABLE :: Bfield_E(:), PrevBField_E(:), dBField_E(:) TYPE(Mesh_t), POINTER :: Mesh REAL(KIND=dp) :: Hc,Kh,Keddy REAL(KIND=dp) :: AxialL LOGICAL :: Found !------------------------------------------------------------------------------- TYPE VarPointer_t TYPE(Variable_t), POINTER :: Var END TYPE VarPointer_t TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff TYPE(Element_t), POINTER :: Element INTEGER :: Nelem !TYPE(VarPointer_t), POINTER :: Ph_E, Peddy_E TYPE(Variable_t), POINTER :: Ph_E, Peddy_E TYPE(Variable_t), POINTER :: Pcoreloss_E, Bx_E, By_E !------------------------------------------------------------------------------- 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 Bx_E => VariableGet( Solver % Mesh % Variables, 'Magnetic Flux Density E 1') IF( .NOT. ASSOCIATED( Bx_E ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Magnetic Flux Density E 1') ) END IF By_E => VariableGet( Solver % Mesh % Variables, 'Magnetic Flux Density E 2') IF( .NOT. ASSOCIATED( By_E ) ) THEN CALL Fatal('CoreLossSolver','Target field not present: '//TRIM('Magnetic Flux Density E 2') ) END IF Nelem = SIZE(Bx_E % Values) ! 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 IF(.NOT. ALLOCATED(PrevBField_E) ) THEN ALLOCATE( BField_E(Nelem), PrevBField_E(Nelem), dBField_E(Nelem)) 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 *, "Nelem =", Nelem !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 ) ! Check for Elemental (Discontinuous Galerkin) Field !------------------------------------------------------------------ !IF(.NOT. ALLOCATED( Ph_E ) ) ALLOCATE( Ph_E ) !IF(.NOT. ALLOCATED( Peddy_E ) ) ALLOCATE( Peddy_E ) !Ph_E => VariableGet( Mesh % Variables,'Variable Core Loss Phys e' ) !Ph_E % Var => VariableGet( Solver % Mesh % Variables,& ! 'Core Loss Phys e' ) !Peddy_E % Var => VariableGet( Solver % Mesh % Variables,& ! 'Core Loss Peddy e' ) Ph_E => VariableGet( Mesh % Variables,'Core Loss Phys e' ) IF(.NOT. ASSOCIATED( Ph_E ))THEN CALL Fatal('CoreLossSolver','Variable Core Loss Phys e does not exist') END IF !Peddy_E % Var => VariableGet( Solver % Mesh % Variables,& ! 'Peddy_e' ) !IF(.NOT. ASSOCIATED( Peddy_E % Var ) )THEN ! CALL Fatal('CoreLossSolver','Variable Core Loss Peddy e does not exist') !END IF !n = MAX(Solver % Mesh % MaxElementDOFs,Solver % Mesh % MaxElementNodes) !n=1 !PRINT *, "n =", SIZE(Ph_E % Values) !Ph_E % Var % Values(EPerm(Element % DGIndexes(1:n))) = FORCE(icomp,1:n) !------------------------------------------------------------------------------ END SUBROUTINE CoreLossSolver !------------------------------------------------------------------------------ !> \}