MagnetoDynamics: JxB Lorentz force for complex fields

Numerical methods and mathematical models of Elmer
Post Reply
vencels
Posts: 30
Joined: 20 Sep 2016, 17:05
Antispam: Yes
Location: Latvia
Contact:

MagnetoDynamics: JxB Lorentz force for complex fields

Post by vencels » 12 Oct 2016, 14:32

Hi!

I am doing frequency-transient simulation with WhitneyAVHarmonicSolver (complex magnetic vector potential), then computing JxB with MagnetoDynamicsCalcFields and using resulting force in FlowSolver (Navier-Stokes), but it seems that Elmer does not support this functionality.

Firstly, it seems that Elmer is not computing current density for complex fields JatIP(2,l): MagnetoDynamics.F90 line 6477:

Code: Select all

           IF (Vdofs == 1) THEN
             ...
           ELSE
              DO l=1,3
                FORCE(p,k+l) = FORCE(p,k+l)+s*SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) )*Basis(p) - &
                  s * SUM( AIMAG(CMat_ip(l,1:3)) * E(2,1:3) ) * Basis(p) + s*REAL(BodyForceCurrDens_ip(l))*Basis(p)
              END DO
              k = k+3
              DO l=1,3
                FORCE(p,k+l) = FORCE(p,k+l)+s*SUM( AIMAG(CMat_ip(l,1:3)) * E(1,1:3) )*Basis(p) + &
                  s * SUM( REAL(CMat_ip(l,1:3)) * E(2,1:3) ) * Basis(p) + s*AIMAG(BodyForceCurrDens_ip(l))*Basis(p)
              END DO
              k = k+3
           END IF
         END IF
I guess it has to be something like this:

Code: Select all

           IF (Vdofs == 1) THEN
             ...
           ELSE
              DO l=1,3
                JatIp(1,l) = SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) ) - &
                             SUM( AIMAG(CMat_ip(l,1:3)) * E(2,1:3) ) + REAL(BodyForceCurrDens_ip(l))
                FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(1,l)*Basis(p)
              END DO
              k = k+3
              DO l=1,3
                JatIp(2,l) = SUM( AIMAG(CMat_ip(l,1:3)) * E(1,1:3) ) + &
                             SUM( REAL(CMat_ip(l,1:3)) * E(2,1:3) ) + AIMAG(BodyForceCurrDens_ip(l))
                FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(2,l)*Basis(p)
              END DO
              k = k+3
Next, I believe JxB is not computed for complex case, line 6500:

Code: Select all

         IF ( ASSOCIATED(JXB).OR.ASSOCIATED(EL_JXB)) THEN
           IF (.NOT. ASSOCIATED(CD) .AND. .NOT. ASSOCIATED(EL_CD)) THEN
             CALL Warn('MagnetoDynamicsCalcFields', 'Cannot Calculate JxB since Current Density is not calculated!')
           ELSE
             IF (Vdofs == 1) THEN
               JXBatIP(1,:) = crossproduct(JatIP(1,:),B(1,:))
               DO l=1,dim
                 FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p)
               END DO
               k = k+3
             END IF
           END IF
         END IF
I guess it has to be something like this:

Code: Select all

         IF ( ASSOCIATED(JXB).OR.ASSOCIATED(EL_JXB)) THEN
           IF (.NOT. ASSOCIATED(CD) .AND. .NOT. ASSOCIATED(EL_CD)) THEN
             CALL Warn('MagnetoDynamicsCalcFields', 'Cannot Calculate JxB since Current Density is not calculated!')
           ELSE
             IF (Vdofs == 1) THEN
               JXBatIP(1,:) = crossproduct(JatIP(1,:),B(1,:))
               DO l=1,dim
                 FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p)
               END DO
               k = k+3
             ELSE
               ! |X1|   |Y1|   |Conj(X2*Y3-X3*Y2)|
               ! |X2| x |Y2| = |Conj(X3*Y1-X1*Y3)|
               ! |X3|   |Y3|   |Conj(X1*Y2-X2*Y1)|

               JXBatIP(1,1) = JatIP(1,2)*B(1,3) - JatIP(1,3)*B(1,2) - JatIP(2,2)*B(2,3) + JatIP(2,3)*B(2,2)
               JXBatIP(1,2) = JatIP(1,3)*B(1,1) - JatIP(1,1)*B(1,3) - JatIP(2,3)*B(2,1) + JatIP(2,1)*B(2,3)
               JXBatIP(1,3) = JatIP(1,1)*B(1,2) - JatIP(1,2)*B(1,1) - JatIP(2,1)*B(2,2) + JatIP(2,2)*B(2,1)
               JXBatIP(2,1) = - JatIP(1,2)*B(2,3) - JatIP(2,2)*B(1,3) + JatIP(1,3)*B(2,2) + JatIP(2,3)*B(1,2)
               JXBatIP(2,2) = - JatIP(1,3)*B(2,1) - JatIP(2,3)*B(1,1) + JatIP(1,1)*B(2,3) + JatIP(2,1)*B(1,3)
               JXBatIP(2,3) = - JatIP(1,1)*B(2,2) - JatIP(2,1)*B(1,2) + JatIP(1,2)*B(2,1) + JatIP(2,2)*B(1,1)
Finally, I modified function LorentzForce in Differentials.F90:

Code: Select all

    TYPE(Variable_t), POINTER :: JxB
    REAL(KIND=dp), POINTER :: JxBvalues(:)
    INTEGER, POINTER :: JxBperm(:)
    INTEGER :: JxBdofs

    JXB => NULL()
    JxB => VariableGet( CurrentModel % Variables, 'JxB' )
    IF ( .NOT.ASSOCIATED( JxB ) ) THEN
      CALL Fatal( 'Differentials', 'Cannot access JxB from MagnetoDynamics.' )
    ELSE
      JxBvalues => JxB % Values
      JxBperm => JxB % Perm
      JxBdofs = JxB % Dofs
      stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, Basis,dBasisdx )
      L(1) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+1 ) )
      L(2) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+2 ) )
      L(3) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+3 ) )
      RETURN
    END IF
My simulation consists of conducting cylinder and one-turn coil around it, then I am comparing simulation results with Comsol. Magnetic field looks right, but I get wrong JxB. Could someone check whether my modifications to the code are correct?

Post Reply