Skip to content

Commit

Permalink
Merge branch 'tensor_reluctivity' into devel
Browse files Browse the repository at this point in the history 
Updates to support reluctivity as a tensor-valued property + some fixes
  • Loading branch information
@mmalinen
mmalinen committed Oct 29, 2020
2 parents 311f84e + 1552543 commit bc9e064
Show file tree
Hide file tree
Showing 3 changed files with 130 additions and 102 deletions.
150 changes: 93 additions & 57 deletions fem/src/modules/MagnetoDynamics/CalcFields.F90
View file
Expand Up @@ -518,24 +518,29 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
! a background element of type 827):
!------------------------------------------------------------------------------
REAL(KIND=dp) :: WBasis(54,3), RotWBasis(54,3), Basis(27), dBasisdx(27,3)
REAL(KIND=dp) :: SOL(2,81), PSOL(81), ElPotSol(1,27), R(27), C(27)
REAL(KIND=dp) :: SOL(2,81), PSOL(81), ElPotSol(1,27), C(27)
REAL(KIND=dp) :: Wbase(27), alpha(27), NF_ip(27,3)
REAL(KIND=dp) :: PR(27), omega_velo(3,27), lorentz_velo(3,27)
COMPLEX(KIND=dp) :: Magnetization(3,27), BodyForceCurrDens(3,27)
COMPLEX(KIND=dp) :: Magnetization(3,27), BodyForceCurrDens(3,27)
COMPLEX(KIND=dp) :: R_Z(27)
!------------------------------------------------------------------------------
REAL(KIND=dp) :: s,u,v,w, Norm
REAL(KIND=dp) :: B(2,3), E(2,3), JatIP(2,3), VP_ip(2,3), JXBatIP(2,3), CC_J(2,3), B2
REAL(KIND=dp) :: detJ, C_ip, R_ip, PR_ip, ST(3,3), Omega, ThinLinePower, Power, Energy(2), w_dens, R_t_ip(3,3)
REAL(KIND=dp) :: B(2,3), E(2,3), JatIP(2,3), VP_ip(2,3), JXBatIP(2,3), CC_J(2,3), HdotB
REAL(KIND=dp) :: detJ, C_ip, PR_ip, ST(3,3), Omega, ThinLinePower, Power, Energy(2), w_dens
REAL(KIND=dp) :: Freq, FreqPower, FieldPower, LossCoeff, ValAtIP
REAL(KIND=dp) :: Freq2, FreqPower2, FieldPower2, LossCoeff2
REAL(KIND=dp) :: ComponentLoss(2,2), rot_velo(3), angular_velo(3)
REAL(KIND=dp) :: Coeff, Coeff2, TotalLoss(3), LumpedForce(3), localAlpha, localV(2), nofturns, coilthickness
REAL(KIND=dp) :: Flux(2), AverageFluxDensity(2), Area, N_j, wvec(3), PosCoord(3), TorqueDeprecated(3)
REAL(KIND=dp) :: R_ip, mu_r

COMPLEX(KIND=dp) :: MG_ip(3), BodyForceCurrDens_ip(3)
COMPLEX(KIND=dp) :: CST(3,3)
COMPLEX(KIND=dp) :: CMat_ip(3,3)
COMPLEX(KIND=dp) :: imag_value, Zs
COMPLEX(KIND=dp), ALLOCATABLE :: Tcoef(:,:,:)
COMPLEX(KIND=dp), POINTER, SAVE :: Reluct_Z(:,:,:) => NULL()
COMPLEX(KIND=dp) :: R_ip_Z, Nu(3,3)

INTEGER, PARAMETER :: ind1(6) = [1,2,3,1,2,1]
INTEGER, PARAMETER :: ind2(6) = [1,2,3,2,3,3]
Expand All @@ -553,11 +558,14 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
CHARACTER(LEN=MAX_NAME_LEN) :: Pname, CoilType, ElectricPotName, LossFile, CurrPathPotName

TYPE(ValueList_t), POINTER :: Material, BC, BodyForce, BodyParams, SolverParams

LOGICAL :: Found, FoundMagnetization, stat, Cubic, LossEstimation, &
CalcFluxLogical, CoilBody, PreComputedElectricPot, ImposeCircuitCurrent, &
ItoJCoeffFound, ImposeBodyForceCurrent, HasVelocity, HasAngularVelocity, &
HasLorenzVelocity, HaveAirGap, UseElementalNF, HasTensorReluctivity, &
ImposeBodyForcePotential, JouleHeatingFromCurrent, HasZirka
LOGICAL :: PiolaVersion, ElementalFields, NodalFields, RealField, SecondOrder
LOGICAL :: CSymmetry, HBCurve, LorentzConductivity, HasThinLines=.FALSE.

TYPE(GaussIntegrationPoints_t) :: IP
TYPE(Nodes_t), SAVE :: Nodes
Expand All @@ -577,15 +585,12 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
REAL(KIND=dp), ALLOCATABLE :: ThinLineCrossect(:),ThinLineCond(:)

REAL(KIND=DP), POINTER :: Cwrk(:,:,:)=>NULL(), Cwrk_im(:,:,:)=>NULL()
COMPLEX(KIND=dp), ALLOCATABLE :: Tcoef(:,:,:)
REAL(KIND=dp), POINTER :: R_t(:,:,:)

LOGICAL :: PiolaVersion, ElementalFields, NodalFields, RealField, SecondOrder
REAL(KIND=dp) :: ItoJCoeff, CircuitCurrent, CircEqVoltageFactor
TYPE(ValueList_t), POINTER :: CompParams
REAL(KIND=dp) :: DetF, F(3,3), G(3,3), GT(3,3)
REAL(KIND=dp), ALLOCATABLE :: EBasis(:,:), CurlEBasis(:,:)
LOGICAL :: CSymmetry, HBCurve, LorentzConductivity, HasThinLines=.FALSE.

REAL(KIND=dp) :: xcoord, grads_coeff, val
TYPE(ValueListEntry_t), POINTER :: HBLst
REAL(KIND=dp) :: HarmPowerCoeff
Expand Down Expand Up @@ -773,7 +778,7 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
IF ( ASSOCIATED(EL_ML2) ) ElementalFields=.TRUE.

n = Mesh % MaxElementDOFs
ALLOCATE( MASS(n,n), FORCE(n,DOFs), Tcoef(3,3,n), RotM(3,3,n), Pivot(n), R_t(3,3,n))
ALLOCATE( MASS(n,n), FORCE(n,DOFs), Tcoef(3,3,n), RotM(3,3,n), Pivot(n))

SOL = 0._dp; PSOL=0._dp

Expand Down Expand Up @@ -812,7 +817,7 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
END IF


C = 0._dp; R=0._dp; PR=0._dp
C = 0._dp; PR=0._dp
Magnetization = 0._dp

Power = 0._dp; Energy = 0._dp
Expand Down Expand Up @@ -976,7 +981,7 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)


!---------------------------------------------------------------------------------------------

R_Z = CMPLX(0.0_dp, 0.0_dp, kind=dp)
HasTensorReluctivity = .FALSE.
CALL GetConstRealArray( Material, HB, 'H-B curve', Found )
IF ( ASSOCIATED(HB) ) THEN
Expand All @@ -1001,17 +1006,22 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
END IF
END IF
ELSE
CALL GetReluctivity(Material,R_t,n,HasTensorReluctivity)
!
! Seek reluctivity as complex-valued: A given reluctivity can be a tensor
!
CALL GetReluctivity(Material,Reluct_Z,n,HasTensorReluctivity)
IF (HasTensorReluctivity) THEN
IF (SIZE(R_t,1)==1 .AND. SIZE(R_t,2)==1) THEN
l = MIN(SIZE(R), SIZE(R_t,3))
R(1:l) = R_t(1,1,1:l)
IF (SIZE(Reluct_Z,1)==1 .AND. SIZE(Reluct_Z,2)==1) THEN
l = MIN(SIZE(R_Z), SIZE(Reluct_Z,3))
R_Z(1:l) = Reluct_Z(1,1,1:l)
HasTensorReluctivity = .FALSE.
ELSE
R = 0.0d0
END IF
R_Z = CMPLX(0.0_dp, 0.0_dp, kind=dp)
END IF
ELSE
CALL GetReluctivity(Material,R,n)
! Seek via a given permeability: In this case the reluctivity will be
! a complex scalar:
CALL GetReluctivity(Material,R_Z,n)
END IF
END IF

Expand Down Expand Up @@ -1338,29 +1348,53 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
END SELECT
END IF


Nu = CMPLX(0.0d0, 0.0d0, kind=dp)
IF ( ASSOCIATED(HB) ) THEN
Babs=SQRT(SUM(B(1,:)**2))
IF (RealField) THEN
Babs=SQRT(SUM(B(1,:)**2))
ELSE
Babs = SQRT(SUM(B(1,:)**2 + B(2,:)**2))
END IF
Babs = MAX(Babs, 1.d-8)
R_ip = InterpolateCurve(HBBval,HBHval,Babs,HBCval)/Babs
w_dens = IntegrateCurve(HBBval,HBHval,HBCval,0._dp,Babs)
DO k=1,3
Nu(k,k) = CMPLX(R_ip, 0.0d0, kind=dp)
END DO
ELSE
R_ip = SUM( Basis(1:n)*R(1:n) )
IF (HasTensorReluctivity) THEN
IF (SIZE(R_t,2) == 1) THEN
R_t_ip = 0.0d0
DO k = 1,3
R_t_ip(k,k) = SUM(Basis(1:n)*R_t(k,1,1:n))
IF (SIZE(Reluct_Z,2) == 1) THEN
DO k = 1, MIN(3, SIZE(Reluct_Z,1))
Nu(k,k) = SUM(Basis(1:n)*Reluct_Z(k,1,1:n))
END DO
ELSE
DO k = 1,3
DO l = 1,3
R_t_ip(k,l) = sum(Basis(1:n)*R_t(k,l,1:n))
DO k = 1, MIN(3, SIZE(Reluct_Z,1))
DO l = 1, MIN(3, SIZE(Reluct_Z,2))
Nu(k,l) = sum(Basis(1:n)*Reluct_Z(k,l,1:n))
END DO
END DO
END IF
w_dens = 0.5*SUM(B(1,:)*MATMUL(R_t_ip,B(1,:)))
R_ip = 0.0d0
ELSE
R_ip_Z = SUM(Basis(1:n)*R_Z(1:n))
DO k=1,3
Nu(k,k) = R_ip_Z
END DO
!
! The calculation of the Maxwell stress tensor doesn't yet support
! a tensor-form reluctivity. Create the scalar reluctivity parameter
! so that the Maxwell stress tensor may be calculated. The complex
! part will be ignored.
!
R_ip = REAL(R_ip_Z)
END IF
IF (RealField) THEN
w_dens = 0.5*SUM(B(1,:)*MATMUL(REAL(Nu), B(1,:)))
ELSE
! This yields twice the time average:
w_dens = 0.5*( SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:)) + &
SUM(MATMUL(REAL(Nu), B(2,:)) * B(2,:)) )
END IF
w_dens = 0.5*R_ip*SUM(B(1,:)**2)
END IF
PR_ip = SUM( Basis(1:n)*PR(1:n) )

Expand All @@ -1385,34 +1419,32 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)

IF (ASSOCIATED(NF).OR.ASSOCIATED(EL_NF)) THEN
NF_ip = 0._dp
B2 = sum(B(1,:)*B(1,:) + B(2,:)*B(2,:))
IF (RealField) THEN
HdotB = SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:))
ELSE
HdotB = SUM(MATMUL(REAL(Nu), B(1,:)) * B(1,:)) + &
SUM(MATMUL(REAL(Nu), B(2,:)) * B(2,:))
END IF
DO k=1,n
DO l=1,3
DO m=1,3
NF_ip(k,l) = NF_ip(k,l) - (R_ip*(B(1,l)*B(1,m)))*dBasisdx(k,m)
END DO
NF_ip(k,l) = NF_ip(k,l) + (R_ip*B2-w_dens)*dBasisdx(k,l)
val = SUM(dBasisdx(k,1:3)*B(1,1:3))
NF_ip(k,l) = NF_ip(k,l) - SUM(REAL(Nu(l,1:3)) * B(1,1:3)) * val + &
(HdotB-w_dens)*dBasisdx(k,l)
END DO
END DO

IF (.NOT. RealField) THEN
DO k=1,n
DO l=1,3
DO m=1,3
NF_ip(k,l) = NF_ip(k,l) - (R_ip*(B(2,l)*B(2,m)))*dBasisdx(k,m)
END DO
val = SUM(dBasisdx(k,1:3)*B(2,1:3))
NF_ip(k,l) = NF_ip(k,l) - SUM(REAL(Nu(l,1:3)) * B(2,1:3)) * val
END DO
END DO
END IF
END IF

Energy(1) = Energy(1) + s*0.5*PR_ip*SUM(E**2)

IF(ASSOCIATED(HB) .AND. RealField) THEN
Energy(2) = Energy(2) + s*w_dens
ELSE
Energy(2) = Energy(2) + s*0.5*R_ip*SUM(B**2)
END IF
Energy(2) = Energy(2) + s*w_dens

DO p=1,n
DO q=1,n
Expand All @@ -1426,13 +1458,20 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)

IF ( (ASSOCIATED(MFS).OR.ASSOCIATED(EL_MFS)) .and. .not. HasZirka) THEN
IF(.NOT. HasZirka) then
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)+s*(R_ip*B(1,:)-REAL(MG_ip))*Basis(p)
k = k+3
IF ( Vdofs>1 ) THEN
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)+s*(R_ip*B(2,:)-AIMAG(MG_ip))*Basis(p)
IF (RealField) THEN
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + &
s * (MATMUL(REAL(Nu), B(1,:)) - REAL(MG_ip)) * Basis(p)
k = k+3
ELSE
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + s * &
(MATMUL(REAL(Nu), B(1,:)) - MATMUL(AIMAG(Nu), B(2,:)) - REAL(MG_ip)) * Basis(p)
k = k+3
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3) + s * &
(MATMUL(AIMAG(Nu), B(1,:)) + MATMUL(REAL(Nu), B(2,:)) - AIMAG(MG_ip)) * Basis(p)
k = k+3
END IF
ELSE
! Never here?
FORCE(p,k+1:k+3) = FORCE(p,k+1:k+3)-s*(REAL(MG_ip))*Basis(p)
END IF
END IF
Expand Down Expand Up @@ -1926,7 +1965,7 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
END IF

IF(ALLOCATED(Gforce)) DEALLOCATE(Gforce)
DEALLOCATE( MASS,FORCE,Tcoef,RotM, R_t )
DEALLOCATE( MASS,FORCE,Tcoef,RotM )

IF (LossEstimation) THEN
CALL ListAddConstReal( Model % Simulation,'res: harmonic loss linear',TotalLoss(1) )
Expand Down Expand Up @@ -2132,10 +2171,10 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
END SELECT
IF (.NOT. ActiveBoundaryElement(Element)) CYCLE

C = GetConstReal(BC, 'Layer Electric Conductivity', Found)
IF (ANY(ABS(C(1:n)) > AEPS)) THEN
R = GetConstReal(BC, 'Layer Relative Permeability', Found)
IF (.NOT. Found) R = 1.0_dp
C_ip = GetConstReal(BC, 'Layer Electric Conductivity', Found)
IF (ABS(C_ip) > AEPS) THEN
mu_r = GetConstReal(BC, 'Layer Relative Permeability', Found)
IF (.NOT. Found) mu_r = 1.0_dp
ELSE
CYCLE
END IF
Expand All @@ -2162,10 +2201,7 @@ SUBROUTINE MagnetoDynamicsCalcFields(Model,Solver,dt,Transient)
CALL GetEdgeBasis(Element, WBasis, RotWBasis, Basis, dBasisdx)
END IF

C_ip = SUM(Basis(1:n) * C(1:n))
R_ip = SUM(Basis(1:n) * R(1:n))
R_ip = 4.0d0 * PI * 1d-7 * R_ip
val = SQRT(2.0_dp/(C_ip * Omega * R_ip)) ! The layer thickness
val = SQRT(2.0_dp/(C_ip * Omega * 4.0d0 * PI * 1d-7 * mu_r)) ! The layer thickness
Zs = CMPLX(1.0_dp, 1.0_dp, KIND=dp) / (C_ip*val)

E(1,:) = Omega * MATMUL(SOL(2,np+1:nd), WBasis(1:nd-np,:)) - MATMUL(SOL(1,1:np), dBasisdx(1:np,:))
Expand Down

0 comments on commit bc9e064

Please sign in to comment.