How to define and update a nonlinear material on subdomains

dolfinx
1

Hi everyone,

I’m aiming to model an electric machine in which only the core material is considered nonlinear.
I have defined the reluctivity of the material as a function, which on some domains is a fixed value, and in other subdomains is defined by the material properties and the solution of the problem, which is hence nonlinear.
The reluctivity is defined by:

W = functionspace(mesh, ("DG", 1))
B_norm = Function(W)
B_expr = Expression(sqrt((Az.dx(1))**2 + (-Az.dx(0))**2),
                                    W.element.interpolation_points())
B_norm.interpolate(B_expr)

for key in material_tags.keys():
    out = self.geometry[key]
    cells = self.domain_data.cell_tags.find(out["geom"]["tag"])
    
    if hasattr(out["material"], "bh_curve"):
        nu_mat = B_norm.x.array[cells] / (out["material"].get_H(B_norm.x.array[cells]) + 1e-12)  # Avoid division by zero
        nu.x.array[cells] = nu_mat
    else:
         nu_mat = 1.0 / (out["material"].mu_0 * out["material"].mu_r)
          nu.x.array[cells] = np.full_like(cells, nu_mat, dtype=default_scalar_type)
nu.x.scatter_forward()

Az is the magnetic vector potential which is the unknown variable that is solved for. and out[“material”].get_H(B) interpolates the given B on a lookup table to get the value for H.

As the machine geometry is reduced, it is a periodic problem so I’m using dolfinx_mpc and the nonlinear solver from: dolfinx_mpc/python/tests/test_nonlinear_assembly.py at v0.9.3 · jorgensd/dolfinx_mpc · GitHub

My question is, how should I modify the newton solver to update the reluctivity in between solving iterations?

2

How about calling your ‘update_nu’ code at line 46 of the NonlinearMPCProblem class in your link?

3

Hi Stein, thanks for you reply.

Implementing your suggestions seems to work in that the reluctivity is updated at every iteration of the solver. I’m not entirely sure this update is actually used in the numerical solution as the solution explodes. I’m still investigating whether the rest of my problem is setup correctly.

Just out of curiosity, why would you call this update in the definition of the Jacobian: J, and not in the definition of the Residual F?

4

Ah, that would make more sense indeed.