Hi everybody,
I am trying to solve a magnetostic problem with the Finite Element Method. But I have a problem applying tangential boundary conditions for the magentic field.
I solve for the vector potential using this equation:
$ \nabla \times (\frac{1}{\mu} \nabla \times \mathbf{A}) = \mu \mathbf{J} $
in 2d this reduces basically to the scalar laplace equation.
I know want to apply tangential boundary conditions, with mean:
$ \mathbf{n} * \mathbf{B} = \mathbf{n}* (\nabla \times \mathbf{A} ) = 0 $
I have the usual boundary integral:
$ \int_{\Gamma} N \left ( \alpha_{x} \frac{\partial U}{\partial x} n_{x} + \alpha_{y} \frac{\partial U}{\partial y} n_{y} \right ) dS $
But how to I incorporate my boundary conditions into this integral ?
Thanks for the help
Best Answer
$\newcommand{\v}[1]{\boldsymbol{#1}}$
First about the appropriate boundary condition we would like to impose: The cross product(scalar-$\mathrm{curl}$) in $\mathbb{R}^2$ is done by we embed a vector into $\mathbb{R}^3$: $\v{v} = \langle v_1,v_2,0\rangle$, then apply the cross product, what we get would be
$$\nabla \times \v{v} = \left\langle 0,0,\left|\begin{matrix}\partial_x&\partial_y\newline v_1&v_2 \end{matrix}\right|\right\rangle $$ Similarly the vector-$\mathbf{curl}$ is we embed a scalar function into the third dimension of $\mathbb{R}^3$: $u = \langle0,0,u\rangle$, and what we get is the rotation of the $\nabla$-operator normally denoted by $\nabla^{\perp}$: $$ \nabla^{\perp} u = \left\langle\frac{\partial u}{\partial y},-\frac{\partial u}{\partial x},0 \right\rangle $$ I read this way of defining 2D-curl operators in a famous computational fluid dynamics book Finite element methods for Navier-Stokes equations: theory and algorithms.
Then in $\mathbb{R}^2$ your problem can be written as: \begin{equation} \nabla^{\perp}(\mu^{-1}\nabla \times \v{A}) = \mu \v{J} \tag{1} \end{equation}
Suppose we multiply a test vector $\v{v}$ vanishing on boundary on both side or the equation and integration by parts:
$$ \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} - \int_{\partial \Omega} \mu^{-1}(\nabla \times \v{A})(\v{v} \times \v{n}) = \int_{\Omega} \mu\v{J}\cdot \v{v} $$ where recall in $\mathbb{R}^2$, $B = \mu^{-1}(\nabla \times \v{A})$ is really a scalar instead of vector, hence for the tangential Neumann boundary condition, doing the $\mathbb{R}^3$-embedding we have: $$ \v{n}\times B = \langle n_2 B, -n_1 B,0 \rangle $$ this is the value of the scalar $B$ along the tangential direction of the domain $\v{t} = \langle n_2, -n_1,0 \rangle$, dot product with $\v{t}$ itself we could get $B$($\v{t}$ is the rotation of a unit outward normal vector), hence the Neumann type boundary condition you really wanna impose here for 2D is: $$ \mu^{-1}(\nabla \times \v{A}) = g $$ where $g$ is a scalar function, and your variational form becomes: $$ \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} = \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) \tag{2} $$ The finite element method is we triangulating the domain, construct a finite dimensional Hilbert space $\v{V}_h$ (usually a subspace of the solution space $\v{V}$ to above variational problem), the boundary integral is treated as:
$$ \int_{\partial \Omega} g(\v{v} \cdot \v{t}) = \sum_{e\subset \partial \Omega} \int_{e} g (\v{v}\cdot \v{t}_e) $$
where $e$ is a boundary edge, and here is the place quadrature kicks in, that you wanna evaluate a line integral on each edge, it depends on what your finite element test function space $\v{V}_h$, nowadays a consensus is to choose Nédélec elements, in lowest order case it is $\v{v}_e = \lambda_A\nabla \lambda_B - \lambda_B\nabla \lambda_A $ for $e = \overrightarrow{AB}$.
Secondly, what is really problematic with your statement here is: for a smooth vector field $\v{A} \in \v{C}^{\infty}(\mathbb{R}^2)$ or $\v{C}^{\infty}_c(\Omega)$, the equation (1) doesn't simply reduce to the scalar Laplace equation as you mentioned in OP, a direct observation would be that homogeneous Neumann problem for $-\Delta u = f$ has a constant as null space because a solution plus a constant is still a solution, while $\v{curl}$-problem has all the gradient vector fields as null space!
Similar thing applies for variational formulation (2), the bilinear form is not coercive, simply plugging a gradient field $\v{A} = \v{v} = \nabla \varphi$ into it. What we numerical analyst usually do here is to find a solution in the quotient space $\v{H}(\mathrm{curl})/\nabla H^1$ by adding a Lagrange multiplier to let the problem fall into mixed finite element(google MFEM if you would like to know more) framework:
$$ \begin{cases} \int_{\Omega}\mu^{-1} \nabla\times \v{A} \cdot \nabla \times \v{v} + \int_{\Omega} \v{v}\cdot \nabla \varphi= \int_{\Omega} \mu\v{J}\cdot \v{v} + \int_{\partial \Omega} g(\v{v} \cdot \v{t}) & \text{ for } \forall \v{v}\in \v{H}(\mathrm{curl})\newline \int_{\Omega} \v{A} \cdot \nabla \psi = 0& \text{ for } \forall \psi\in H^1/\mathbb{R} \end{cases} $$ to make sure problem (2) is solvable using finite element, not only we have to use Lagrange multiplier to circumvent the "null space" problem, we also have to guarantee that the right side is in the range of the curl-operator, here plug $\v{v} = \nabla \varphi$ into (2): $$ 0=\int_{\Omega} \mu\v{J}\cdot \nabla \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) = -\int_{\Omega} \nabla \cdot(\mu\v{J}) \varphi + \int_{\partial \Omega} (\mu\v{J} \cdot \v{n}) \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \v{t}) $$ in your case, $g=0$, also by the arbitrariness of $\varphi$, for eg we choose $\nabla \varphi \cdot \v{t} = 0$ for $\varphi\in H^1_0$, we would get so called "compatibility condition" in finite element method for Neumann problem: $$ \nabla \cdot(\mu\v{J}) = 0\text{ in }\Omega\quad \text{ and }\quad \mu\v{J}\cdot \v{n} = 0\text{ on }\partial\Omega $$
Lastly about the implementation of finite element method for this equation, if you have a pure Neumann problem with homogeneous boundary condition(derivative on boundary is zero), then you do not have to do anything during assembling of the stiffness matrix or right hand side integral.