The boundary condition of Hall-MHD model

General Discussion
1

Hello, everyone.

Now I’m using freefem++ to solve the incompressible hall-MHD equations, but I’ve encountered a problem regarding the boundary conditions of the magnetic field.

\mu\nabla\times (\nabla\times\mathbf{B})+(\mathbf{u}^n\cdot\nabla)\mathbf{B}^n-(\mathbf{B}^n\cdot\nabla)\mathbf{u}^n+\eta\nabla\times\left((\nabla\times\mathbf{B}^n)\times\mathbf{B}^n\right)=\mathbf{g}
\nabla\cdot\mathbf{B}=0
Before solving, assume that the velocity field \mathbf{u}^n,\mathbf{B}^n is already known (I decoupled the equation).

I use the boundary condition of magnetic
\mathbf{B}\times \mathbf{n} = \mathbf{0} or \mathbf{B}\times \mathbf{n} = \mathbf{B}_e , where \mathbf{B}_e is the analytical solution that I set, and \mathbf{n} is the unit external normal vector of the boundary. The domain is [0,1]^3.

solve Magnetic(Bx1, By1, Bz1, wx, wy, wz, solver = GMRES) 
        = int3d(Th)( mu*curl(Bx1, By1, Bz1)'*curl(wx, wy, wz) + mu*div(Bx1, By1, Bz1)*div(wx, wy, wz) )
        + int3d(Th)( UgradV(ux0, uy0, uz0, Bx0, By0, Bz0)'*[wx, wy, wz] - UgradV(Bx0, By0, Bz0, ux0, uy0, uz0)'*[wx, wy, wz] 
            + eta*crossProduct(dy(Bz0)-dz(By0), dz(Bx0)-dx(Bz0), dx(By0)-dy(Bx0), Bx0, By0, Bz0)'*curl(wx, wy, wz) 
            - (g1t*wx + g2t*wy + g3t*wz))
    ;

with

    macro Grad(ux, uy, uz) [(dx(ux)), (dy(ux)), (dz(ux)), (dx(uy)), (dy(uy)), (dz(uy)), (dx(uz)), (dy(uz)), (dz(uz))] //
    macro div(ux, uy, uz) ( dx(ux) + dy(uy) + dz(uz) ) // 
    macro UgradV(ux, uy, uz, vx, vy, vz) [(ux)*(dx(vx)) + (uy)*(dy(vx)) + (uz)*(dz(vx)), (ux)*(dx(vy)) + (uy)*(dy(vy)) + (uz)*(dz(vy)), (ux)*(dx(vz)) + (uy)*(dy(vz)) + (uz)*(dz(vz))] //
    macro curl(Bx, By, Bz) [dy(Bz) - dz(By), dz(Bx) - dx(Bz), dx(By) - dy(Bx)] //
    macro crossProduct(ux, uy, uz, vx, vy, vz) [(uy)*(vz) - (uz)*(vy), (uz)*(vx) - (ux)*(vz), (ux)*(vy) - (uy)*(vx)] //

But, the result is wrong. Maybe I should add a stabilization term, but I do not how to do

2

To be consistent with the boundary condition B\times n=0 you should add a penalty boundary term
+int2d(Th)(1.e15*crossProduct(Bx1, By1, Bz1, N.x, N.y, N.z)'*crossProduct(wx, wy, wz, N.x, N.y, N.z))
For a non homogeneous condition B\times n = B_e you should add also
-int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z))

3

I understand. But is there any requirement for the penalty parameter? I see you gave 1.e15

4

The value is not important, it needs just to be much larger than the other terms, from 1.e8 to 1.e30 should do the job. The standards of FFEM uses tgv=1.e30.

5

OK, I understan.

Finally, I need to solve a Poisson equation for pressure

solve poisson(p, pp, solver = CG) = int3d(Th)( dx(p)*dx(pp) + dy(p)*dy(pp) + dz(p)*dz(pp) )
        - int3d(Th)( (dx(p0)*dx(pp) + dy(p0)*dy(pp) + dz(p0)*dz(pp)) - div(ux, uy, uz)*pp/dt )
        + int2d(Th)( 1.e-8*p*pp );

Can I set the parameter to 1.e-8?

6

Yes it should work with 1.e8 (not 1.e-8). Alternatively you can use the built-in command
+on(1,2,3,4,5,6,p=0)
1,2,3,4,5,6 refers to the boundary labels. This built-in command is fully equivalent to
+int2d(Th,1,2,3,4,5,6)(tgv*p*pp)

7

It should be 1e-8 (I try to set 1.e8, but it is wrong)

截屏2024-09-14 20.58.49
截屏2024-09-14 20.58.49968×140 5.54 KB

This is my equation, which has the artificial Newman boundary condition.

Now there is another problem.The convergence order of the L2 norm of the magnetic field I calculated is not correct (as the space is divided finer, the order decreases), but the convergence order of the H1 semi norm is no problem. I don’t know if this is due to the code.

h = [1./2, 1./4, 1./8, 1./12, 1./14, 1./16, 1./18, 1./20, 1./22, 1./24];

Convergence B_L2: 0     2.937716417     3.003305649     2.935914756     2.751860597
        2.933047063     2.438786889     1.90737733      4.622726824     0.3392563845
Convergence B_H1_semi: 0     1.885937569     1.974759866     1.985809507     1.985307756
        1.983228133     1.998874224     1.943471302     2.039907908     2.000346527
8

The 1.e8 is for homogeneous Dirichlet boundary condition. If you want Neumann, you just don’t put the penalty term!
About the order of convergence I don’t know, but your results look pretty good.

9

Dear Professor,
For the case non homegenous case B \times n=B_e, why should add both

+int2d(Th)(1.e15*crossProduct(Bx1, By1, Bz1, N.x, N.y, N.z)'*crossProduct(wx, wy, wz, N.x, N.y, N.z))

-int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z)) ?
Why is adding only -int2d(Th)(1.e15*[Bex,Bey,Bez]'*crossProduct(wx, wy, wz, N.x, N.y, N.z)) not enough?

10

Because both terms together give
10^{15}*\int_\Gamma(B\times N-B_e)\cdot(w\times N)
Hence putting this in the variational formulation implies that
\int_\Gamma(B\times N-B_e)\cdot(w\times N)\simeq 0
for all test w.
We have to assume that B_e\cdot N=0 otherwise the boundary condition is impossible to satisfy. Then since w\times N is an arbitrary vector on the boundary, orthogonal to N, it gives that B\times N-B_e\simeq 0 on \Gamma.
After that it remains the variational formulation (interior terms) for all test w
such that w\times N=0 on \Gamma.

11

Dear Professor,
What should we do for the case n \times B =n \times B_e?

12

Just replace B_e by B_e\times N in the code!