Good morning community;

is there a way to implment bh curve using ngsolve :

i do somthing like

def nuf(u):

return Norm(grad(u))

nu=CoefficientFunction(nuf)

a+=nu*grad(u)*grad(v)*dx

but always i got nan values !

any tricks ?

Dear Lahoussaine,

what is the actual model you want to implement?

âNorm(grad(u))â does not seem to a suitable model in first place, as it leads to a singular matrix for grad(u) = 0.

Best,

Matthias

I want something like

1/(B**3+B**2+B )

or :

3.8*exp(2.17*B*B)+396.2

I ve tried those expression but it gives always NAN;

Thankes for your reply

TBH, I have never seen any magnetic model like that (which might not mean a lotâŠ)

I canât see how the first one could give something reasonable for B=0.

Before using a nonlinear model, did you have a working version for a constant ânuâ?

I fear, we need much more information to be able to help you.

Best

def nuf(u):

IfPos(2.28-B,(1/(-0.003995*B**3+0.005956*B**2+0.00646*B + 0.001925)), ((1902.211*B-4266.42343)/0.016694669))

Here is my formule, but it doesnât work

The problem is whenever i put norm(B) it returns NAN

even if the formula is well defined( 1/(B**2+1) for example )

Thank you and Sorry to not be claire

Unfortunately, I still miss some information: how is âBâ defined?

Where does the âNormâ enter in "1/(B**2 + 1)

Which command did you use for assembly? Where exactly do you get NaN? Already in the matrix or only after solving, i.e. in the solution?

Could you maybe post a complete minimal example that reproduces your problem?

Here is the code :

u, v = fes.TnT()

mu_r = {âairâ: 1, âmagnetâ: 1.1, âfer1â: 1000, âfer2â: 1000}

# Definition des permeabilites -------------------------

def nuf(u):

return IfPos(B-2.28,(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)), ((1902.211579*BÂ±4266.42343)/0.016694669))

B = Norm(grad(u))

nu_test = [nuf(B) if mat==âfer1â else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]

nu = CoefficientFunction(nu_test)

a = BilinearForm(fes, symmetric=False)

a += nu *(grad(u)*grad(v))*dx

âŠ

You still use Norm(grad(u)) in the bilinear form instead of nu.

Is this an oversight?

Did you call a.Assemble or a.AssembleLinearization further doen behind the dots?

Here is the reste of code

a += nu*(grad(u)*grad(v))*dx

# f= LinearForm(fes)

j = CoefficientFunction(10*1e3)

a+=- (j*v) * dx(âmagnetâ)

def SimpleNewtonSolve(gfu, a, tol=1e-8, maxits=30): res = gfu.vec.CreateVector()

du = gfu.vec.CreateVector()

fes = gfu.space

for it in range(maxits):

print(âIteration {:3} â.format(it), end=ââ)

a.Apply(gfu.vec, res)

a.AssembleLinearization(gfu.vec)

du.data = a.mat.Inverse(fes.FreeDofs()) * res

gfu.vec.data -= du

```
# stopping criteria
stopcritval = sqrt(abs(InnerProduct(du, res)))
print("<A u", it, ", A u", it, ">_{-1}^0.5 = ", stopcritval)
if stopcritval < tol:
break
```

# Solve the linear system

gfa = GridFunction(fes)

SimpleNewtonSolve(gfa, a)

Okay, there a couple of issues:

First of all, using grad(u) and grad(v) suggests that you are using the (auxiliary) scalar magnetic potential (H=-grad(u)). However, everywhere else you work like youâd use the vector potential (B=curl(A)), i.e. using nu=1/mu and j*v as loading term. So, to me, the whole formulation seems to be inconsistent. Could you recheck, please?

Second, when using a.AssembleLinearization(âŠ), ngsolve will linearize the expressions it finds in the bilinear form in some way. However, the nu*grad(u)*grad(v) is a linearized form already.

When trying to linearize this term, you get NaNâs because of a division by zero (or a very small number) resulting from the norm (at B=0).

When you work with CFs representing already âlinearizedâ forms,

they must not depend on trial functions but grid functions only.

So nu must not depend on âuâ but âgfuâ in this case.

Please note that I did not try to fix the example, so there may be some things that I have overlooked (also in my own explanations). I hope they are still useful though.

Best,

Matthias

Thank you for your time, I will just try to clarify my problem one last time.

You will find below the code of the Linear solution.

Concerning the linear solution, the weak formulation has been simplified, so the equation might be not clear, sorry for that. However, this simulation has been validated both with a commercial solver and another FEA api.

My problem is not about the weak formulation but the way to carry on a NonLinear simulation, the NonLinearity coming from the nu expression which depend on B and thus on u.

The Nonlinear simulation is thus the same as the linear one but with nu(B)

The question is : How would you perform a Non linear algorithm iterating on nu ?

Youâll find below the linear code that worked for me âŠ

stator = Circle(center=(0, 0), radius=10, bc=âbc_rect_extâ )

rotor = Circle(center=(0, 0), radius=6)

circle_airgap = Circle(center=(0, 0), radius=7.5)

magnet1 = Circle(center=(0, 0), radius=7)

magnet_f = magnet1-rotor

airgap = circle_airgap-(magnet1)

stator_f = stator-circle_airgap

# change domain name and maxh

magnet_f.Mat(âmagnetâ).Maxh(0.5)

rotor.Mat(âfer1â).Maxh(0.5)

stator_f.Mat(âfer2â).Maxh(0.5)

airgap.Mat(âairâ).Maxh(0.5)

geo.Add(magnet_f)

geo.Add(rotor)

geo.Add(stator_f)

geo.Add(airgap)

mesh = Mesh(geo.GenerateMesh(maxh=0.1))

Draw(mesh)

# FEM -------------------------------------------------

fes = H1(mesh, order=3, dirichlet=âbc_rect_extâ)

u, v = fes.TnT()

mu_r = {âairâ: 1, âmagnetâ: 1.1, âfer1â: 1000, âfer2â: 1000}

# Definition des permeabilites -------------------------

mu0 = 4*pi*1e-7

nu_coef = [1/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]
nu = CoefficientFunction(nu_coef)
a = BilinearForm(fes)
a += nu*(grad(u)

*grad(v))*

Mx = CoefficientFunction(Mcos(atan2(y,x))*1e0)

*dx*

M=80000

f = LinearForm(fes)

j = CoefficientFunction(101e3)M=80000

f = LinearForm(fes)

j = CoefficientFunction(10

Mx = CoefficientFunction(M

My = CoefficientFunction(0)

f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx(âmagnetâ)

with TaskManager():

a.Assemble()

f.Assemble()

# Solve the linear system

##
gfa = GridFunction(fes)

gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

In my previous post, I might have had a confusion of the meaning of nu. Is the your nu âH(B)â or âdH/dBâ evaluated for some B? I am used to work with âincrementalâ coefficients, ie âdH/dBâ. So *if you do not so*, then my second âissueâ in my previous response is not valid.

Assuming that your nu if such that H = nu B, then a.AssembleLinearization(gfa.vec) will linearized it for you.

However, linearizing Norm(B) when B is zero, results in a division by zero. You can work around that as follows (see comments in the code below, in particular around âdef nufâ):

fes = H1(mesh, order=3, dirichlet=âbc_rect_extâ)

u, v = fes.TnT()

mu_r = {âairâ: 1, âmagnetâ: 1.1, âfer1â: 1000, âfer2â: 1000}

# Definition des permeabilites

mu0 = 4*pi*1e-7

def nuf(u):

# You need to avoid division by zero in the linearization! Either by perturbing the norm of by providing

# the linearized norm for the norm being close to zero.

```
# Perturbation
#B = Norm(grad(u)+CF((1e-6, 1e-6)))
# manual linearization
B = Norm(grad(u))
return IfPos(B-2.28,
(1/(-0.0039954722*B**3+0.00595660*B**2+0.00646251*B + 0.001925625)),
IfPos(B - 1e-8,
(1902.211579*B+-4266.42343)/0.016694669,
1902.211579/0.016694669) # linearized nu (dh/dB) for small B
)
```

nu_coef = [nuf(u) if mat==âfer1â else 1.0/(mu0*mu_r[mat]) for mat in mesh.GetMaterials()]

nu = CoefficientFunction(nu_coef)

a = BilinearForm(fes)

a += nu*(grad(u)*grad(v)) dx
M=80000
f = LinearForm(fes)
j = CoefficientFunction(101e3)
Mx = CoefficientFunction(M*cos(atan2(y,x))*1e0)

My = CoefficientFunction(0)

f += (Mx*grad(v)[1]-My*grad(v)[0]) * dx(âmagnetâ)

### âSIMULATINGâ the first Newton step only

# define gfa upfront

gfa = GridFunction(fes)

with TaskManager():

# check that AssembleLinearization() does not give NaNs

```
a.AssembleLinearization(gfa.vec)
f.Assemble()
```

# Solve the linear system

gfa.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Hope that helps!

Dear community,

This topic is very important in the field of electromagnetic analysis.

However, it appears that the full script has not yet been explained in this topic.

I believe that NGSolve may solve the nonlinear problem of magnetic permeability.

Because there exist simulation results of transformers calculated by Prof. Schoeberl.Gallery - Wiki (tuwien.ac.at)

I am trying to adapt the following tutorial to nonlinear problems.

Symbolic definition of forms : magnetic field â NGS-Py 6.2.2307 documentation (ngsolve.org)

However, the task has not yet been completed.

This is because I donât know how to set the permeability for each element instead of the material.

I would be grateful if someone would give me on how to resolve this issue.

Best,

Daigo

Dear community,

I post code for the modified relaxation method whose convergence time is inferior to the Newtonâs method.

Modified relaxation method:

Flux Distribution in DC Machines On-Load and Overloads | IEEE Journals & Magazine | IEEE Xplore

```
def ModifiedRelaxationMethod(fes, u, v, f, funcNu, omega=0.3, tol=1e-3, maxits=20):
fesH1 = H1(fes.mesh, order=1, dirichlet="outer", nograds=True)
gfu = GridFunction(fes)
nu = makeCf(fes.mesh, makeGf(fesH1, funcNu(makeVec(fes.mesh, curl(gfu).Norm()))))
print (f"Initial value calculating...")
gfu, a = solve(u, v, f, nu)
for it in range(maxits):
print (f"Iteration {it+1:3} ",end="")
nu = makeCf(fes.mesh, makeGf(fesH1, funcNu(makeVec(fes.mesh, curl(gfu).Norm()))))
_gfu, _a = solve(u, v, f, nu)
du = omega * (_gfu.vec.FV().NumPy()[:] - gfu.vec.FV().NumPy()[:])
gfu.vec.FV().NumPy()[:] += du
#stopping criteria
stopcritval = np.sum(np.abs(du))/np.sum(np.abs(gfu.vec.FV().NumPy()[:]))
print (f'err:{stopcritval:5.3f}, {curl(gfu).Norm()(fes.mesh(0,0,0.5)):.2f}T')
if stopcritval < tol:
break
return gfu
```

All codes are shown in the colab below.

We know that Newtonâs method is shorter in convergence time than the modified relaxation method.

Newtonâs method:

Efficient Techniques for Finite Element Analysis of Electric Machines | IEEE Journals & Magazine | IEEE Xplore

I believe the following content may shed some light on this question.

3.7 Nonlinear problems â NGS-Py 6.2.1808 documentation (ngsolve.org)

I have a question.

How would I be able to get \delta A(u^i) in Maxwellâs Equations?

Best,

Daigo

I only tried with two dimensinal problem, but Newtom method implemented in NGSolve works with magnetic problems as well.

I assumed magnetic reclutivity \nu _r = 1 + \left| B \right|^2 .

```
import os, sys
from netgen.occ import *
from numpy import *
from ngsolve import *
d = 0.005
h = 0.002
xp = [-d/2, -d/2, 0, d/2, d/2, 0]
yp = [ h/2,-h/2, -h/2, -h/2, h/2, h/2]
face1 = MoveTo(xp[0], yp[0]).LineTo(xp[1], yp[1]).LineTo(xp[2], yp[2]).LineTo(xp[5], yp[5]).Close().Face()
face2 = MoveTo(xp[2], yp[2]).LineTo(xp[3], yp[3]).LineTo(xp[4], yp[4]).LineTo(xp[5], yp[5]).Close().Face()
face = Glue([face1, face2])
face.faces.name = "air"
face.edges.Nearest((-d/2, 0, 0)).name = "left"
face.edges.Nearest(( d/2, 0, 0)).name = "right"
face.edges.Nearest((-d/4, h/2, 0)).name = "top_1"
face.edges.Nearest((-d/4,-h/2, 0)).name = "bot_1"
face.edges.Nearest(( d/4, h/2, 0)).name = "top_2"
face.edges.Nearest(( d/4,-h/2, 0)).name = "bot_2"
face.edges.Nearest(( 0, 0, 0)).name = "center"
face.vertices.Nearest((0, -h/2, 0)).name = "gnd_1"
face.vertices.Nearest((0, h/2, 0)).name = "gnd_2"
face.edges[1].Identify(face.edges[3], "periodicy", IdentificationType.PERIODIC)
face.edges[4].Identify(face.edges[6], "periodicy", IdentificationType.PERIODIC)
from ngsolve import *
from ngsolve import Mesh
geo = OCCGeometry(face, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=h/10, quad_dominated=True) )
fes = H1(mesh, order=2, dirichlet='center')
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=True)
a += (1.0 + Norm(grad(u))**2)*(grad(u)*grad(v) - 1e6*v*x)*dx
c = Preconditioner(a, type="bddc")
gfAz = GridFunction(fes)
gfAz.Set(-1e6*(1/6*x**3 - 1/2*(d/2)**2*x))
res = gfAz.vec.CreateVector()
w = gfAz.vec.CreateVector()
for n in range(7):
a.Apply(gfAz.vec, res)
a.AssembleLinearization(gfAz.vec)
solvers.CG(sol=w.data, rhs=res.data, mat=a.mat, pre=c.mat, printrates='', tol=1e-16, maxsteps=10000)
# w.data = a.mat.Inverse(fes.FreeDofs()) * res
gfAz.vec.data -= w
```