Solver Problem in Nonlinear Magnetostatics


291
views
1
5 months ago by

Hi everyone,

I'm currently working on nonlinear magnetostatics. I'm using a convenient scalar potential formulation.
I want to use a constitutive model which takes into account the effect of the magnetic saturation, but I'm facing  problems trying to solve the nonlinear variational problem.
Here is a short working example to reproduce the error:
from dolfin import *

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["quadrature_degree"] = 4

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Dirichlet boundary::
bcl = DirichletBC(V, Constant(0), left)
bcr = DirichletBC(V, Constant(100e+3), right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            
v  = TestFunction(V)             
u  = Function(V)                 

# Kinematics
H = -grad(u)

# Invariant
IN = dot(H, H)


# parameters
mu0= 4.0*pi * 1e-13
phi=10
ms= 1e+6
a=mu0*ms*ms/phi
b=phi/ms

# Magnetic induction
B=mu0*H + a*b*tanh(b*(IN)**0.5)/((IN)**0.5)*H

# Variational Formulation
F = dot(B, grad(v))*dx

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['linear_solver'] = 'mumps'

solver.solve()

# Save solution in VTK format
file = File("potential.pvd");
file << u;
​

It seems like there is a linearization problem with the expression accounting for the saturation of the magnetization.

I will be grateful for any help.

Best regards,

Dustin
Community: FEniCS Project

2 Answers


4
5 months ago by
Emek  
Hi there,

please be cautious about dividing by the invariant, which is zero (initially). Use something like: IN=dot(H,H)+1E-10 in order to circumvent this issue.

Best, Emek
Could you avoid this by using a non-zero initial guess for u?
written 5 months ago by Nate  
0
5 months ago by
Hi Emek,

thanks for the little hint. Of course this was quite obvious ...


Best regards,

Dustin
Please login to add an answer/comment or follow this question.

Similar posts:
Search »