WhitneyAVHarmonicSolver - "unphysical" solution

[gbassi]

Hi all,
I am currently experimenting with Elmer to simulate PCB traces in different frequency conditions to get quantities such as S-parameters.
As a first step I took the "Magnetic field induced by harmonic current in a wire" example from Elmer tutorials and changed the geometry to be somehow closer to an actual trace with a rectangular section (1mm x 1mm x 10cm).

strip_geometry.png (25.32 KiB) Viewed 348 times

From geometry to simulation I performed the following steps (please find the relevant files at this link):
- designed the geometry in FreeCAD and exported .step file
- imported .step file in Salome, meshed the geometry and exported the .unv file
- opened the .unv file in ElmerGUI and created the mesh.* files
- run the simulation using the case.sif file
- checked the output .vtu file with Paraview

With a Linear System Convergence Tolerance = 1.0e-10, the simulation converges in about 150 steps. However from Paraview I can see that the current density is concentrated close to the input port and the current does not flow through the entire trace length.

Can someone please help me in finding the source of this issue?
And maybe a more general question: are the tools that I am currently using (WhitneyAVHarmonicSolver, MagnetoDynamicsCalcFields and SaveScalars with boundary int mean) the correct ones to simulate the frequency response of PCB traces?

Any comment/suggestion is very welcome as I started learning Elmer recently.

Thank you

Giovanni

[Tapegoji]

Hi,

I am doing a similar thing with the harmonic solver. The first thing I noticed is that you don't have Permeability and Permittivity defined in your materials.
Here are the material data I used. Note that it's in yaml format. Also for the boundary setup, here are my setup and I think they are working correctly.

Try these and see if that makes any difference. I did not quite see the problem. maybe you can send a picture of what you think is not right?

Furthermore, I noticed that you have some other boundaries that you have not defined set how Elmer needs to deal with them. if you look at the magnetic harmonic wire example, there are only three boundaries and they are defined. What I do is to make the common faces between copepr and the surrounding air to be the same. in Gmsh it's called fragmentation. It's even the same in Ansys. you need to define what should happen from copper to air at those boundaries (surfaces). when you make those common surfaces to be the same you are essentially telling the solver thay these two surfaces have tight coupling.

Try these suggestions and let me know how it goes.

Code: Select all

air:
  Density: 1.1885  # 20°C
  Electric Conductivity: 0.0 
  Heat Capacity: 1006.4  # 20°C
  Heat Conductivity: 0.025873  # 20°C
  Relative Permeability: 1
  Relative Permittivity: 1

copper:
  Density: 8960.0  # 0°C
  Electric Conductivity: 32300000  # 200°C
  Emissivity: 0.012  # 327°C
  Heat Capacity: 415.0  # 200°C
  Heat Conductivity: 401.0  # 0°C
  Relative Permeability: 1
  Relative Permittivity: 1
  
! pos
Boundary Condition 1
  Target Boundaries(1) = 3
  Name = "VoltageP"
  AV im {e} = 0
  AV re {e} = 0
  AV im = 0
  AV re = 0.01  
End

! neg
Boundary Condition 2
  Target Boundaries(1) = 4
  Name = "VoltageN"
  AV im {e} = 0
  AV re {e} = 0
  AV im = 0
  AV re = 0.0
End

! air_boundary
Boundary Condition 3
  Target Boundaries(1) = 5
  Name = "AxialField"
  AV re {e} = 0
  AV im {e} = 0
End
  

[gbassi]

Hi,
Thanks a lot for for your reply and detailed instructions!
Here is the new .sif file where I applied the suggestions (both copper and air now have Relative Permeability and Relative Permittivity set and boundary conditions are now the same as in the provided example).

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 6
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals(1) = 1
  Coordinate Scaling = 1.0e-3
  Angular Frequency = $ 2*pi*0.001e9
  Solver Input File = case.sif
  Post File = case.vtu
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.670374419e-08
  Permittivity of Vacuum = 8.85418781e-12
  Permeability of Vacuum = 1.25663706e-6
  Boltzmann Constant = 1.380649e-23
  Unit Charge = 1.6021766e-19
End

Body 1
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1
End

Body 2
  Target Bodies(1) = 2
  Name = "Body 2"
  Equation = 1
  Material = 2
End

Solver 1
  Equation = MgHarm
  Procedure = "MagnetoDynamics" "WhitneyAVHarmonicSolver"
  Electrodynamics Model = Logical True
  Exec Solver = Always
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStabl
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-7
  BiCGstabl polynomial degree = 4
  Linear System Preconditioning = none
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Solver 2
  Equation = MgDynPost
  Calculate Joule Heating = True
  Calculate Magnetic Field Strength = True
  Calculate Electric Field = Logical True
  Calculate Electric Field Strength = Logical True
  Calculate Current Density = True
  Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
  Discontinuous Bodies = True
  Exec Solver = Before Saving
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-7
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = ILU0
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Solver 3
  Equation = Result Output
  Output Format = Vtu
  Output File Name = case_0.001_GHz
  Binary Output = False
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Exec Solver = after All
End

Solver 4
  Filename = "data_0.001_GHz.txt"
  Exec Solver = after All
  Equation = SaveScalars
  Procedure = "SaveData" "SaveScalars"

  Variable 1 = electric field re 1
  Operator 1 = String "boundary int mean"
  mask name 1 = inport
  Variable 2 = electric field re 2
  Operator 2 = String "boundary int mean"
  mask name 2 = inport
  Variable 3 = electric field re 3
  Operator 3 = String "boundary int mean"
  mask name 3 = inport
  Variable 4 = electric field im 1
  Operator 4 = String "boundary int mean"
  mask name 4 = inport
  Variable 5 = electric field im 2
  Operator 5 = String "boundary int mean"
  mask name 5 = inport
  Variable 6 = electric field im 3
  Operator 6 = String "boundary int mean"
  mask name 6 = inport

  Variable 7 = electric field re 1
  Operator 7 = String "boundary int mean"
  mask name 7 = outport
  Variable 8 = electric field re 2
  Operator 8 = String "boundary int mean"
  mask name 8 = outport
  Variable 9 = electric field re 3
  Operator 9 = String "boundary int mean"
  mask name 9 = outport
  Variable 10 = electric field im 1
  Operator 10 = String "boundary int mean"
  mask name 10 = outport
  Variable 11 = electric field im 2
  Operator 11 = String "boundary int mean"
  mask name 11 = outport
  Variable 12 = electric field im 3
  Operator 12 = String "boundary int mean"
  mask name 12 = outport
End

Equation 1
  Name = "Equation 1"
  Active Solvers(4) = 1 2 3 4
End

Material 1
  Name = "Copper (generic)"
  Heat Conductivity = 401.0
  Electric Conductivity = 59.59e6
  Poisson ratio = 0.34
  Relative Permeability = 0.999994
  Youngs modulus = 115.0e9
  Heat expansion Coefficient = 16.5e-6
  Density = 8960.0
  Sound speed = 3810.0
  Relative Permittivity = 1.0
  Heat Capacity = 385.0
End

Material 2
  Name = "Air (room temperature)"
  Relative Permeability = 1.00000037
  Viscosity = 1.983e-5
  Heat expansion Coefficient = 3.43e-3
  Heat Conductivity = 0.0257
  Density = 1.205
  Sound speed = 343.0
  Relative Permittivity = 1.00059
  Heat Capacity = 1005.0
End

Boundary Condition 1
  Target Boundaries(1) = 3
  Name = "Ground"
  AV re {e} = 0.0
  AV im {e} = 0.0
  AV re = 0.0
  AV im = 0.0
  outport = Logical True
End

Boundary Condition 2
  Target Boundaries(1) = 1
  Name = "Voltage"
  AV re {e} = 0.0
  AV im {e} = 0.0
  AV re = 0.01
  AV im = 0.0
  inport = Logical True
End

Boundary Condition 3
  Target Boundaries(3) = 2 4 5
  Name = "AxialField"
  AV re {e} = 0
  AV re {e} = 0
End

I also modified the geometry to have common contact surfaces between copper and air. Please find the updated files here .
Despite the changes I still observe the "unphysical" behaviour as no current is flowing out of the output port whereas at the input port I observe current flowing in, as shown in the attached screenshots from Paraview.

outport_current.png

(20.39 KiB) Not downloaded yet

inport_current.png

(30.75 KiB) Not downloaded yet

Please note that with the current mesh and .sif configuration the simulation does not converge with a Linear System Convergence Tolerance = 1.0e-10 so I had to increase this value to 1.0e-7.

I have the feeling that I am doing something wrong in the mesh side but I haven't figured out what exactly yet.

Would you have any suggestion on how to solve this issue?

Thank you

Giovanni

[kevinarden]

Showing a mid-plane cut with a plot along the length.

plot.png

(83.56 KiB) Not downloaded yet

[gbassi]

Hi,
thank you for the plot!
I think this shows more clearly that there is something wrong. The equivalent plot for the "Magnetic field induced by harmonic current in a wire" example shows a constant current along the length of the wire, as expected.

Do you have an idea of what might be the cause of this issue?

Thank you

Giovanni

[gbassi]

After some research I found that there is the possibility to define surface impedance boundary conditions that, according to this article, should be applicable to my case (high frequency). So I modified the .sif file as follows

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 6
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals(1) = 1
  Coordinate Scaling = 1.0e-3
  Angular Frequency = $ 2*pi*0.001e9
  Solver Input File = case_v2.sif
  Post File = case_v2.vtu
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.670374419e-08
  Permittivity of Vacuum = 8.85418781e-12
  Permeability of Vacuum = 1.25663706e-6
  Boltzmann Constant = 1.380649e-23
  Unit Charge = 1.6021766e-19
End

Body 1
  Target Bodies(1) = 1
  Name = "Wire"
  Equation = 1
  Material = 1
End

Body 2
  Target Bodies(1) = 2
  Name = "Air"
  Equation = 1
  Material = 2
End

Solver 1
  Equation = MgHarm
  Procedure = "MagnetoDynamics" "WhitneyAVHarmonicSolver"
  Use Piola Transform = True
  Electrodynamics Model = Logical True
  Exec Solver = Always
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStabl
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-7
  BiCGstabl polynomial degree = 4
  Linear System Preconditioning = none
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
End

Solver 2
  Equation = MgDynPost
  Calculate Joule Heating = True
  Calculate Magnetic Field Strength = True
  Calculate Electric Field = Logical True
  Calculate Electric Field Strength = Logical True
  Calculate Current Density = True
  Calculate Magnetic Vector Potential = True
  Calculate Nodal Heating = True
  Calculate Nodal Current = Logical True
  Procedure = "MagnetoDynamics" "MagnetoDynamicsCalcFields"
  Discontinuous Bodies = True
  Exec Solver = Before Saving
  Stabilize = True
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 20
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations = 500
  Linear System Convergence Tolerance = 1.0e-7
  BiCGstabl polynomial degree = 2
  Linear System Preconditioning = ILU0
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = False
  Linear System Residual Output = 10
  Linear System Precondition Recompute = 1
  
  Calculate Nodal Fields = Logical False
  Impose Body Force Potential = Logical True
  Impose Body Force Current = Logical True

  Discontinuous Bodies = True
  calculate harmonic peak power = logical False
End

Solver 3
  Equation = Result Output
  Output Format = Vtu
  Output File Name = case_0.001_GHz
  Binary Output = False
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Exec Solver = after All
  Discontinuous Bodies = Logical True
  Save Geometry IDs = Logical True
End

Solver 4
  Filename = "data_0.001_GHz.txt"
  Exec Solver = after All
  Equation = SaveScalars
  Procedure = "SaveData" "SaveScalars"

  Variable 1 = electric field re 1
  Operator 1 = String "boundary int mean"
  mask name 1 = inport
  Variable 2 = electric field re 2
  Operator 2 = String "boundary int mean"
  mask name 2 = inport
  Variable 3 = electric field re 3
  Operator 3 = String "boundary int mean"
  mask name 3 = inport
  Variable 4 = electric field im 1
  Operator 4 = String "boundary int mean"
  mask name 4 = inport
  Variable 5 = electric field im 2
  Operator 5 = String "boundary int mean"
  mask name 5 = inport
  Variable 6 = electric field im 3
  Operator 6 = String "boundary int mean"
  mask name 6 = inport

  Variable 7 = electric field re 1
  Operator 7 = String "boundary int mean"
  mask name 7 = outport
  Variable 8 = electric field re 2
  Operator 8 = String "boundary int mean"
  mask name 8 = outport
  Variable 9 = electric field re 3
  Operator 9 = String "boundary int mean"
  mask name 9 = outport
  Variable 10 = electric field im 1
  Operator 10 = String "boundary int mean"
  mask name 10 = outport
  Variable 11 = electric field im 2
  Operator 11 = String "boundary int mean"
  mask name 11 = outport
  Variable 12 = electric field im 3
  Operator 12 = String "boundary int mean"
  mask name 12 = outport
End

Equation 1
  Name = "Equation 1"
  Active Solvers(4) = 1 2 3 4
End

Material 1
  Name = "Copper (generic)"
  Heat Conductivity = 401.0
  Electric Conductivity = 0.0 !59.59e6
  Poisson ratio = 0.34
  Relative Permeability = 1 !0.999994
  Youngs modulus = 115.0e9
  Heat expansion Coefficient = 16.5e-6
  Density = 8960.0
  Sound speed = 3810.0
  !Relative Permittivity = 1.0
  Permittivity = 0.0
  Heat Capacity = 385.0
End

Material 2
  Name = "Air (room temperature)"
  Relative Permeability = 1.0 !1.00000037
  Viscosity = 1.983e-5
  Heat expansion Coefficient = 3.43e-3
  Heat Conductivity = 0.0257
  Density = 1.205
  Sound speed = 343.0
  !Relative Permittivity = 1.00059
  Permittivity = 0.0
  Heat Capacity = 1005.0
  Electric Conductivity = 0.0
End

Boundary Condition 1
  Target Boundaries(1) = 3
  Name = "Ground"
  !AV im {e} 2 = 0
  !AV re {e} 2 = 0
  !AV im {e} 3 = 0
  !AV re {e} 3 = 0
  AV re = 0.0
  AV im = 0.0
  outport = Logical True
End

Boundary Condition 2
  Target Boundaries(1) = 1
  Name = "Voltage"
  !AV im {e} 2 = 0
  !AV re {e} 2 = 0
  !AV im {e} 3 = 0
  !AV re {e} 3 = 0
  AV re = 0.01
  AV im = 0.0
  inport = Logical True
End

Boundary Condition 3
  Target Boundaries(1) = 6
  Name = "WireSurface"
  Layer Electric Conductivity = Real 58e6
  Layer Relative Permeability = Real 1
End

Boundary Condition 4
  Target Boundaries(3) = 2 4 5
  Name = "AxialField"
  AV im {e} 1 = 0
  AV re {e} 1 = 0
  AV im {e} 2 = 0
  AV re {e} 2 = 0
  AV im {e} 3 = 0
  AV re {e} 3 = 0
End

and got some more realistic results, where at least the current reaches the output side.

surface_conduction.png

(91.09 KiB) Not downloaded yet

The mesh and project files can be found at the following link.

On the still-to-investigate side:
- even with this "simple" geometry around 1000 steps are needed to converge at 1.0e-7 level.
- starting with the available quantities I still don't know how to compute s-parameters.

As always suggestions are very welcome

Giovanni

[Tapegoji]

Glad you figured it out. I was trying to recreate your case by using gmsh to mesh it. However, I had convergance issues.

One thing I noticed was that, your frequency was really high, 1 MHz. You had fine mesh on the surfaces but was it fine enough for Elmer to give you good results? I don't know.

I want to do some similar work for PCB impedances as well. Let me know if you like to chat and see how we collaborate.

[kevinarden]

Mesh Quality looks like an improvement on the original project_v2

hexamesh.png

(67.6 KiB) Not downloaded yet

ends.png

(483.65 KiB) Not downloaded yet

[gbassi]

Hi,
Definitely much better!
Would you be so kind to share the instructions on how to create this mesh?

Thank you

Giovanni

[Tapegoji]

In Gmsh this can be generated using transfinite method. it would be easy in this simple geometry you have.
However, since you want to do pcb simulation, transfinite would be hard.
You mean want to look into sweep function for doing pcb simulation. this might give you a good quality mesh. I suggest making a 2D quadratic surface mesh and then sweep it for the thickness of the pcb. you can specify how many layers you want for in that depth.

Regards,