Adding air gap to ferromagnetic material in A formulation

[Takala]

Hello,

I'm struggling a little bit with a very simple problem, please help me out:

In "The finite element method for electromagnetic modeling", G. Meunier, p. 249, there is a following description for adding an air gap for formulation A to magnetic material:

Let us suppose there is an integrating term "F" in "omega" that is constant along some path "e" then the volume integral

Code: Select all

int_omega(F)

could be replaced with a surface integral

Code: Select all

int_gamma(int_e(F))

where gamma is the surface of the thin region.

Now when int_e(F) is obviously eF (F constant along e). Thus we have transformed the volume elements into surface elements:

Code: Select all

int_gamma(eF)

In magnetostatics A formulation can be described with:
F=(1/mu * rot A, rot w)

Now when I try to write this in to the code:

Instead of:

Code: Select all

    CALL LocalMatrixBC(STIFF,FORCE,LOAD,Acoef,Element,n,nd )

I write this:

Code: Select all

     AirGapLength=GetConstReal( BC, 'Air Gap Length', Found)
     if (Found) then
       CALL LocalMatrixAirGapBC(STIFF,FORCE,LOAD,AirGapLength,Element,n,nd )
     else
       CALL LocalMatrixBC(STIFF,FORCE,LOAD,Acoef,Element,n,nd )
     end if

and in LocalMatrixAirGapBC, I write:

Code: Select all

!------------------------------------------------------------------------------
  SUBROUTINE LocalMatrixAirGapBC(  STIFF, FORCE, LOAD, AirGapL, Element, n, nd )
!------------------------------------------------------------------------------
    REAL(KIND=dp) :: LOAD(:,:), AirGapL(:)
    REAL(KIND=dp) :: STIFF(:,:), FORCE(:)
    INTEGER :: n, nd
    TYPE(Element_t), POINTER :: Element, Parent, Edge
!------------------------------------------------------------------------------
    REAL(KIND=dp) :: Basis(n),dBasisdx(n,3),DetJ,L(3),Normal(3)
    REAL(KIND=dp) :: WBasis(nd,3), RotWBasis(nd,3), A, F, TC, mu
    LOGICAL :: Stat
    INTEGER, POINTER :: EdgeMap(:,:)
    TYPE(GaussIntegrationPoints_t) :: IP
    INTEGER :: t, i, j, k, ii,jj, np, p, q

    TYPE(Nodes_t), SAVE :: Nodes
!------------------------------------------------------------------------------
    CALL GetElementNodes( Nodes )

    STIFF = 0.0_dp
    FORCE = 0.0_dp
    MASS  = 0.0_dp

    mu = 4 * pi * 1d-7

    ! Numerical integration:
    !-----------------------
    IP = GaussPoints(Element)

    np = n*MAXVAL(Solver % Def_Dofs(:,1))
    DO t=1,IP % n
       stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
                  IP % W(t), detJ, Basis, dBasisdx )

       CALL GetEdgeBasis(Element, WBasis, RotWBasis, Basis, dBasisdx)

       A  = SUM(Basis(1:n) * AirGapL(1:n))

       ! Compute element stiffness matrix:
       !---------------------------------------------------


       DO i = 1,nd-np
         DO j = 1,nd-np
           q = j+np
           STIFF(p,q) = STIFF(p,q) + A / mu * &
              SUM(RotWBasis(i,:)*RotWBasis(j,:))*detJ*IP%s(t)
         END DO
       END DO  
    END DO
!------------------------------------------------------------------------------
  END SUBROUTINE LocalMatrixAirGapBC
!------------------------------------------------------------------------------

In sif file I say

Boundary Condition 5
Target Boundaries(1) = $ airgap

Air Gap Length = Real 1000

End

Could you please tell me, why this does not have any effect on the result what so ever?

Best Regards,

Eelis

[Takala]

A brain fart: I forgot the force term. It was as simple as that.

Anyway, thanks for nothing

-E