Saturated BH curve

I am solving nonlinear Eddy Current Problem. The Iron reluctivity nu dependents on the flux density norm B0.

Without saturation:
nu =CoefficientFunction([3.8*exp(2.17*B0*B0)+396.2 if mat== "magnet" else 1.0/(1.0*mu0)
for mat in mesh.GetMaterials()])

The magnetization curve saturates, e.g., at 1.8 T. How to set this?

Thanks,
Xiaotao

I try this:
def nuf(B0):
if B0 >=1.8:
return 3.8*exp(2.17*1.8*1.+396.2
else:
return 3.8*exp(2.17*B0*B0)+396.2
nu =CoefficientFunction([nuf(B0) if mat== "inner" else 1.0/(1.0*mu0)
for mat in mesh.GetMaterials()])
# TypeError '>=' not supported between instances of 'ngsolve.fem.CoefficientFunction' and 'float'



Here is B0.

mesh is 2D
fes = HCurl(mesh, order=2, dirichlet="b|r|t|l", nograds=True)
gfu = GridFunction(fes)
gfu.Set(...) # gfu is an in-plane vector
B0=curl(gfu) # B0 has one component

Hi Xiaotao,

I think the "IfPos" function might help you. It checks if a CoefficientFunction is positive at the evaluation points. It is documented here .

By transforming "B0 >= 1.8" to "B0-1.8 >= 0", we can write the definition for the `nuf` function as:

Thank you, Bernd.
Now the function nuf works perfectly well.