How to output user defined variable to the vtu file?

Post processing utility for Elmer
raback
Site Admin
Posts: 3240
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

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

Post by 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

Jenwel
Posts: 19
Joined: 22 Jun 2017, 07:07
Antispam: Yes

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

Post by 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: 19
Joined: 22 Jun 2017, 07:07
Antispam: Yes

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

Post by 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

raback
Site Admin
Posts: 3240
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

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

Post by 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

kimmysawi
Posts: 1
Joined: 09 Jul 2018, 14:03
Antispam: Yes

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

Post by kimmysawi » 13 Jul 2018, 17:16

Jenwel wrote:
31 Aug 2017, 06:03
Hi Peter,
Thank you again for your replay.
For the 2D magnetic field, actually I only need the norm value (scalar value) of magnetic flux density of each nodal.
Now, I tried to find corresponding variable Bx and By in MagnetoDynamics2D.F90 or CalcFields.F90, but I am not quite sure about that I have found them.

Can you give me any hints? see reviews of phenq diet pills Maybe I can directly calculate the core loss after the flux density field has been solved.

Now the solvers I used in sif file includes:

Code: Select all

Solver 2
Exec Solver = Always
Equation = MgDyn2D
Variable = A
Procedure = "MagnetoDynamics2D" "MagnetoDynamics2D"

Stabilize = False
  
Partition Local Constraints = Logical True
  
Nonlinear System Compute Change in Scaled System = Logical True
  
Nonlinear System Convergence Measure = Residual

Nonlinear System Convergence Tolerance = 1e-6
  
Nonlinear System Max Iterations = 20
  
Nonlinear System Min Iterations = 1
  
Nonlinear System Relaxation Factor = 1.0
  
Nonlinear System Newton After Iterations = 7

  
Export Lagrange Multiplier = Logical True

  
Linear System Abort Not Converged = False
  
Linear System Solver = Iterative
  
Linear System Iterative Method = GCR 
  
Linear System GCR Restart = 500
  
Bicgstabl Polynomial Degree = 4

  
Linear System Preconditioning = ILU4
  
Linear System Max Iterations = 1500
  
Linear System Residual Output = Integer 20
  
Linear System Convergence Tolerance = 1e-7 ! 2.0e-6

  
Mortar BCs Additive = Logical True

End

Solver 3
Exec Solver = Always
Equation = CalcFields
Potential Variable = "A"
Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
  
Calculate Nodal Forces = Logical True
  
Calculate Magnetic Torque = Logical True

  
Calculate Magnetic Force = Logical True

  
Calculate Current Density = Logical True 
  
Calculate Magnetic Vector Potential = Logical True

Frequency = Real $ f

Linear System Solver = Direct
Linear System Iterative Method = MUMPS
End

I am not sure the solver "Bsolver" in MagnetoDynamics2D.F90 is called by the solver "MagnetoDynamicsCalcFields" or other solver?
If not, why the magnetic flux density still can be calculated and saved in the vtu file?

Best Regards
Jenwel
All was meant to be those variable and defined directly to the tape

Post Reply