Mixed finite element methods for the ferrofluid model with magnetization paralleled to the magnetic field

Yongke Wu and Xiaoping Xie

Abstract.

Mixed finite element methods are considered for a ferrofluid flow model with magnetization paralleled to the magnetic field. The ferrofluid model is a coupled system of the Maxwell equations and the incompressible Navier-Stokes equations. By skillfully introducing some new variables, the model is rewritten as several decoupled subsystems that can be solved independently. Mixed finite element formulations are given to discretize the decoupled systems with proper finite element spaces. Existence and uniqueness of the mixed finite element solutions are shown, and optimal order error estimates are obtained under some reasonable assumptions. Numerical experiments confirm the theoretical results.

Key words and phrases:

ferrofluid flow, decoupled system, mixed finite element method, error estimate

2010 Mathematics Subject Classification:

65N55; 65F10; 65N22; 65N30;

Y. Wu was supported by the National Natural Science Foundation of China (11971094 and 11501088), and X. Xie was supported by the National Nature Science Foundation of China (12171340)

Y. Wu, School of Mathematical Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China. Email: [email protected]

X. Xie (Corresponding author), School of Mathematical, Sichuan University, Chengdu 610064, China. Email: [email protected]

1. Introduction

Ferrofluids are colloidal liquids consisting of non-conductive nanoscale ferromagnetic or ferrimagnetic particles suspended in carrier fluids, and have extensive applications in many technology and biomedicine fields [22, 31]. There are two main ferrohydrodynamics (FHD) models which treat ferrofluids as homogeneous monophase fluids: the Rosensweig’s model [24, 23] and the Shliomis’ model [25, 26]. The main difference between these two models lies in that the former one considers the internal rotation of the nanoparticles, while the latter deals with the rotation as a magnetic torque. We refer to [1, 4, 2, 3, 21] for some existence results of solutions to the two FHD models.

The FHD models are coupled nonlinear systems of the Maxwell equations and the incompressible Navier-Stokes equations. There are limited works in the literature on the numerical analysis of the FHD models. In [27, 17, 16, 30], several numerical schemes were applied to solve reduced two dimensional FHD models where some nonlinear terms of the original models are dropped. Nochetto et al. [20] showed the energy stability of the Rosensweig’s model, proposed an energy-stable numerical scheme using finite element methods, and gave the existence and convergence of the numerical solutions. Recently, Wu and Xie [29] developed a class of energy-preserving mixed finite element methods for the Shliomis’ FHD model, and derived optimal error estimates for both the the semi- and fully discrete schemes. We also note that in [32] an unconditionally energy-stable fully discrete finite element method was presented for a two-phase FHD model.

In this paper, we consider the Shliomis’ FHD model with the assumption that the magnetization field is parallel to the magnetic field. Under this assumption, the magnetization equation in the Shliomis’ model degenerates to the Langevin magnetization law [17, 24, 23]. We introduce some new variables to transform the model into two main decoupled subsystems, i.e., a nonlinear elliptic equation and the incompressible Navier-Stokes equations. We apply proper finite element spaces to discretize the nonlinear decoupled systems, prove the existence and, under some reasonable assumptions, uniqueness of the finite element solutions, and derive optimal error estimates. We also show that our scheme preserves the ferrofluids’ nonconductive property ${\rm curl\,}\boldsymbol{H}=0$ exactly.

The rest of this paper is organized as follows. In section 2, we introduce several Sobolev spaces, give the governing equation of the FHD model with magnetization paralleled to magnetic field, reform the FHD model equivalently, and construct the weak formulations. In section 3, we recall the finite element spaces, show the existence and uniqueness of solutions for the constructed finite element methods, and give the optimal order error estimates. In section 4, some numerical experiments will be given to verify our theoretical results.

2. Preliminary

In this section, we introduce several Sobolev spaces, give the governing equations of the FHD model with magnetization paralleled to magnetic field, derive the equivalent decoupled systems, and present the weak formulations.

2.1. Sobolev spaces

Let $\Omega\subset\mathbb{R}^{d}$ ( $d=2,3$ ) be a bounded and simply connected convex domain with Lipschitz boundary $\partial\Omega$ , and $\boldsymbol{n}$ be the unit outward normal vector on $\partial\Omega$ .

For any $p\geq 1$ , we denote by $L^{p}(\Omega)$ the space of all power- $p$ integrable functions on $\Omega$ with norm $\|\cdot\|_{L^{p}}$ . For any nonnegative integer $m$ , denote by $H^{m}(\Omega)$ the usual $m$ -th order Sobolev space with norm $\|\cdot\|_{m}$ and semi-norm $|\cdot|_{m}$ . In particular, $H^{0}(\Omega)=L^{2}(\Omega)$ denotes the space of square integrable functions on $\Omega$ , with the inner product $(\cdot,\cdot)$ and norm $\|\cdot\|$ . For the vector spaces $\boldsymbol{H}^{m}(\Omega):=(H^{m}(\Omega))^{d}$ and $\boldsymbol{L}^{2}(\Omega):=(L^{2}(\Omega))^{d}$ , we use the same notations of norm, semi-norm and inner product as those for the scalar cases.

We further introduce the Sobolev spaces

\boldsymbol{H}({\rm curl\,}):=\left\{\begin{array}[]{ll}\{\boldsymbol{v}\in(L^{2}(\Omega))^{2}:~{}~{}{\rm curl\,}\boldsymbol{v}\in L^{2}(\Omega)\}&\text{if }d=2,\\ \{\boldsymbol{v}\in(L^{2}(\Omega))^{3}:~{}~{}{\rm curl\,}\boldsymbol{v}\in(L^{2}(\Omega))^{3}\}&\text{if }d=3\end{array}\right.

and

\boldsymbol{H}(\operatorname{div}):=\{\boldsymbol{v}\in(L^{2}(\Omega))^{d}:~{}~{}\operatorname{div}\boldsymbol{v}\in L^{2}(\Omega)\},

and set

	$\displaystyle H_{0}^{1}(\Omega):=\{v\in H^{1}(\Omega):~{}v=0\text{ on }\partial\Omega\},$
	$\displaystyle\boldsymbol{H}_{0}({\rm curl\,}):=\{\boldsymbol{v}\in\boldsymbol{H}({\rm curl\,}):~{}\boldsymbol{v}\times\boldsymbol{n}=0\text{ on }\partial\Omega\},$
	$\displaystyle\boldsymbol{H}_{0}(\operatorname{div}):=\{\boldsymbol{v}\in\boldsymbol{H}(\operatorname{div}):~{}\boldsymbol{v}\cdot\boldsymbol{n}=0\text{ on }\partial\Omega\},$
	$\displaystyle L_{0}^{2}(\Omega):=\{v\in L^{2}(\Omega):~{}\int_{\Omega}v\,{\rm d}\Omega=0\},$

where

	$\displaystyle{\rm curl\,}\boldsymbol{v}:=\left\{\begin{array}[]{ll}\partial_{x}v_{2}-\partial_{y}v_{1}&\text{if }\boldsymbol{v}=(v_{1},~{}v_{2})^{\intercal},\\ (\partial_{y}v_{3}-\partial_{z}v_{2},~{}\partial_{z}v_{1}-\partial_{x}v_{3},~{}\partial_{x}v_{2}-\partial_{y}v_{1})^{\intercal}&\text{if }\boldsymbol{v}=(v_{1},~{}v_{2},~{}v_{3})^{\intercal}.\end{array}\right.$
	$\displaystyle\operatorname{div}\boldsymbol{v}:=\left\{\begin{array}[]{ll}\partial_{x}v_{1}+\partial_{y}v_{2}&\text{if }\boldsymbol{v}=(v_{1},v_{2})^{\intercal},\\ \partial_{x}v_{1}+\partial_{y}v_{2}+\partial_{z}v_{3}&\text{if }\boldsymbol{v}=(v_{1},v_{2},v_{3})^{\intercal},\end{array}\right.$

and $(x,y)$ and $(x,y,z)$ are the Cartesian coordinates in two and three dimensions, respectively.

For any Sobolev space $S$ with norm $\|\cdot\|_{S}$ , we use $S^{\prime}$ to denote the dual space of $S$ , and $\langle\cdot,\cdot\rangle$ to denote the dual product between $S^{\prime}$ and $S$ . For any $f\in S^{\prime}$ , the operator norm of $f$ is defined as $\|f\|_{S^{\prime}}=\sup\limits_{0\neq v\in S}\frac{\langle f,v\rangle}{\|v\|_{S}}$ .

2.2. Governing equations of the ferrofluid flow

We consider the domain $\Omega$ filled with ferrofluid flow. On the macroscopic level, the mathematical model for describing the interactions between magnetic fields and ferrofluids consists of the Maxwell’s equations and the Navier-Stokes equations [23, 24, 25, 26]. Since the ferrofluid flow is nonconductive, the corresponding Maxwell’s equations read as:

(1)			$\displaystyle{\rm curl\,}\boldsymbol{H}=0\qquad\qquad\qquad\,\text{in}\quad\Omega,$
(2)			$\displaystyle\operatorname{div}\boldsymbol{B}=\operatorname{div}\boldsymbol{H}_{e}\qquad\qquad\text{in}\quad\Omega,$

with the magnetic field $\boldsymbol{H}$ and the magnetic induction $\boldsymbol{B}$ satisfying the relation

(3)

\boldsymbol{B}=\mu_{0}(\boldsymbol{H}+\boldsymbol{M})\qquad\qquad\text{in}\quad\Omega,

where $\mu_{0}>0$ is the permeability constant, $\boldsymbol{M}$ is the magnetization, and $\boldsymbol{H}_{e}$ is the known external magnetic field that satisfies $\boldsymbol{H}_{e}\cdot\boldsymbol{n}=0$ on $\partial\Omega$ .

Under the assumption that the magnetization $\boldsymbol{M}$ of the ferrofluid flow is parallel to the magnetic field $\boldsymbol{H}$ , it follows the nonlinear Langevin magnetization law [23, 24, 17], i.e.,

(4)

\boldsymbol{M}(\boldsymbol{H})=M_{s}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{\boldsymbol{H}}{H},

with the saturation magnetization $M_{s}>0$ , the Lagevin parameter $\gamma=3\chi_{0}/M_{s}$ , the initial susceptibility $\chi_{0}>0$ , $H:=|\boldsymbol{H}|=(\boldsymbol{H}\cdot\boldsymbol{H})^{1/2}>0$ , and $\coth(x)=\frac{e^{x}+e^{-x}}{e^{x}-e^{-x}}$ .

The hydrodynamic properties of the ferrofluid flow are described by the incompressible Navier-Stokes equations

(5)			$\displaystyle\rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\eta\Delta\boldsymbol{u}+\nabla p=\boldsymbol{f}+\mu_{0}(\boldsymbol{M}\cdot\nabla)\boldsymbol{H}\qquad\quad\text{in}\quad\Omega,$
(6)			$\displaystyle\operatorname{div}\boldsymbol{u}=0\,\,\,\,\,\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\text{in}\quad\Omega,$

where $\rho$ denotes the fluid density, $\boldsymbol{u}$ the velocity field of the flow, $\eta$ the dynamic viscosity, and $\boldsymbol{f}$ the known volume force.

We consider the following homogenous boundary conditions for equation (1)-(6):

(7)

\boldsymbol{u}=0\quad\text{ and }\quad\boldsymbol{H}\times\boldsymbol{n}=0\qquad\text{on}\quad\partial\Omega.

Using the fact that

(\boldsymbol{H}\cdot\nabla)\boldsymbol{H}=\frac{1}{2}\nabla H^{2}-\boldsymbol{H}\times({\rm curl\,}\boldsymbol{H})=\frac{1}{2}\nabla H^{2}

and the Langevin magnetization law (4), we have

(8)

(\boldsymbol{M}\cdot\nabla)\boldsymbol{H}=\frac{M_{s}}{H}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{1}{2}\nabla H^{2}=\frac{M_{s}}{\gamma}\nabla\ln\frac{\sinh(\gamma H)}{H},

where $\sinh(x)=\frac{e^{x}-e^{-x}}{2}$ . Denote

(9)

\beta(x):=\frac{M_{s}}{\gamma}\ln\frac{\sinh(\gamma x)}{x},

and introduce two variables

(10)

\psi:=\beta(H)

and

(11)

\tilde{p}:=p-\mu_{0}\psi,

then equation (5) becomes

(12)

\rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\eta\Delta\boldsymbol{u}+\nabla\tilde{p}=\boldsymbol{f}\qquad\qquad\qquad\text{in }\Omega.

Equation (1) and the boundary condition $\boldsymbol{H}\times\boldsymbol{n}|_{\partial\Omega}=0$ in (7) imply that (cf. [6]) there exists $\phi\in H_{0}^{1}(\Omega)$ with

(13)

\boldsymbol{H}=\nabla\phi,

then combining (2), (3) and (4) leads to

(14)

\nabla\cdot\left(\alpha(|\nabla\phi|)\nabla\phi\right)=\frac{1}{\mu_{0}}\operatorname{div}\boldsymbol{H}_{e}=:g\qquad\text{in}~{}~{}\Omega,

with

(15)

\alpha(x):=1+\frac{M_{s}}{x}\left(\coth(\gamma x)-\frac{1}{\gamma x}\right).

The above equivalence transformation shows that the FHD model (1) - (6) with the boundary conditions (7) can be equivalently written as follows: Find $\phi\in H_{0}^{1}(\Omega)$ , $\boldsymbol{H}\in\boldsymbol{H}_{0}({\rm curl\,})$ , $\boldsymbol{M}\in\boldsymbol{H}_{0}({\rm curl\,})$ , $\boldsymbol{u}\in\boldsymbol{H}_{0}^{1}(\Omega)$ , $\tilde{p}\in L_{0}^{2}(\Omega)$ , $\psi\in L_{0}^{2}(\Omega)$ , and $p\in L_{0}^{2}(\Omega)$ such that

(16)

\left\{\begin{array}[]{ll}\nabla\cdot(\alpha(|\nabla\phi|)\nabla\phi)=g&\text{in }\Omega,\\ \boldsymbol{H}-\nabla\phi=0&\text{in }\Omega,\\ \boldsymbol{M}=M_{s}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{\boldsymbol{H}}{H}=\big{(}\alpha(H)-1\big{)}\boldsymbol{H}&\text{in }\Omega,\\ \rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\eta\Delta\boldsymbol{u}+\nabla\tilde{p}=\boldsymbol{f}&\text{in }\Omega,\\ \nabla\cdot\boldsymbol{u}=0&\text{in }\Omega,\\ \psi=\frac{M_{s}}{\gamma}\ln\frac{\sinh(\gamma H)}{H}=\beta(H)&\text{in }\Omega,\\ p=\tilde{p}+\mu_{0}\psi&\text{in }\Omega.\end{array}\right.

Remark 2.1.

The system (16) is a decoupled system. Firstly, we can solve the first nonlinear elliptic equation of (16) to get $\phi$ , and solve the Navier-Stokes equations, i.e. the fourth and fifth equations, to get $\boldsymbol{u}$ and $\tilde{p}$ . Secondly, we can get $\boldsymbol{H}$ from the second equation of (16). Finally, we can obtain $\boldsymbol{M}$ , $\psi$ , and $p$ from the third, the sixth, and the seventh equations, respectively.

2.3. Preliminary estimates for nonlinear functions

Notice that the decoupled system (16) involves the nonlinear functions $\alpha(x)$ and $\beta(x)$ defined in (15) and (9), respectively. The following basic estimates of these two functions will be used in later analysis.

Lemma 2.2.

(i)

There exist a positive constants $C_{1}$ and $C_{\alpha}$ such that for any $x>0$ , there hold

$1<\alpha(x)\leq C_{1},\qquad|\alpha^{\prime}(x)|\leq C_{\alpha}\ ;$
(ii)

For any $x>0$ , there holds

$|\beta^{\prime}(x)|\leq M_{s}.$

Proof.

(i) We first show $1<\alpha(x)\leq C_{1}.$ On one hand, the L’Hopital law implies that

\lim\limits_{y\rightarrow 0^{+}}y\coth(y)=1,

which, together with the fact that

(y\coth(y))^{\prime}=\coth(y)-y{\rm csch}^{2}(y)=\frac{e^{2y}-e^{-2y}-4y}{(e^{y}-e^{-y})^{2}}>0\quad\forall~{}~{}y>0,

shows

y\coth(y)>1\qquad\forall~{}~{}y>0.

Therefore,

\alpha(x)=1+\gamma M_{s}\frac{y\coth(y)-1}{y^{2}}>1\qquad\text{with }y=\gamma x,\text{ for }\forall~{}~{}x>0.

On the other hand, we easily see that

\lim\limits_{x\rightarrow 0^{+}}\alpha(x)=1+\frac{\gamma M_{s}}{3}\qquad\text{and}\qquad\lim\limits_{x\rightarrow+\infty}\alpha(x)=1,

which, together the continuity of $\alpha(x)$ for $x>0$ , indicate that there exists a positive constant $C_{1}$ such that

\alpha(x)\leq C_{1}\qquad\forall~{}~{}x>0.

Second, let us show $|\alpha^{\prime}(x)|\leq C_{\alpha}.$ It is easy to get

\begin{split}\alpha^{\prime}(x)&=\frac{2M_{s}}{\gamma x^{3}}-\frac{M_{s}}{x^{2}}\coth(\gamma x)-\frac{\gamma M_{s}}{x}{\rm csch}^{2}(\gamma x)\\ &=\frac{\gamma^{2}M_{s}}{y^{3}}(2-y^{2}{\rm csch}^{2}(y)-y\coth(y))\qquad\qquad\text{with }y=\gamma x.\end{split}

The L’Hopital law implies that

\lim\limits_{x\rightarrow 0^{+}}\alpha^{\prime}(x)=0\qquad\text{and}\qquad\lim\limits_{x\rightarrow+\infty}\alpha^{\prime}(x)=0.

Since $\alpha^{\prime}(x)$ is continuous for $x>0$ , we know that there exists a positive constant $C_{\alpha}$ such that

|\alpha^{\prime}(x)|\leq C_{\alpha}\qquad\forall x>0.

As a result, the conclusions of (i) follow.

(ii) It is easy to find that

\beta^{\prime}(x)=M_{s}\left(\coth(\gamma x)-\frac{1}{\gamma x}\right)=M_{s}\left(\coth(y)-\frac{1}{y}\right)

with $y=\gamma x>0$ . The L’Hopital law implies that

\lim\limits_{x\rightarrow 0^{+}}\beta^{\prime}(x)=0\quad\text{and}\quad\lim\limits_{x\rightarrow+\infty}\beta^{\prime}(x)=M_{s}.

The fact that

\left(\coth(y)-\frac{1}{y}\right)^{\prime}=\frac{1-y^{2}{\rm csch}^{2}(y)}{y^{2}}=\frac{1+y{\rm csch}(y)}{y^{2}}\frac{e^{y}-e^{-y}-2y}{e^{y}-e^{-y}}>0\qquad\forall~{}~{}y>0

implies that $\beta^{\prime}(x)$ is a monotonically increasing function on the interval $(0,+\infty)$ . We conclude that

0<\beta^{\prime}(x)\leq M_{s}.

This finishes the proof. ∎

2.4. Weak formulations

Based on the FHD model (16) and Remark 2.1, we consider the following weak formulations: Find $\phi\in H_{0}^{1}(\Omega)$ , $\boldsymbol{H}\in\boldsymbol{H}_{0}({\rm curl\,})$ , $\boldsymbol{u}\in\boldsymbol{H}_{0}^{1}(\Omega)$ and $\tilde{p}\in L_{0}^{2}(\Omega)$ such that

(17)			$\displaystyle a(\phi;\phi,\tau)=-(g,\tau)\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \,\forall~{}~{}\tau\in H_{0}^{1}(\Omega),$
(18)			$\displaystyle(\boldsymbol{H},\boldsymbol{C})-(\nabla\phi,\boldsymbol{C})=0\quad\qquad\qquad\qquad\qquad\qquad\qquad\forall~{}~{}\boldsymbol{C}\in\boldsymbol{H}_{0}({\rm curl\,}),$
(19)			$\displaystyle(\rho(\boldsymbol{u}\cdot\nabla)\boldsymbol{u},\boldsymbol{v})+\eta(\nabla\boldsymbol{u},\nabla\boldsymbol{v})-(\tilde{p},\nabla\cdot\boldsymbol{v})=(\boldsymbol{f},\boldsymbol{v})\qquad\forall~{}~{}\boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega),$
(20)			$\displaystyle(\nabla\cdot\boldsymbol{u},q)=0\,\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\forall~{}~{}q\in L_{0}^{2}(\Omega),$

where $a(\cdot;\cdot,\cdot):~{}H^{1}(\Omega)\times H^{1}(\Omega)\times H^{1}(\Omega)\rightarrow\mathbb{R}$ is defined by

a(w;\phi,\tau):=\int_{\Omega}\alpha(|\nabla w|)\nabla\phi\cdot\nabla\tau\,{\rm d}\Omega\qquad\forall~{}~{}w,~{}\phi,~{}\tau\in H^{1}(\Omega).

Remark 2.3.

As shown in Remark 2.1, once the variables $\phi$ , $\boldsymbol{H}$ , $\boldsymbol{u}$ and $\tilde{p}$ are solved, the other variables, i.e. $\boldsymbol{M}$ , $\psi$ and $p$ , can immediately be obtained by the third, the sixth and the seventh equations of (16), respectively.

In what follows, we discuss the existence and uniqueness of the solutions to the weak formulations (17) - (20).

We first consider the nonlinear equation (17). It is easy to see that, for any given $w\in H^{1}(\Omega)$ , $a(w;\cdot,\cdot)$ is a bilinear form on $H^{1}(\Omega)$ . Recall that the Poincaré inequality

\|v\|\leq C_{p}\|\nabla v\|\qquad\forall~{}~{}v\in H_{0}^{1}(\Omega),

with $C_{p}>0$ being a constant depending only on $\Omega$ , means that the semi-norm $\|\nabla(\cdot)\|$ is also a norm on $H_{0}^{1}(\Omega)$ . Then, from Lemma 2.2 (i), we easily obtain the following uniform coercivity and continuity results for $a(w;\cdot,\cdot)$ :

Lemma 2.4.

For any $w\in H_{0}^{1}(\Omega)$ and $\phi,\tau\in H_{0}^{1}(\Omega)$ , we have

a(w;\phi,\phi)\geq\|\nabla\phi\|^{2}

and

a(w;\phi,\tau)\leq C_{1}\|\nabla\phi\|\|\nabla\tau\|.

We have the following wellposedness results for equation (17).

Theorem 2.5.

The nonlinear equation (17) has at least one solution $\phi\in H_{0}^{1}(\Omega)$ , and there holds

(21)

\|\nabla\phi\|\leq\frac{1}{\mu_{0}}\|\boldsymbol{H}_{e}\|.

Moreover, (17) admits at most one solution $\phi\in H_{0}^{1}(\Omega)$ satisfying

(22)

\|\nabla\phi\|_{L^{\infty}}<1/C_{\alpha}.

Proof.

The existence of a solution $\phi\in H_{0}^{1}(\Omega)$ follows from Lemma 2.4 and [8, Page 332, Theorem 2] directly.

Recalling that $g=\frac{1}{\mu_{0}}\operatorname{div}\boldsymbol{H}_{e}$ , from (17) and the above lemma we have

\|\nabla\phi\|^{2}\leq a(\phi;\phi,\phi)=-(g,\phi)=\frac{1}{\mu_{0}}(\boldsymbol{H}_{e},\nabla\phi)\leq\frac{1}{\mu_{0}}\|\boldsymbol{H}_{e}\|\|\nabla\phi\|.

This yields the estimate (21).

The thing left is to show the uniqueness of the solution under condition (22). In fact, let $\phi$ and $\tilde{\phi}$ be any two solutions of (17) satisfying (22), we have

a(\tilde{\phi};\tilde{\phi},\tau)=a(\phi;\phi,\tau)\qquad\forall~{}~{}\tau\in H_{0}^{1}(\Omega),

which yields

\displaystyle a(\tilde{\phi};\tilde{\phi}-\phi,\tau)=\int_{\Omega}\big{(}\alpha(|\nabla\phi|)-\alpha(|\nabla\tilde{\phi}|)\big{)}\nabla\phi\cdot\nabla\tau\,{\rm d}\Omega\qquad\forall~{}~{}\tau\in H_{0}^{1}(\Omega).

Taking $\tau=\tilde{\phi}-\phi$ in the above equation and using Lemma 2.2 give

	$\displaystyle\\|\nabla(\tilde{\phi}-\phi)\\|^{2}$	$\displaystyle\leq a(\tilde{\phi};\tilde{\phi}-\phi,\tilde{\phi}-\phi)$
		$\displaystyle=a(\phi;\phi,\tilde{\phi}-\phi)-a(\tilde{\phi};\phi,\tilde{\phi}-\phi)$
		$\displaystyle\leq C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}\\|\nabla(\tilde{\phi}-\phi)\\|^{2},$

which, together with (22), implies $\tilde{\phi}=\phi$ . This finishes the proof. ∎

For equation (18), we have the following conclusion:

Theorem 2.6.

For any given $\phi\in H_{0}^{1}(\Omega)$ , equation (18) admits a unique solution $\boldsymbol{H}=\nabla\phi\in\boldsymbol{H}_{0}({\rm curl\,})$ , which means ${\rm curl\,}\boldsymbol{H}=0.$

Proof.

The de Rham complex [6] on $2D$ and $3D$ domain implies that $\nabla\phi\in\boldsymbol{H}_{0}({\rm curl\,})$ . Thus, taking $\boldsymbol{C}=\boldsymbol{H}-\nabla\phi\in\boldsymbol{H}_{0}({\rm curl\,})$ in (18) yields $\boldsymbol{H}=\nabla\phi.$ ∎

Since (19)-(20) are the weak formulations of the Navier-Stokes equations, the following existence and uniqueness results are standard (cf. [13, Page 285-287, Theorems 2.1 and 2.2]).

Theorem 2.7.

Given $\boldsymbol{f}\in(H^{-1}(\Omega))^{d}$ , there exists at least one pair $(\boldsymbol{u},\tilde{p})\in\boldsymbol{H}_{0}^{1}(\Omega)\times L_{0}^{2}(\Omega)$ satisfying (19) - (20), and holds

\|\nabla\boldsymbol{u}\|\leq\frac{1}{\eta}\|\boldsymbol{f}\|_{-1},

where $\|\boldsymbol{f}\|_{-1}:=\sup\limits_{\boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega)}\frac{\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{v}\,{\rm d}\Omega}{\|\nabla\boldsymbol{v}\|}.$ In addition, if

(\rho\mathcal{N}/\eta^{2})\|\boldsymbol{f}\|_{-1}<1\quad\text{with}\quad\mathcal{N}:=\sup\limits_{\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}\in\boldsymbol{H}_{0}^{1}(\Omega)}\frac{\int_{\Omega}(\boldsymbol{u}\cdot\nabla)\boldsymbol{v}\cdot\boldsymbol{w}\,{\rm d}\Omega}{\|\nabla\boldsymbol{u}\|\|\nabla\boldsymbol{v}\|\|\nabla\boldsymbol{w}\|},

then the solution pair $(\boldsymbol{u},\tilde{p})$ is unique.

3. Finite element methods

3.1. Finite element spaces

Let $\mathcal{T}_{h}$ be a quasi-uniform simplicial decomposition of $\Omega$ with mesh size $h:=\max\limits_{K\in\mathcal{T}_{h}}h_{K}$ , where $h_{K}$ denotes the diameter of $K$ for any $K\in\mathcal{T}_{h}$ .

For an integer $l\geq 0$ , let $\mathbb{P}_{l}(K)$ denote the set of polynomials, defined on $K\in\mathcal{T}_{h}$ , of degree no more than $l$ . We introduce the following finite element spaces:

•
$\boldsymbol{V}_{h}=(V_{h})^{d}\ \big{(}\subset\text{or }\not\subset(H_{0}^{1}(\Omega))^{d}\big{)}$ is the Lagrange element space [13, 9, 14] for the velocity $\boldsymbol{u}$ , with $\mathbb{P}_{l+1}(K)\subset V_{h}|_{K}$ for any $K\in\mathcal{T}_{h}$ . In particular, for the nonconforming case that $V_{h}\not\subset H_{0}^{1}(\Omega)$ , each $v_{h}\in V_{h}$ is required to satisfy the following conditions:
- (i)
  
  $v_{h}$ vanishes at all the nodes on $\partial\Omega$ ;
- (ii)
  
  $|v_{h}|_{1,h}:=(\sum\limits_{K\in\mathcal{T}_{h}}\|\nabla v_{h}\|_{K}^{2})^{1/2}$ is a norm.
Note that the classical nonconforming Crouzeix-Raviart (CR) finite element [12] is corresponding to the nonconforming case with $l=0$ . In the nonconforming case, the gradient and divergence operators, $\nabla$ and $\nabla\cdot$ , in the finite element scheme (28) - (29) will be understood as $\nabla_{h}$ and $\nabla_{h}\cdot$ , respectively, which denote respectively the piecewise gradient and divergence operators acting on element by element in $\mathcal{T}_{h}$ .
•

$W_{h}\subset L_{0}^{2}(\Omega)$ is the piecewise polynomial space for the ‘pressure’ variable $\tilde{p}$ , with $\mathbb{P}_{l}(K)\subset W_{h}|_{K}$ for any $K\in\mathcal{T}_{h}$ .
•

$S_{h}\subset H_{0}^{1}(\Omega)$ is the Lagrange element space [11] for the new variable $\phi$ , with $\mathbb{P}_{l+1}(K)\subset S_{h}|_{K}$ for any $K\in\mathcal{T}_{h}$ .
•

$\boldsymbol{U}_{h}\subset\boldsymbol{H}_{0}({\rm curl\,})$ is the edge element space [18, 19] for the magnetic field $\boldsymbol{H}$ , with $(\mathbb{P}_{l}(K))^{d}\subset\boldsymbol{U}_{h}|_{K}$ and $(\mathbb{P}_{l}(K))^{2d-3}\subset{\rm curl\,}\boldsymbol{U}_{h}|_{K}$ for any $K\in\mathcal{T}_{h}$ .

In addition, we make the following assumptions for the above finite element spaces.

(A1)

There holds the inf-sup condition

(23)

\inf\limits_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{(\nabla\cdot\boldsymbol{v}_{h},q_{h})}{\|\boldsymbol{v}_{h}\|_{1}}\geq\beta_{0}\|q_{h}\|\qquad\forall~{}~{}q_{h}\in W_{h},

where $\beta_{0}>0$ is a constant independent of $h$ ;

(A2)

The diagram

(24)

\begin{CD}H_{0}^{1}@>{{\rm grad\,}}>{}>\boldsymbol{H}_{0}({\rm curl\,})\\ @V{}V{\pi_{h}^{s}}V@ VV\pi^{c}_{h}V\\ \ S_{h}@>{{\rm grad\,}}>{}>\boldsymbol{U}_{h}\end{CD}

is a commutative sequence. Here $\pi_{h}^{s}:~{}H_{0}^{1}(\Omega)\cap C^{0}(\overline{\Omega})\rightarrow S_{h}$ and $\pi_{h}^{c}:~{}\boldsymbol{H}_{0}({\rm curl\,})\rightarrow\boldsymbol{U}_{h}$ are the classical interpolation operators, and ${\rm grad\,}$ denotes the gradient operator. Note that the diagram (24) also indicates that

{\rm grad\,}S_{h}\subset\boldsymbol{U}_{h}.

We recall the following inverse inequality:

(25)

\|\nabla s_{h}\|_{\infty}\leq C_{inv}h^{-d/2}\|\nabla s_{h}\|\qquad\forall s_{h}\in S_{h},

where $C_{inv}>0$ is a constant independent of $h$ .

Remark 3.1.

There are many finite element spaces that satisfy (A1) and (A2). For example, we can choose $\boldsymbol{V}_{h}$ and $W_{h}$ as the Taylor-Hood element pairs. And the spaces $S_{h}$ and $\boldsymbol{U}_{h}$ can be respectively chosen as the lowest order Lagrange finite element space and the lowest order Nédélec edge element space [5, 6, 7].

Remark 3.2.

In this paper, we consider the conforming finite element spaces for the Navier-Stokes equation for simply of notations. In fact the nonconforming finite element spaces which satisfying the inf-sup condition (A1), with replace the global differential operators to element wise, are also work for the Navier-Stokes equations. For example the $CR$ - $\mathbb{P}_{0}$ finite element spaces are work for the Navier-Stokes equations.

3.2. Finite element scheme

In view of (17)-(20), we consider the following finite element scheme for the FHD model: Find $\phi_{h}\in S_{h}$ , $\boldsymbol{H}_{h}\in\boldsymbol{U}_{h}$ , $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ and $\tilde{p}_{h}\in W_{h}$ , such that

(26)	$\displaystyle a(\phi_{h};\phi_{h},\tau_{h})=-(g,\tau_{h})$	$\displaystyle\forall~{}~{}\tau_{h}\in S_{h},$
(27)	$\displaystyle(\boldsymbol{H}_{h},\boldsymbol{C}_{h})-(\nabla\phi_{h},\boldsymbol{C}_{h})=0$	$\displaystyle\forall~{}~{}\boldsymbol{C}_{h}\in\boldsymbol{U}_{h},$
(28)	$\displaystyle\eta(\nabla\boldsymbol{u}_{h},\nabla\boldsymbol{v}_{h})+b(\boldsymbol{u}_{h};\boldsymbol{u}_{h},\boldsymbol{v}_{h})-(\tilde{p}_{h},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f},\boldsymbol{v}_{h})$	$\displaystyle\forall~{}~{}\boldsymbol{v}_{h}\in\boldsymbol{V}_{h},$
(29)	$\displaystyle(\nabla\cdot\boldsymbol{u}_{h},q_{h})=0\$	$\displaystyle\forall~{}~{}q_{h}\in W_{h},$

where the trilinear form $b(\cdot;\cdot,\cdot):\boldsymbol{H}^{1}(\Omega)\times\boldsymbol{H}^{1}(\Omega)\times\boldsymbol{H}^{1}(\Omega)\rightarrow\mathbb{R}$ is defined by

b(\boldsymbol{w};\boldsymbol{u},\boldsymbol{v}):=\frac{\rho}{2}[((\boldsymbol{w}\cdot\nabla)\boldsymbol{u},\boldsymbol{v})-((\boldsymbol{w}\cdot\nabla)\boldsymbol{v},\boldsymbol{u})].

Remark 3.3.

Similar to (17)-(20), the finite element scheme (26) - (29) is a decoupled system. We can first solve the nonlinear equation (26) to get $\phi_{h}$ , and solve the Navier-Stokes system (28) - (29) to get $\boldsymbol{u}_{h}$ and $\tilde{p}_{h}$ . Then we can get $\boldsymbol{H}_{h}$ from (27). Finally, we can get $\boldsymbol{M}_{h}\in\boldsymbol{U}_{h}$ , the approximation of the magnetization $\boldsymbol{M}$ , from

(30)

(\boldsymbol{M}_{h},\boldsymbol{F}_{h})=\big{(}(\alpha(H_{h})-1)\boldsymbol{H}_{h},\boldsymbol{F}_{h}\big{)}\,\,\,\,\,\qquad\quad\forall~{}~{}\boldsymbol{F}_{h}\in\boldsymbol{U}_{h}

with $H_{h}:=|\boldsymbol{H}_{h}|$ , get $\psi_{h}\in W_{h}$ , the approximation of $\psi$ , from

(31)

(\psi_{h},\chi_{h})=\big{(}\beta(H_{h}),\chi_{h}\big{)}\ \ \qquad\qquad\qquad\forall~{}~{}\chi_{h}\in W_{h},

and get $p_{h}\in W_{h}$ , the approximation of the pressure $p$ , from

(32)

p_{h}=\tilde{p}_{h}+\mu_{0}\psi_{h}.

Remark 3.4.

It is easy to see that the trilinear form $b(\cdot;\cdot,\cdot)$ is skew-symmetric with respect to the last two variables, i.e.,

(33)

b(\boldsymbol{w};\boldsymbol{u},\boldsymbol{v})=-b(\boldsymbol{w};\boldsymbol{v},\boldsymbol{u})\qquad\forall~{}\boldsymbol{w},~{}\boldsymbol{u},~{}\boldsymbol{v}\in\boldsymbol{H}^{1}(\Omega),

and that

(34)

b(\boldsymbol{w};\boldsymbol{u},\boldsymbol{v})=\rho\big{(}(\boldsymbol{w}\cdot\nabla)\boldsymbol{u},\boldsymbol{v}\big{)}+\frac{\rho}{2}\big{(}(\nabla\cdot\boldsymbol{w})\boldsymbol{u},\boldsymbol{v}\big{)}\quad\forall~{}\boldsymbol{w},~{}\boldsymbol{u}\in\boldsymbol{H}^{1}(\Omega),\ \ \forall\ \boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega).

As a result, the following two relations hold:

(35)

b(\boldsymbol{w};\boldsymbol{v},\boldsymbol{v})=0\qquad\forall~{}\boldsymbol{w},\boldsymbol{v}\in\boldsymbol{H}^{1}(\Omega),\\

(36)

b(\boldsymbol{w};\boldsymbol{u},\boldsymbol{v})=\rho\big{(}(\boldsymbol{w}\cdot\nabla)\boldsymbol{u},\boldsymbol{v}\big{)}\quad\forall~{}\boldsymbol{w},~{}\boldsymbol{u}\in\boldsymbol{H}^{1}(\Omega)\text{ with }\operatorname{div}\boldsymbol{w}=0,\ \forall\ \boldsymbol{v}\in\boldsymbol{H}_{0}^{1}(\Omega).

Theorems 3.5-3.7 show the wellposedness of the finite element scheme (26) - (29).

Theorem 3.5.

The nonlinear discrete equation (26) has at least one solution $\phi_{h}\in S_{h}$ , and there holds

(37)

\|\nabla\phi_{h}\|\leq\frac{1}{\mu_{0}}\|\boldsymbol{H}_{e}\|.

Furthermore, if the external magnetic filed $\boldsymbol{H}_{e}$ satisfies

(38)

\|\boldsymbol{H}_{e}\|<\mu_{0}C_{inv}^{-1}C_{\alpha}^{-1}h^{d/2},

then (26) admits a unique solution $\phi_{h}\in S_{h}$ , and there holds

(39)

\|\nabla\phi_{h}\|_{\infty}<1/C_{\alpha}.

Proof.

Define an operator $A:~{}S_{h}\rightarrow S_{h}^{\prime}$ by

\langle A(\phi_{h}),\tau_{h}\rangle=a(\phi_{h};\phi_{h},\tau_{h})\qquad\forall~{}~{}\phi_{h},~{}\tau_{h}\in S_{h},

and define $\Phi:~{}S_{h}\rightarrow S_{h}^{\prime}$ by

\Phi(\phi_{h})=A(\phi_{h})+Q_{h}^{\prime}g\qquad\forall~{}~{}\phi_{h}\in S_{h},

where $Q_{h}^{\prime}:~{}(H_{0}^{1}(\Omega))^{\prime}\rightarrow S_{h}^{\prime}$ is given by

(40)

\langle Q_{h}^{\prime}g,\tau_{h}\rangle=\langle g,\tau_{h}\rangle=-\frac{1}{\mu_{0}}(\boldsymbol{H}_{e},\nabla\tau_{h})\qquad\forall~{}~{}\tau_{h}\in S_{h}.

Thus, (26) is equivalent to the equation

\Phi(\phi_{h})=0.

Lemma 2.4 implies

\langle A(\phi_{h}),\phi_{h}\rangle=a(\phi_{h};\phi_{h},\phi_{h})\geq\|\nabla\phi_{h}\|^{2}\qquad\forall~{}~{}\phi_{h}\in S_{h}.

The definition (40) and the Cauchy-Schwarz inequality imply

\langle Q_{h}^{\prime}g,\phi_{h}\rangle=\langle g,\phi_{h}\rangle\leq\frac{1}{\mu_{0}}\|\boldsymbol{H}_{e}\|\|\nabla\phi_{h}\|\qquad\forall~{}~{}\phi_{h}\in S_{h}.

Therefore, we have

\langle\Phi(\phi_{h}),\phi_{h}\rangle=\langle A(\phi_{h})+Q_{h}^{\prime}g,\phi_{h}\rangle\geq 0,\qquad\forall~{}~{}\phi_{h}\in S_{h}\text{ with }\|\nabla\phi_{h}\|=\frac{1}{\mu_{0}}\|\boldsymbol{H}_{e}\|.

For any given $0\neq\phi_{h0}\in S_{h}$ and $\epsilon>0$ , denote $\delta:=\min\left\{\frac{\epsilon}{2C_{1}},~{}\frac{\epsilon}{2C_{\alpha}\|\nabla\phi_{h0}\|_{\infty}}\right\}$ . Then for any $\phi_{h}\in S_{h}$ satisfying $\|\nabla(\phi_{h}-\phi_{h0})\|<\delta$ , by Lemmas 2.2 and 2.4 we have

		$\displaystyle\\|\Phi(\phi_{h})-\Phi(\phi_{h0})\\|_{S_{h}^{\prime}}=\\|A(\phi_{h})-A(\phi_{h0})\\|_{S_{h}^{\prime}}$
	$\displaystyle=$	$\displaystyle\sup\limits_{\tau_{h}\in S_{h}}\frac{\langle A(\phi_{h})-A(\phi_{h0}),\tau_{h}\rangle}{\\|\nabla\tau_{h}\\|}$
	$\displaystyle=$	$\displaystyle\sup\limits_{\tau_{h}\in S_{h}}\frac{a(\phi_{h};\phi_{h},\tau_{h})-a(\phi_{h0};\phi_{h0},\tau_{h})}{\\|\nabla\tau_{h}\\|}$
	$\displaystyle=$	$\displaystyle\sup\limits_{\tau_{h}\in S_{h}}\frac{a(\phi_{h};\phi_{h},\tau_{h})-a(\phi_{h};\phi_{h0},\tau_{h})+a(\phi_{h};\phi_{h0},\tau_{h})-a(\phi_{h0};\phi_{h0},\tau_{h})}{\\|\nabla\tau_{h}\\|}$
	$\displaystyle\leq$	$\displaystyle C_{1}\\|\nabla(\phi_{h}-\phi_{h0})\\|+C_{\alpha}\\|\nabla\phi_{h0}\\|_{\infty}\\|\nabla(\phi_{h}-\phi_{h0})\\|$
	$\displaystyle<$	$\displaystyle\epsilon.$

This means that $\Phi$ is continuous on the set $S_{h}\setminus\{0\}$ . Notice that Lemma 2.4 implies $\Phi$ is continuous at the point $0$ . Hence, $\Phi$ is continuous on $S_{h}$ .

On the other hand, the Riesz representation theorem implies that the spaces $S_{h}$ and $S_{h}^{\prime}$ are isometry. As a result, by [13, Page 279, Corollary 1.1] equation (26) admits at least one solution $\phi_{h}\in S_{h}$ satisfiying (37).

If $\boldsymbol{H}_{e}$ satisfies assumption (38), then from (37) and the inverse inequality (25) we get

\|\nabla\phi_{h}\|_{\infty}\leq C_{inv}h^{-d/2}\|\nabla\phi_{h}\|<1/C_{\alpha},

i.e. (39) holds. Thus, by following the same line as in the proof of the uniqueness of the weak solution $\phi$ in Theorem 2.5, we can easily obtain the uniqueness of the discrete solution $\phi_{h}$ . ∎

Theorem 3.6.

For any given $\phi_{h}\in S_{h}$ , equation (27) admits a unique solution $\boldsymbol{H}_{h}=\nabla\phi_{h}\in\boldsymbol{U}_{h}$ , which means ${\rm curl\,}\boldsymbol{H}_{h}=0.$

Proof.

For any given $\phi_{h}\in S_{h}$ , assumption (A2) on the finite element spaces implies that $\nabla\phi_{h}\in\boldsymbol{U}_{h}$ . Thus, taking $\boldsymbol{C}_{h}=\boldsymbol{H}_{h}-\nabla\phi_{h}\in\boldsymbol{U}_{h}$ in (27) yields $\boldsymbol{H}_{h}=\nabla\phi_{h}.$ ∎

The following well-posedness results for the finite element scheme (28) - (29) of the Navier-Stokes equations are standard (cf. [13]).

Theorem 3.7.

The finite element scheme (28) - (29) has at least one solution $(\boldsymbol{u}_{h},\tilde{p}_{h})\in\boldsymbol{V}_{h}\times W_{h}$ , and there holds

\|\nabla\boldsymbol{u}_{h}\|\leq\frac{1}{\eta}\|\boldsymbol{f}\|_{-1}.

Further more, if

(41)

(\tilde{\mathcal{N}}/\eta^{2})\|\boldsymbol{f}\|_{-1}<1\quad\text{with}\quad\tilde{\mathcal{N}}=\sup\limits_{\boldsymbol{w}_{h},\boldsymbol{u}_{h},\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{b(\boldsymbol{w}_{h};\boldsymbol{u}_{h},\boldsymbol{v}_{h})}{\|\nabla\boldsymbol{w}_{h}\|\|\nabla\boldsymbol{u}_{h}\|\|\nabla\boldsymbol{v}_{h}\|},

then the finite element solution $(\boldsymbol{u}_{h},\tilde{p}_{h})$ is unique.

3.3. Error analysis

In this subsection, we will give some error estimates for the finite element scheme (26) - (29), under some rational assumptions.

Firstly, we have the following error estimate for the discrete solution $\phi_{h}$ .

Theorem 3.8.

Let $\phi\in H_{0}^{1}(\Omega)\cap H^{l+2}(\Omega)$ be the solution of (17), and let $\phi_{h}\in S_{h}$ be the solution of (26). Assume that $\phi$ satisfies (22), i.e. $\|\nabla\phi\|_{L^{\infty}}<1/C_{\alpha},$ then there exists a positive constant $C_{s}$ , independent of $h$ , such that

\|\nabla(\phi-\phi_{h})\|\leq\frac{C_{s}(1+C_{1})}{1-C_{\alpha}\|\nabla\phi\|_{L^{\infty}}}h^{l+1}\|\phi\|_{l+2}.

Proof.

Let $\phi_{h}^{*}\in S_{h}$ be the classical interpolation of $\phi$ , satisfying the estimate

(42)

\|\nabla(\phi-\phi_{h}^{*})\|\leq C_{s}h^{l+1}\|\phi\|_{l+2},

where $C_{s}>0$ is a constant independent of $h$ and $\phi$ .

In view of the inequality

\|\nabla(\phi-\phi_{h})\|\leq\|\nabla(\phi-\phi_{h}^{*})\|+\|\nabla(\phi_{h}^{*}-\phi_{h})\|,

it suffices to estimate the term $\|\nabla(\phi_{h}^{*}-\phi_{h})\|$ . To this end, we apply Lemma 2.4 and the relation

a(\phi;\phi,\tau_{h})=a(\phi_{h};\phi_{h},\tau_{h}),\quad\forall\ \tau_{h}\in S_{h},

to get

	$\displaystyle\\|\nabla(\phi_{h}^{*}-\phi_{h})\\|^{2}$	$\displaystyle\leq$	$\displaystyle a(\phi_{h};\phi_{h}^{}-\phi_{h},\phi_{h}^{}-\phi_{h})$
		$\displaystyle=$	$\displaystyle a(\phi_{h};\phi_{h}^{}-\phi,\phi_{h}^{}-\phi_{h})+a(\phi_{h};\phi-\phi_{h},\phi_{h}^{*}-\phi_{h}).$

Since

|a(\phi_{h};\phi_{h}^{*}-\phi,\phi_{h}^{*}-\phi_{h})|\leq C_{1}\|\nabla(\phi-\phi_{h}^{*})\|\|\nabla(\phi_{h}^{*}-\phi_{h})\|

and

	$\displaystyle a(\phi_{h};\phi-\phi_{h},\phi_{h}^{*}-\phi_{h})=$	$\displaystyle a(\phi_{h};\phi,\phi_{h}^{}-\phi_{h})-a(\phi;\phi,\phi_{h}^{}-\phi_{h})$
	$\displaystyle=$	$\displaystyle\left(\big{(}\alpha(\|\nabla\phi_{h}\|)-\alpha(\|\nabla\phi\|)\big{)}\nabla\phi,\nabla(\phi_{h}^{*}-\phi_{h})\right)$
	$\displaystyle\leq$	$\displaystyle C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}\\|\nabla(\phi-\phi_{h})\\|\\|\nabla(\phi_{h}^{*}-\phi_{h})\\|$
	$\displaystyle\leq$	$\displaystyle C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}\big{(}\\|\nabla(\phi-\phi_{h}^{})\\|+\\|\nabla(\phi_{h}^{}-\phi_{h})\\|\big{)}\\|\nabla(\phi_{h}^{*}-\phi_{h})\\|,$

we have

\displaystyle\|\nabla(\phi_{h}^{*}-\phi_{h})\|\leq C_{1}\|\nabla(\phi-\phi_{h}^{*})\|+C_{\alpha}\|\nabla\phi\|_{L^{\infty}}\big{(}\|\nabla(\phi-\phi_{h}^{*})\|+\|\nabla(\phi_{h}^{*}-\phi_{h})\|\big{)},

which, together with (22), implies

\displaystyle\|\nabla(\phi_{h}^{*}-\phi_{h})\|\leq\frac{C_{1}+C_{\alpha}\|\nabla\phi\|_{L^{\infty}}}{1-C_{\alpha}\|\nabla\phi\|_{L^{\infty}}}\|\nabla(\phi-\phi_{h}^{*})\|.

This inequality, together with (42), indicates the desired conclusion. ∎

In light of Theorems 2.6 and 3.6, we know that

\boldsymbol{H}=\nabla\phi,\quad\boldsymbol{H}_{h}=\nabla\phi_{h}\quad\text{and}\quad{\rm curl\,}\boldsymbol{H}={\rm curl\,}\boldsymbol{H}_{h}=0.

Thus, by Theorem 3.8, we immediately get the following error estimate for the discrete solution $\boldsymbol{H}_{h}$ of (27).

Theorem 3.9.

Let $\boldsymbol{H}\in\boldsymbol{H}_{0}({\rm curl\,})$ be the solution of (18), and let $\boldsymbol{H}_{h}\in\boldsymbol{U}_{h}$ be the solution of (27). Then, under the same conditions as in Theorem 3.8, there holds

\|\boldsymbol{H}-\boldsymbol{H}_{h}\|+\|{\rm curl\,}(\boldsymbol{H}-\boldsymbol{H}_{h})\|\leq\frac{C_{s}(1+C_{1})}{1-C_{\alpha}\|\nabla\phi\|_{L^{\infty}}}h^{l+1}\|\phi\|_{l+2}.

The following error estimate for the Navier-Stokes system (28)-(29) is standard (cf. [13]).

Theorem 3.10.

Let $\boldsymbol{u}\in\boldsymbol{H}_{0}^{1}(\Omega)\cap\boldsymbol{H}^{l+2}(\Omega)$ and $\tilde{p}\in L_{0}^{2}(\Omega)\cap H^{l+1}(\Omega)$ be the solution of the Navier-Stokes system (19)-(20), and let $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ and $\tilde{p}_{h}\in W_{h}$ be the solution of (28)-(29). There exists a positive constant $C$ , independent of $h$ , such that

\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}+\|\tilde{p}-\tilde{p}_{h}\|\leq Ch^{l+1}(\|\boldsymbol{u}\|_{l+2}+\|\tilde{p}\|_{l+1}).

As shown in Remark 3.3, for the the magnetization $\boldsymbol{M}=(\alpha(H)-1)\boldsymbol{H}$ , the new variable $\psi=\beta(H)$ , and the pressure $p=\tilde{p}+\mu_{0}\psi$ , we can obtain their approximations $\boldsymbol{M}_{h},\psi_{h}$ and $p_{h}$ from (30), (31) and (32), respectively. In view of Theorems 3.9 and 3.10, we easily get the following error estimates for these three approximation solutions.

Corollary 3.11.

Assume that $\boldsymbol{M}\in\boldsymbol{H}^{l+1}(\Omega)$ and $\psi\in H^{l+1}(\Omega)$ . Under the same conditions as in Theorems 3.8 and 3.10, there exists positive constants $C_{M}$ and $C_{\psi}$ , independent of $h$ , such that

(43)		$\displaystyle\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|\leq C_{M}h^{l+1}\left(\\|\boldsymbol{M}\\|_{l+1}+\frac{(1+C_{1})(1+C_{1}+C_{\alpha}\\|\nabla\phi\\|_{\infty})}{1-C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}}\\|\phi\\|_{l+2}\right),$
(44)		$\displaystyle\\|\psi-\psi_{h}\\|\leq\left(C_{\psi}\\|\psi\\|_{l+1}+\frac{M_{s}C_{s}(1+C_{1})}{1-C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}}\\|\phi\\|_{l+2}\right)h^{l+1},$
(45)		$\displaystyle\\|p-p_{h}\\|\leq\left(C\\|\boldsymbol{u}\\|_{l+2}+C\\|\tilde{p}\\|_{l+1}+\mu_{0}C_{\psi}\\|\psi\\|_{l+1}+\frac{\mu_{0}M_{s}C_{s}(1+C_{1})}{1-C_{\alpha}\\|\nabla\phi\\|_{L^{\infty}}}\\|\phi\\|_{l+2}\right)h^{l+1}.$

Proof.

Equation (30) implies that

\boldsymbol{M}_{h}=Q_{h}^{c}((\alpha(H_{h})-1)\boldsymbol{H}_{h}),

where $Q_{h}^{c}:~{}\boldsymbol{L}^{2}(\Omega)\rightarrow\boldsymbol{U}_{h}$ is the $L^{2}$ orthogonal projection operator. Using the Langevin law (4), we have

\displaystyle\boldsymbol{M}-\boldsymbol{M}_{h}=\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}+Q_{h}^{c}\big{(}(\alpha(H)-1)\boldsymbol{H}-(\alpha(H_{h})-1)\boldsymbol{H}_{h}\big{)}.

Thus, by Lemma 2.2 (i), the boundedness of projection $Q_{h}^{c}$ , and the relation $\boldsymbol{H}=\nabla\phi$ , we get

	$\displaystyle\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|$	$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+\\|(\alpha(H)-\alpha(H_{h}))\boldsymbol{H}\\|+\\|(\alpha(H_{h})-1)(\boldsymbol{H}-\boldsymbol{H}_{h})\\|$
		$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+(C_{\alpha}\\|\boldsymbol{H}\\|_{\infty}+C_{1}+1)\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|$
		$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+(C_{\alpha}\\|\nabla\phi\\|_{\infty}+C_{1}+1)\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|,$

which, together with Theorem 3.9 and the standard error estimation of the projection, gives the desired estimate for $\boldsymbol{M}_{h}$ .

Similarly, equation (31) implies that

\boldsymbol{\psi}_{h}=Q_{h}^{w}(\beta(H_{h})),

where $Q_{h}^{w}:~{}L^{2}(\Omega)\rightarrow W_{h}$ is the $L^{2}$ orthogonal projection. Then by Lemma 2.2 (ii) we obtain

	$\displaystyle\\|\psi-\psi_{h}\\|$	$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+\\|Q_{h}^{w}\big{(}\beta(H)-\beta(H_{H})\big{)}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+\\|\big{(}\beta(H)-\beta(H_{h})\big{)}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+M_{s}\\|H-H_{h}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+M_{s}\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|,$

which, together with Theorem 3.9 and the projection property, yields the desired estimate (44).

Finally, (32) means that

\|p-p_{h}\|\leq\|\tilde{p}-\tilde{p}_{h}\|+\mu_{0}\|\psi-\psi_{h}\|,

which, together with Theorem 3.10 and estimate (44), indicates the desired result (45). This completes the proof. ∎

4. Numerical experiments

This section is devoted to three numerical examples to verify the performance of the mixed finite element methods. In all the examples, we solve the nonlinear system (26) - (32) by using the iFEM package [10] and Algorithm 4.1.

Algorithm 4.1.

Given $\phi_{h}^{0}\in S_{h}$ and $\boldsymbol{u}_{h}^{0}\in\boldsymbol{V}_{h}$ , find $\phi_{h}\in S_{h}$ , $\boldsymbol{H}_{h}\in\boldsymbol{U}_{h}$ , $\boldsymbol{M}_{h}\in\boldsymbol{U}_{h}$ , $\boldsymbol{u}_{h}\in\boldsymbol{V}_{h}$ , $\tilde{p}_{h}\in W_{h}$ , $\psi_{h}\in W_{h}$ , and $p_{h}\in W_{h}$ through five steps:

Step 1.

For $n=1,2,\dots,L$ do

(1)

Solving the nonlinear system (26) as: Find $\phi_{h}^{n}\in S_{h}$ such that

$a(\phi_{h}^{n-1};\phi_{h}^{n},\tau_{h})=-(g,\tau_{h})\qquad\forall~{}\tau_{h}\in S_{h}.$

(2)

Solving the nonlinear saddle point system (28)-(29) as: Find $\boldsymbol{u}_{h}^{n}\in\boldsymbol{V}_{h}$ and $\tilde{p}_{h}^{n}\in W_{h}$ such that

	$\displaystyle\eta(\nabla\boldsymbol{u}_{h}^{n},\nabla\boldsymbol{v}_{h})+b(\boldsymbol{u}_{h}^{n-1};\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{h})-(\tilde{p}_{h}^{n},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f},\boldsymbol{v}_{h})$	$\displaystyle\forall~{}\boldsymbol{v}_{h}\in\boldsymbol{V}_{h},$
	$\displaystyle(\nabla\cdot\boldsymbol{u}_{h}^{n},q_{h})=0$	$\displaystyle\forall~{}q_{h}\in W_{h}.$

Step 2.

Let $\phi_{h}=\phi_{h}^{L}$ , $\boldsymbol{u}_{h}=\boldsymbol{u}_{h}^{L}$ and $\tilde{p}_{h}=\tilde{p}_{h}^{L}$ .
Step 3.

Solving equation (27) to get $\boldsymbol{H}_{h}\in\boldsymbol{U}_{h}$ .
Step 4.

Solving equations (30) and (31) to get $\boldsymbol{M}_{h}\in\boldsymbol{U}_{h}$ and $\psi_{h}\in W_{h}$ , respectively.
Step 5.

Substituting $\tilde{p}_{h}$ and $\psi_{h}$ into (32) to get $p_{h}\in W_{h}$ .

Remark 4.2.

From the convergence theory of Newton-type methods [15, 28], we know that the iterative solution of Algorithm 4.1 will converge to the exact solution, provided that the iteration number $L$ is big enough and the initial guess is nearby the exact solution. In fact, in all the subsequent numerical examples we choose the initial guess $\phi_{h}^{0}$ as the corresponding finite element solution of the Poisson equation

\Delta\tilde{\phi}=g\qquad\text{in }\ \Omega

with the same boundary condition as that of the exact solution $\phi$ , and choose the initial guess $\boldsymbol{u}_{h}^{0}$ as the corresponding finite element solution of the Stokes equations

\left\{\begin{array}[]{ll}-\eta\Delta\tilde{\boldsymbol{u}}+\nabla p=\boldsymbol{f}&\text{in }\Omega,\\ \nabla\cdot\tilde{\boldsymbol{u}}=0&\text{in }\Omega\end{array}\right.

with the same boundary condition as that of the exact solution $\boldsymbol{u}$ . In so doing, the choice $L=2$ works well in the algorithm.

Refer to caption — Figure 1. The domain $\Omega=(0,1)^{2}$ : $4\times 4$ (left) and $8\times 8$ (right).

Example 4.3.

This is a 2D test. We take $\Omega=(0,1)^{2}$ , and use $N\times N$ uniform triangular mesh (cf. Figure 1) with $N=4,~{}8,~{}16,~{}32,~{}64,~{}128$ . The exact solution of the FHD model (16) is given by

	$\displaystyle\phi(x,y)=(x^{2}-x)(y^{2}-y),$
	$\displaystyle\boldsymbol{H}(x,y)=((2x-1)(y^{2}-y),(x^{2}-x)(2y-1))^{\intercal},$
	$\displaystyle\boldsymbol{M}(x,y)=M_{s}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{\boldsymbol{H}}{H},$
	$\displaystyle\boldsymbol{u}(x,y)=(\sin(\pi y),\sin(\pi x))^{\intercal},$
	$\displaystyle p(x,y)=60x^{2}y-20y^{3}-5,$

with the parameters $M_{s}$ , $\rho$ , $\eta$ , $\mu_{0}$ and $\gamma$ all being chosen as $1$ .

We use the conforming linear ( $\mathbb{P}_{1}$ ) element to discretize the variable $\phi$ , the lowest order edge element $NE_{0}$ to discretize the variables $\boldsymbol{H}$ and $\boldsymbol{M}$ , the $CR$ (nonconforming- $\mathbb{P}_{1}$ ) element to discretize the variables $\boldsymbol{u}$ , and the piecewise constant ( $\mathbb{P}_{0}$ ) element to discretize the variables $\tilde{p}$ , $\psi$ and $p$ . Note that such a combination of finite element spaces corresponds to $l=0$ , then we easily see from Theorems 3.8 - 3.10 and Corollary 3.11 that the theoretical accuracy of the scheme is $\mathcal{O}(h)$ . Numerical results are listed in Table 1.

Table 1. Numerical results for Example 4.3.

$N$	$\frac{\\|\nabla(\phi-\phi_{h})\\|}{\\|\nabla\phi\\|}$	$\frac{\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|_{c}}{\\|\boldsymbol{H}\\|_{c}}$	$\frac{\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|}{\\|\boldsymbol{M}\\|}$	$\frac{\\|\boldsymbol{u}-\boldsymbol{u}_{h}\\|_{1,h}}{\\|\boldsymbol{u}\\|_{1}}$	$\frac{\\|p-p_{h}\\|}{\\|p\\|}$	$\\|{\rm curl\,}\boldsymbol{H}_{h}\\|_{\infty}$
$4$	$0.3943$	$0.3943$	$0.3941$	$0.7424$	$0.3083$	$0.0003e-12$
$8$	$0.2023$	$0.2023$	$0.2023$	$0.4385$	$0.1524$	$0.0013e-12$
$16$	$0.1018$	$0.1018$	$0.1018$	$0.2352$	$0.0717$	$0.0083e-12$
$32$	$0.0510$	$0.0510$	$0.0510$	$0.1207$	$0.0342$	$0.0315e-12$
$64$	$0.0255$	$0.0255$	$0.0255$	$0.0609$	$0.0167$	$0.1485e-12$
$128$	$0.0128$	$0.0128$	$0.0128$	$0.0305$	$0.0083$	$0.6481e-12$
order	$0.9900$	$0.9900$	$0.9899$	$0.9207$	$1.0439$	–

The other two experiments, Examples 4.4 and 4.5, are 3D tests, with $\Omega=(0,1)^{3}$ and use $N\times N\times N$ uniform tetrahedral meshes (cf. Figure 2) with $N=4,~{}8,~{}16$ .

Example 4.4 (The lowest order approximation).

In this 3D test, the exact solution of the FHD model (16) is given by

	$\displaystyle\phi=\sin(\pi x)\sin(\pi y)\sin(\pi z),$
	$\displaystyle\boldsymbol{H}=\pi\begin{pmatrix}\cos(\pi x)\sin(\pi y)\sin(\pi z)\\ \sin(\pi x)\cos(\pi y)\sin(\pi z)\\ \sin(\pi x)\sin(\pi y)\cos(\pi z)\end{pmatrix},$
	$\displaystyle\boldsymbol{M}=M_{s}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{\boldsymbol{H}}{H},$
	$\displaystyle\boldsymbol{u}=\begin{pmatrix}\sin(\pi y),&\sin(\pi z),&\sin(\pi x)\end{pmatrix}^{\intercal},$
	$\displaystyle p=120x^{2}yz-40y^{3}z-40yz^{3},$

with the parameters $M_{s}$ , $\gamma$ , $\rho$ and $\eta$ all being choosen as $1$ and $\mu_{0}=10$ .

We use the the conforming $\mathbb{P}_{1}$ -element to discretize the variable $\phi$ , the lowest order Nédélec finite element space $NE_{0}$ [18, 19] to discretize the variables $\boldsymbol{H}$ and $\boldsymbol{M}$ , the nonconforming $CR$ element [12] to discretize the variable $\boldsymbol{u}$ , and the piecewise constant $\mathbb{P}_{0}$ element to discretize the variables $\tilde{p}$ , $\psi$ and $p$ . With these settings, from Theorems 3.8 - 3.10 and Corollary 3.11 we see that the theoretical accuracy of the scheme is $\mathcal{O}(h)$ . Numerical results are given in Table 2.

Table 2. Numerical results for Example 4.4.

$N$	$\frac{\\|\nabla(\phi-\phi_{h})\\|}{\\|\nabla\phi\\|}$	$\frac{\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|_{c}}{\\|\boldsymbol{H}\\|_{c}}$	$\frac{\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|}{\\|\boldsymbol{M}\\|}$	$\frac{\\|\boldsymbol{u}-\boldsymbol{u}_{h}\\|_{1,h}}{\\|\boldsymbol{u}\\|_{1}}$	$\frac{\\|p-p_{h}\\|}{\\|p\\|}$	$\\|{\rm curl\,}\boldsymbol{H}_{h}\\|_{\infty}$
$4$	$0.4739$	$0.4739$	$0.4755$	$0.2066$	$0.3665$	$0.0213e-12$
$8$	$0.2491$	$0.2491$	$0.2506$	$0.1082$	$0.1803$	$0.1048e-12$
$16$	$0.1262$	$0.1262$	$0.1271$	$0.0551$	$0.0835$	$0.6821e-12$
order	$0.9545$	$0.9545$	$0.9515$	$0.9529$	$1.0671$	$--$

Example 4.5 (A higher order approximation).

In this 3D test, the exact solution of the FHD model (16) is given by

	$\displaystyle\phi(x,y,z)=(x^{2}+y^{2}+z^{2})(x^{2}-x)(y^{2}-y)(z^{2}-z),$
	$\displaystyle\boldsymbol{H}(x,y,z)=\begin{pmatrix}2x(x^{2}-x)(y^{2}-y)(z^{2}-z)+(2x-1)(x^{2}+y^{2}+z^{2})(y^{2}-y)(z^{2}-z)\\ 2y(x^{2}-x)(y^{2}-y)(z^{2}-z)+(2y-1)(x^{2}+y^{2}+z^{2})(x^{2}-x)(z^{2}-z)\\ 2z(x^{2}-x)(y^{2}-y)(z^{2}-z)+(2z-1)(x^{2}+y^{2}+z^{2})(x^{2}-x)(y^{2}-y)\end{pmatrix},$
	$\displaystyle\boldsymbol{M}(x,y,z)=M_{s}\left(\coth(\gamma H)-\frac{1}{\gamma H}\right)\frac{\boldsymbol{H}}{H},$
	$\displaystyle\boldsymbol{u}(x,y,z)=\begin{pmatrix}\sin(\pi y),&\sin(\pi z),&\sin(\pi x)\end{pmatrix}^{\intercal},$
	$\displaystyle p(x,y,z)=120x^{2}yz-40y^{3}z-40yz^{3},$

where the parameters $M_{s}$ , $\gamma$ , $\rho$ , $\eta$ and $\mu_{0}$ are all choosen as $1$ .

We use the conforming quadratic ( $\mathbb{P}_{2}$ ) element to discretize the variable $\phi$ , the first order Nédélec edge element $NE_{1}$ [18, 19] to discretize the variables $\boldsymbol{H}$ and $\boldsymbol{M}$ , the Taylor-Hood $\mathbb{P}_{2}$ - $\mathbb{P}_{1}$ element to discretize the variables $\boldsymbol{u}$ and $\tilde{p}$ , and the continuous linear ( $\mathbb{P}_{1}$ ) element to discretize the variables $\psi$ and $p$ . With these settings, we see that the theoretical accuracy of the scheme is $\mathcal{O}(h^{2})$ . We list numerical results in Table 3.

Table 3. Numerical results for Example 4.5.

$N$	$\frac{\\|\nabla(\phi-\phi_{h})\\|}{\\|\nabla\phi\\|}$	$\frac{\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|_{c}}{\\|\boldsymbol{H}\\|_{c}}$	$\frac{\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|}{\\|\boldsymbol{M}\\|}$	$\frac{\\|\nabla(\boldsymbol{u}-\boldsymbol{u}_{h})\\|}{\\|\nabla\boldsymbol{u}\\|}$	$\frac{\\|p-p_{h}\\|}{\\|p\\|}$	$\\|{\rm curl\,}\boldsymbol{H}_{h}\\|$
$4$	$0.1369$	$0.1369$	$0.1369$	$0.0450$	$0.0535$	$0.2970e-08$
$8$	$0.0372$	$0.0372$	$0.0372$	$0.0081$	$0.0135$	$0.7640e-08$
$16$	$0.0095$	$0.0095$	$0.0095$	$0.0016$	$0.0034$	$1.6060e-08$
order	$1.9245$	$1.9245$	$1.9245$	$2.4069$	$1.9880$	$--$

From Examples 4.3, 4.4 and 4.5, we have the following observations:

•

From Tables 1 and 2, we see that the errors in $H^{1}$ semi-norm for $\phi$ , $\boldsymbol{H}({\rm curl\,})$ norm for $\boldsymbol{H}$ , $L^{2}$ norm for $\boldsymbol{M}$ and $p$ , and discrete $H^{1}$ norm for $\boldsymbol{u}$ , all have the first (optimal) order rate.
•

From Table 4.5, we see that the errors in $H^{1}$ semi-norm for $\phi$ and $\boldsymbol{u}$ , $\boldsymbol{H}({\rm curl\,})$ norm for $\boldsymbol{H}$ , $L^{2}$ norm for $\boldsymbol{M}$ and $p$ , all have the second (optimal) order rate.
•

From Tables 1 - 3, we can see that the numerical scheme preserves ${\rm curl\,}\boldsymbol{H}_{h}=0$ exactly.

References

[1] Y. Amirat and K. Hamdache. Global weak solutions to a ferrofluid flow model. Mathematical Methods in the Applied Sciences, 31(2):123–151, 2008.
[2] Y. Amirat and K. Hamdache. Strong solutions to the equations of a ferrofluid flow model. Journal of Mathematical Analysis and Applications, 353(1):271–294, 2009.
[3] Y. Amirat and K. Hamdache. Unique solvability of equations of motion for ferrofluids. Nonlinear Analysis Theory Methods and Applications, 73(2):471–494, 2010.
[4] Y. Amirat, K. Hamdache, and F. Murat. Global weak solutions to equations of motion for magnetic fluids. Journal of Mathematical Fluid Mechanics, 10(3):326–351, 2008.
[5] D. N. Arnold, R. S. Falk, and J. Gopalakrishnan. Mixed finite element approximation of the vector laplacian with dirichlet boundary conditions. Mathematical Models and Methods in Applied Sciences, 22(09):1250024, 2012.
[6] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homological techniques, and applications. Acta Numerica, 15:1–155, 2006.
[7] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American Mathematical Society, 47(2):281–354, 2010.
[8] L. Boccardo, F. Murat, and J.-P. Puel. $L^{\infty}$ estimate for some nonlinear elliptic partial differential equations and application to an existence result. SIAM J. Math. Anal., 23(2): 326–333, 1992.
[9] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods. Springer-Verlag, New York, 1991.
[10] L. Chen. ifem: an integrated finite element methods package in matlab. Technical report, University of California at Irvine, 2009.
[11] P. Ciarlet. The finite element method for elliptic problems. Amsterdam: North-Holland, 1978.
[12] M. Crouzeix and P. A. Raviart. Conforming and nonconforming finite element methods for solving the stationary stokes equations i. Revue française d automatique informatique recherche opérationnelle Mathématique, 7(3), 1973.
[13] V. Girault and P. A. Raviart. Finite element methods for Navier-Stokes equations. Springer-Verlag, New York, 1986.
[14] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations: theory and algorithms, volume 5. Springer Science & Business Media, 2012.
[15] A. Gil, J. Segura, and N. M. Temme, Numerical methods for special functions, Society for Industrial and Applied Mathematics, 2007.
[16] P. Knobloch, S. S, and L. Tobiska. Numerical treatment of a free surface problem in ferrohydrodynamics. Pamm, 10(1):573–574, 2010.
[17] O. Lavrova, G. Matthies, T. Mitkova, V. Polevikov, and L. Tobiska. Numerical treatment of free surface problems in ferrohydrodynamics. Journal of Physics Condensed Matter, 18(38):2657–2669, 2006.
[18] J.-C. Nédélec. Mixed finite elements in $\mathbb{R}^{3}$ . Numerische Mathematik, 35(3):315–341, 1980.
[19] J.-C. Nédélec. A new family of mixed finite elements in $r^{3}$ . Numerische Mathematik, 50(1):57–81, 1986.
[20] R. Nochetto, A. Salgado, and I. Tomas. The equations of ferrohydrodynamics: modeling and numerical methods. Mathematical Models and Methods in Applied Sciences, 26(13), 2015.
[21] R. Nochetto, A. Salgado, and I. Tomas. On the Dynamics of Ferrofluids: Global Weak Solutions to the Rosensweig System and Rigorous Convergence to Equilibrium. SIAM Journal on Mathematical Analysis, 51(6), 4245-4286, 2019.
[22] Q. Pankhurst, J. Connolly, S. Jones, and J. Dobson. Applications of magnetic nanoparticles in biomedicine. Journal of Physics D-Applied Physics, 36: R167-R181, 2003.
[23] R. Rosensweig. Magnetic fluids. Ann. Rev. Fluid Mech., 19: 437-463, 1987.
[24] R. Rosensweig. Ferrohydrodynamics. Cambridge University Press, Cambridge, UK, 1985.
[25] M. I. Shliomis. Effective viscosity of magnetic suspensions. Sov. Phys. jetp, 34(6):1291–1294, 1972.
[26] M. I. Shliomis. Ferrofluids: Magnetically controllable fluids and their applications. Lecture Notes in Physics, 2002.
[27] S. M. Snyder, T. Cader, and B. A. Finlayson. Finite element model of magnetoconvection of a ferrofluid. Journal of Magnetism and Magnetic Materials, 262(2):269–279, 2003.
[28] E. Sǔli and D. Mayers, An introduction to Numerical Analysis, Cabridge University Press, 2003.
[29] Y. Wu and X. Xie, Energy-preserving mixed finite element methods for a ferrofluid flow model, submitted; arXiv: 2206.03129, 2022.
[30] G. Yoshikawa, K. Hirata, F. Miyasaka, and O. Yu. Numerical analysis of transitional behavior of ferrofluid employing mps method and fem. IEEE Transaction on Magnetics, 47(5):1370–1373, 2010.
[31] M. Zahn. Magnetic fluid and nanoparticle applications to nanotechnology. Journal of Nanoparticle Research, 3(73):73–78, 2001.
[32] G. Zhang, X. He, and X. Yang. Decoupled, Linear, and Unconditionally energy stable fully discrete finite element numerical scheme for a two-phase ferrohydrodynamics model. SIAM Journal on Scientific Computing, 43(1): B167-B193, 2021.

	$\displaystyle\\|\nabla(\phi_{h}^{*}-\phi_{h})\\|^{2}$	$\displaystyle\leq$	$\displaystyle a(\phi_{h};\phi_{h}^{}-\phi_{h},\phi_{h}^{}-\phi_{h})$
		$\displaystyle=$	$\displaystyle a(\phi_{h};\phi_{h}^{}-\phi,\phi_{h}^{}-\phi_{h})+a(\phi_{h};\phi-\phi_{h},\phi_{h}^{*}-\phi_{h}).$

	$\displaystyle\\|\boldsymbol{M}-\boldsymbol{M}_{h}\\|$	$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+\\|(\alpha(H)-\alpha(H_{h}))\boldsymbol{H}\\|+\\|(\alpha(H_{h})-1)(\boldsymbol{H}-\boldsymbol{H}_{h})\\|$
		$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+(C_{\alpha}\\|\boldsymbol{H}\\|_{\infty}+C_{1}+1)\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|$
		$\displaystyle\leq\\|\boldsymbol{M}-Q_{h}^{c}\boldsymbol{M}\\|+(C_{\alpha}\\|\nabla\phi\\|_{\infty}+C_{1}+1)\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|,$

	$\displaystyle\\|\psi-\psi_{h}\\|$	$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+\\|Q_{h}^{w}\big{(}\beta(H)-\beta(H_{H})\big{)}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+\\|\big{(}\beta(H)-\beta(H_{h})\big{)}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+M_{s}\\|H-H_{h}\\|$
		$\displaystyle\leq\\|\psi-Q_{h}^{w}\psi\\|+M_{s}\\|\boldsymbol{H}-\boldsymbol{H}_{h}\\|,$