How to output user defined variable to the vtu file?

Post processing utility for Elmer

Re: How to output user defined variable to the vtu file?

Postby raback » 08 Sep 2017, 18:56

Hi

Elemental fields are not constant here. Within each element there is a different value at each node. So this is as in "discontinuous galerkin" method. This way the code can in principle maintain some accuracy when using higher order edge elements. Of course for lowest order edge elements the \curl of the vector field is more or less constant within element so the only purpose of the elemental DG field is that it maintains the discontinuity of the data.

-Peter
raback
Site Admin
 
Posts: 3111
Joined: 22 Aug 2009, 11:57
Location: Espoo, Finland

Re: How to output user defined variable to the vtu file?

Postby Jenwel » 09 Sep 2017, 09:17

Hi, Peter

Thank you so much for the answers! I understand now!

I revised the code as:
Code: Select all
DO i = 1, GetNOFActive()
      Element => GetActiveElement(i)
      Bnode1 = sqrt(EL_MFD % Values(EPerm(Element % DGIndexes(1)) * 3 - 2)**2 + &
      EL_MFD % Values(EPerm(Element % DGIndexes(1)) * 3 - 1)**2)
      Bnode2 = sqrt(EL_MFD % Values(EPerm(Element % DGIndexes(2)) * 3 - 2)**2 + &
      EL_MFD % Values(EPerm(Element % DGIndexes(2)) * 3 - 1)**2)
      Bnode3 = sqrt(EL_MFD % Values(EPerm(Element % DGIndexes(3)) * 3 - 2)**2 + &
      EL_MFD % Values(EPerm(Element % DGIndexes(3)) * 3 - 1)**2)
      ! 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")
    IF (CurrentTimeStep>1) THEN
      IF ((Element % BodyId == 19) .OR. (Element % BodyId == 21)) THEN

         dB = Bnode1 - PrevB_E(EPerm(Element % DGIndexes(1)))
         PrevB_E(EPerm(Element % DGIndexes(1))) = Bnode1
         Phys_E % Values(EPerm(Element % DGIndexes(1))) = Hc * ABS(dB/dt) + Kh * Babs * ABS(dB/dt)
         Peddy_E % Values(EPerm(Element % DGIndexes(1))) = Keddy * (dB/dt)**2
         Pcoreloss_E % Values(EPerm(Element % DGIndexes(1))) = Phys_E % Values(EPerm(Element % DGIndexes(1))) +&
          Peddy_E % Values(EPerm(Element % DGIndexes(1)))

         dB = Bnode2 - PrevB_E(EPerm(Element % DGIndexes(2)))
         PrevB_E(EPerm(Element % DGIndexes(2))) = Bnode2
         Phys_E % Values(EPerm(Element % DGIndexes(2))) = Hc * ABS(dB/dt) + Kh * Babs * ABS(dB/dt)
         Peddy_E % Values(EPerm(Element % DGIndexes(2))) = Keddy * (dB/dt)**2
         Pcoreloss_E % Values(EPerm(Element % DGIndexes(2))) = Phys_E % Values(EPerm(Element % DGIndexes(2))) +&
          Peddy_E % Values(EPerm(Element % DGIndexes(2)))

         dB = Bnode3 - PrevB_E(EPerm(Element % DGIndexes(3)))
         PrevB_E(EPerm(Element % DGIndexes(3))) = Bnode3
         Phys_E % Values(EPerm(Element % DGIndexes(3))) = Hc * ABS(dB/dt) + Kh * Babs * ABS(dB/dt)
         Peddy_E % Values(EPerm(Element % DGIndexes(3))) = Keddy * (dB/dt)**2
         Pcoreloss_E % Values(EPerm(Element % DGIndexes(3))) = Phys_E % Values(EPerm(Element % DGIndexes(3))) +&
          Peddy_E % Values(EPerm(Element % DGIndexes(3)))
      END IF
    END IF
END DO


Now only the transient loss of elements belong to electrical steel are solved, as shown in the following figures.
CoreLoss_ElementField.jpg
CoreLoss_ElementField
(119.26 KiB) Not downloaded yet

Fig 01

CoreLoss_ElementField_AirGap.jpg
CoreLoss_ElementField_AirGap
(137.95 KiB) Not downloaded yet

Fig 02

Fig02 shows that the core loss in the air gap has been removed.

Now I would continue multiplying the core loss per volume with the area of each element.

Thank you again!

BR
Jenwel
Jenwel
 
Posts: 16
Joined: 22 Jun 2017, 07:07

Re: How to output user defined variable to the vtu file?

Postby Jenwel » 15 Sep 2017, 16:15

Hi, Peter

I want to ask that how can I obtain the area of each element in "CalcFields.F90"? Is there a variable that denotes the area of the element, in m^2 for the 2-D element?

Thank you in advance!

BR
Jenwel
Jenwel
 
Posts: 16
Joined: 22 Jun 2017, 07:07

Re: How to output user defined variable to the vtu file?

Postby raback » 15 Sep 2017, 17:10

Hi

The total Joule heating calculate below illustrates how to integrate over the mesh:
Code: Select all
Loop over elements
  Loop over integration points
       s = IP % s(j) * detJ
       ! The Joule heating power per unit volume: J.E = (sigma * E).E
       Coeff = SUM( MATMUL( REAL(CMat_ip(1:3,1:3)), TRANSPOSE(E(1:1,1:3)) ) * &
                   TRANSPOSE(E(1:1,1:3)) ) * Basis(p) * s
       Power = Power + Coeff

So the contribution in each Gaussian integration point is multiplied by the value of Basis function and integration weight times the element metric detJ. For linear element detJ is the area/volume of the element. This code was taken from CalcFields routine where you can see these lines in context.

-Peter
raback
Site Admin
 
Posts: 3111
Joined: 22 Aug 2009, 11:57
Location: Espoo, Finland

Previous

Return to ElmerPost

Who is online

Users browsing this forum: No registered users and 1 guest