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
How to output user defined variable to the vtu file?
Re: How to output user defined variable to the vtu file?
Hi, Peter
Thank you so much for the answers! I understand now!
I revised the code as:
Now only the transient loss of elements belong to electrical steel are solved, as shown in the following figures.
Fig 01
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
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
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
Re: How to output user defined variable to the vtu file?
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 2D element?
Thank you in advance!
BR
Jenwel
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 2D element?
Thank you in advance!
BR
Jenwel

 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?
Hi
The total Joule heating calculate below illustrates how to integrate over the mesh:
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
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
Peter
Re: How to output user defined variable to the vtu file?
All was meant to be those variable and defined directly to the tapeJenwel wrote: ↑31 Aug 2017, 06:03Hi 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 = 1e6 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 = 1e7 ! 2.0e6 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