1 Introduction

The flow problem of viscous, incompressible, electrically conducting fluids in the channels and ducts under a uniform oblique magnetic field is of great interest, because it has many practical applications in the field of magnetohydrodynamics (MHD), such as the blood flow control and measurements, MHD flowmeters, MHD power generation, and accelerators. Only for some very special cases, the problem can be exactly solved. Therefore, for the sake of application, an effective numerical method to solve the MHD flow problem becomes an very important research topic.

It is well-known that, for large values of the Hartmann number, the MHD flow problem is convection-dominated such that the exact solution may display localized phenomena, such as interior or boundary layers. For solving the MHD flow problem, there have been many research works using various numerical methods, such as the finite element, finite difference, and boundary element methods. The readers car refer to Refs. [1]-[2] and many references cited therein. Unfortunately, because of the lack of stability or accuracy, most of conventional numerical methods cannot solve the layer problems efficiently. In Ref. [1], it was pointed out that the common deficiency among the existing numerical methods is that they produce accurate results only in several special configurations for the MHD duct flow problems, but the Hartmann number Ha cannot exceed about 10². However, in many industrial applications, the important range of the Hartmann number is 10² < Ha < 10⁶. For this reason, some seemingly robust numerical methods with large Hartmann numbers were suggested, such as the stabilized finite element method (FEM)^[1] and least-squares FEM^[2] using residual-free bubble functions. Even so, for large Hartmann numbers, some spurious oscillations still appear near the boundary layer regions when these two bubble-stabilized methods are used (see Refs. [1] and [2] or also the discussion on this in Ref. [3]). Later, a tailored finite point method was proposed^[3], and high accuracy was achieved even for large Hartmann numbers. Beyond that, the boundary elements method^[4] and the two-level element free Galerkin method^[5] were also presented for solving the MHD duct flow problem at high Hartmann numbers, respectively. The theoretical error analysis of these methods above remains seriously poor except for the least-squares FEM^[2]. In addition, those methods were only investigated on isotropic meshes where the aspect ratio of mesh elements was uniformly bounded.

For such layer problems, a special anisotropic mesh adaptivity is more desirable, i.e., the mesh elements should adapt in both size and shape in order to approximate the solution better. There have been a well-developed literature on a posteriori error estimation but on isotropic meshes (see Refs. [6] and [7] and the references therein). Some types of the posteriori error estimate have already been extended to anisotropic meshes with large aspect ratio by Kunert^[8-10], but no anisotropic mesh refinements were made. Due to large aspect ratio, there appears the alignment measure in a posteriori error estimation to quantify the alignment between the solution and the mesh, which is one of the most important contributions of Kunert.

Based on a priori error estimate, the concept of metrics was mostly used in anisotropic mesh adaptation^[11-12]. A priori error analysis on anisotropic meshes can be found in Refs. [13] and [14] and the references therein. Moreover, a posteriori error indicator was used in anisotropic mesh adaptation^[15]. The concept of metrics was combined with the solution of a dual problem and a posteriori error analysis in order to generate an adaptive mesh^[16]. The metric field was directly obtained from a posteriori error estimate and used in the anisotropic mesh adaptation^[17].

In this paper, we present an anisotropic adaptive FEM to solve the MHD duct flow problem in a straight channel of uniform cross-section. By the analytic tools of Kunert^[9], a residual error estimator is proposed and provides two-sided bounds on the error on anisotropic meshes. However, the error estimator is not available for anisotropic adaptive refinement, because it cannot represent the contributions to the whole error in different directions, and the alignment measure, in which the exact solution is contained, enters the upper error bound. In order to use the alignment measure as a guide to design a proper anisotropic mesh, we reformulate its definition and use the well-known Zienkiewicz-Zhu (Z-Z) estimates^[6]. Then, a computable anisotropic error indicator is derived, and the corresponding locally error indicators are used to adjust the mesh sizes in different directions. Note that an analogous error indicator was proposed for anisotropic adaptive refinement for elliptic and parabolic problems in Ref. [15]. Finally, we implement an adaptive algorithm using the software FREEFEM++^[18] with the mesh generator Bamg^[19] based on a right metric field obtained by the locally error estimates and Z-Z estimates of the error gradient matrix. By this method, the layer information from some directions can be well captured such that the number of mesh vertices is dramatically reduced for a given level of accuracy. Numerical results confirm the error analysis.

The paper is organized as follows. We first specify our notation and give some preliminary results in Section 2. Section 3 is devoted to analytical tools, such as inverse inequalities, alignment measure, and anisotropic interpolation estimates. An error estimator is proposed for the MHD problem, and the two-sided bounds on the error are derived on anisotropic meshes in Section 4. In Section 5, we derive an anisotropic error indicator and implement an adaptive algorithm. Therein, a series of numerical simulations for three test problems at different Hartmann numbers are carried out to demonstrate effectiveness of the anisotropic adaptive FEM.

2 Preliminaries

In this section, we give some notation and the anisotropic mesh description, present some preliminary results, and finally introduce the model problem and its conforming discretization.

2.1 Notation

For a bounded domain Ω with the Lipschitz continuous boundary ∂Ω, let S be any given open subset of Ω. For the spaces L²(S) and L²(S)², the integral inner product and the norm are denoted by (·, ·)_S and ||·||_S, respectively. If S = Ω, the subscript will be omitted. Let |S| be the Lebesgue measure of S, and in particular, |s| is the length of a segment s. P_k(S) is the space of polynomials of order k or less, where k is a given non-negative integer. For simplicity, instead of x ≤ cy and c₁x ≤ y ≤ c₂x, we use the abbreviated notation x y and x ~ y, where the constants c, c₁, and c₂ are independent of x, y, the Hartmann number, and the mesh.

By F = {T_h}, we denote a family of triangulations T_h of Ω, where any two triangles are either disjoint or share a common vertex or edge. Let E_h be the set of all edges of T_h, let E_h^int be the set of interior, and let E_h^ext be the set of boundary. We denote by N_h (N_h^int, N_h^ext) the set of all (interior, boundary) vertices of T_h. Define for a∈N_h, σ∈E_h, and K∈T_h, T_a := {L∈T_h; a∈L}, T_σ := {L∈T_h; σ ⊂ L}, and T_K := {L∈T_h; L ∩ K ≠ ∅}, respectively. Furthermore, three auxiliary subdomains need to be defined by

(1)

Next, the notation n_K always denotes the exterior unit normal vector for any given K∈T_h, and n_σ denotes the unit normal vector for any given edge σ∈E_h, of which the orientation is arbitrarily chosen but coinciding with the exterior normal of the domain Ω for boundary edges and fixed for interior edges.

For a function ϕ and an edge σ∈E_h^int shared by two triangles K and L, where n_σ points from K to L, we define the jump operator through σ by

(2)

For any σ∈E_h^ext, set .

2.2 Anisotropic meshes

For an arbitrary (anisotropic) triangle K∈T_h, we enumerate its vertices such that P₀P₁ is the longest edge and |P₁P₂| ≥ |P₀P₂|. Moreover, we define two orthogonal vectors p_i with the length h_{i, K} := |p_i| (i = 1, 2), as shown in Fig. 1. Notice that h_{1, K} ≥ h_{2, K}, set h_{min, K} := h_{2, K}, and h_{max, K} := h_{1, K}. For σ ⊂ ∂K, let

(3)

Define two 2 × 2 matrices A_K and C_K by

(4)

Additionally, we require that the mesh T_h satisfies the following requirements:

(i) The number of triangles contained in T_a is bounded uniformly for each fixed vertex a∈N_h.

(ii) The dimension of adjacent triangles is close to each other, i.e.,

(5)

For convenience, for the common edge σ shared by two triangles K and L, set

(6)

An advantage of doing this is that they are no longer related to either of K or L but only to σ. For the boundary edges, the definitions are changed in the obvious way. Obviously, they satisfy h_σ ~ h_{σ, K} ~ h_{σ, L} and h_{min, σ} ~ h_{min, K} ~ h_{min, L}. For more details on anisotropic meshes, see Refs. [8] and [9].

2.3 Definition of problem and its conforming discretization

In a straight channel of a uniform cross-section Ω, we assume that there exists a laminar fully developed flow of a viscous, incompressible, and electrically conducting fluid. Let u be the velocity, and let B be the induced magnetic field. Here, the direction of the uniform transverse applied magnetic field B₀ may be arbitrary to the x-axis, and the fields u and B are parallel to the z-axis. Further, The fluid is driven down by a constant pressure gradient. The governing equations for the above duct flow in a dimensionless form with suitable boundary conditions can be posed as follows^{[4, 20]}:

(1)

where Ha = B₀l(δ/ν)^1/2 is the Hartmann number, l is the characteristic length of the duct, B₀ is the intensity of the external magnetic field, and ν and δ are the viscosity coefficient and electric conductivity of the fluid, respectively. a = (cos α, sin α)^T, where α is the angle between the x-axis and the externally applied magnetic field B₀. ∂Ω = Γ_D ∪ Γ_N, where Γ_D ∩ Γ_N = ∅ and Γ_D has a positive measure. We call Γ_N the conducting part and Γ_D the insulated part of the boundary ∂Ω.

Let W_0B(Ω) := H₀¹(Ω) × H_B¹ (Ω), where two particular subspaces of H¹(Ω) are

(7)

Due to the Poincaré inequality, the space W_0B(Ω) will be equipped with the product norm,

(8)

The weak form of the problem (1) is to find (u, B)∈W_0B(Ω) such that

(2)

In order to define the finite element approximation, let V^h(Ω) be a subspace of H¹(Ω) consisting of continuous piecewise affine functions on the mesh T_h, and

(9)

Then, the standard FEM is to find (u_h, B_h)∈W_0B^h (Ω) such that

(3)

The existence and uniqueness of solutions to (2) and (3) can be easily obtained by the wellknown Lax-Milgram lemma. A priori error analysis can be found in Ref. [21].

3 Analytical tools

In order to treat anisotropic elements, some analytical tools in the anisotropic setting have to be introduced here, which are taken from Refs. [8] and [9].

3.1 Inverse inequalities

Inverse inequalities for bubble functions are very important in deriving lower error bounds. As usual, we first introduce the bubble functions^[7]. For an arbitrary triangle K, denote the barycentric coordinates by λ_{K, 1}, λ_{K, 2}, and λ_{K, 3}. We define the element bubble function b_K by

(10)

Let σ be an inner edge shared by K₁ and K₂. We here enumerate all the vertices of K₁ and K₂ such that the two vertices of σ are first numbered. We define the edge bubble function b_σ by

(11)

Obvious modifications need to be done for boundary edges. For simplicity, b_K and b_σ are assumed to be extended by zero outside their original domain of definition.

For a given edge σ of a triangle K, an extension operator F_ext: P₀(σ) → P₀(K) is defined as follows:

(12)

By standard scaling arguments, we can easily derive the following anisotropic inverse inequalities.

Lemma 1 Assume that ϕ_K∈P₀(K) and ϕ_σ∈P₀(σ), where σ ⊂ ∂K. Then,

(13)

3.2 Alignment measure and anisotropic interpolation estimates

From a heuristic point of view, if the (directional) derivative of the solution displays little change in some direction, in that direction, we should stretch the mesh elements. The better the alignment is between the anisotropic mesh and the anisotropy of solution, the more accurate the error estimates can be. All practical applications intuitively follow this concept. Next, we should introduce the alignment measure m₁(v, T_h) in order to quantify this alignment.

For v∈H¹(Ω) and a family of triangulations {T_h}, the alignment measure m₁: H¹(Ω) × {T_h} → R is defined by

(4)

The alignment measure has the following property:

(14)

The above property implies that, if the mesh T_h is well aligned with an anisotropic function v, it leads to a small alignment measure. In practice, for sensible anisotropic meshes, one almost always obtains m₁(v, T_h) ~ 1. In this paper, we should apply this property to anisotropic mesh adaption (see numerical experiments in Section 5).

Consider a node a∈N_h. The local L²-projection P_a: H¹(ω_a) → P₀(ω_a) is uniquely defined by . Then, the Clément interpolation operator I₀: H₀¹(Ω) → V^h(Ω) ∩ H₀¹(Ω) is defined by

(15)

and I_B: H_B¹ (Ω) → V^h(Ω) ∩ H_B¹ (Ω) by

(16)

where λ_a is the (piecewise linear) basis function related to the vertex a, and N_N^ext := {a∈N_h^ext; a∈Γ_N}.

Finally, we state the anisotropic interpolation error estimates based on the alignment measure.

Lemma 2 For all functions v∈H₀¹(Ω) and Q∈H_B¹ (Ω), there hold

(5)

(6)

(7)

(8)

where E_N^ext := {σ∈E_h^ext; σ ⊂ Γ_N}.

4 Residual error estimation

In this section, we present the residual error estimator for the conforming approximation (3) and show that it provides two-sided bounds on the error ||(u − u_h, B − B_h)||_W on anisotropic meshes.

Let (u_h, B_h) be the conforming finite element approximation. Here, we define the element residuals over an element K by

(17)

and edge residuals by

(18)

Obviously, R_{1, K} = 1+Ha a·▽B_h and R_{2, K} = Ha a·▽u_h hold for piecewise linear functions as considered here.

The local residual error estimators η_{1, K} and η_{2, K} for an element K are defined by

(9)

(10)

and the global terms are

(19)

separately.

Theorem 3 Let (u, B) be the weak solution of problem (1), and let (u_h, B_h) be the corresponding finite element approximation defined by (3). Then, the error is bounded globally from above by

(11)

Proof For the convenience of expression, we set e_u = u − u_h and e_B = B − B_h. Since (a · ▽e_B, e_u) = −(a · ▽e_u, e_B), due to integration by parts, the Galerkin orthogonality and Lemma 2 yield

(20)

which concludes the proof.

Theorem 4 Under the assumption of Theorem 3, the error is bounded locally from below by

(12)

(13)

for all K∈T_h.

Proof The proof of lower error bounds needs to use the bubble functions and corresponding anisotropic inverse inequalities (see Lemma 1). Here, we only need to prove the lower bound (12). The remaining one can be derived analogously.

Observe R_{1, K} = 1 + Δu_h + Ha a · ▽B_h∈P₀(K) and set

(21)

The Cauchy-Schwarz inequality and integration by parts yield

(14)

The inverse inequalities of Lemma 1 lead to

(22)

which, together with (14), result in

(15)

For an inner edge σ∈E_h^int, set

(23)

Integration by parts and (2) for (v, Q) = (w_σ, 0) implies that

(16)

Due to Lemma 1, we have the following relations:

(24)

which, together with (15) and (16), lead to

(17)

For an edge on the Dirichlet boundary, nothing needs to be done because R_{1, σ} = 0. Hence, combining (15) and (17), we conclude the lower error bound (12).

Remark 1 From Theorems 3 and 4, the lower bounds contain additional L²-error terms Ha||u−u_h||_ωK and Ha||B−B_h||_ωK that do not appear in the upper error bound (11). Therefore, the two-sided bounds on the error do not correspond completely. We point out that a very similar situation can be found for error estimators for the convection-diffusion problem^{[10, 22-23]}. Moreover, only if the L²-error terms is dominated by the H¹-error terms ||▽(u − u_h)||_ωK and ||▽(B − B_h)||_ωK, the upper and lower bounds on the error will be of the same quality. Hence, the additional L²-error terms are mainly due to the H¹-error terms, which is confirmed by our numerical experiments in the next section.

5 Numerical experiments 5.1 Anisotropic error indicator

Due to the large aspect ratio of anisotropic mesh, the alignment measure m₁(·, T_h) enters the upper error bound in order to measure how well the alignment between the solution and the mesh does in the sense: the worse the alignment is, the larger the alignment measure is. In other words, a sufficiently good mesh makes sure that the alignment measure is O(1), i.e.,

(18)

In turn, the relation (18) can be used as a guide to design a proper anisotropic mesh. To this end, the definition (4) of the alignment measure m₁(v, T_h) is reformulated by

(19)

where r₁ and r₂ are the corresponding unit vectors of p₁ and p₂ on each K∈T_h, respectively, and so C_K = (p₁, p₂) = (h_{max, Kr1}, h_{min, Kr2}). Therefore, the relation (18) holds, provided that

(20)

This formulation can be used to indicate the direction of adaptive mesh refinement. More specifically, due to (19), the upper error bound (11) holds as long as we have the following estimate:

(21)

where

(25)

The relation (20) implies η⁽¹⁾ ~ η⁽²⁾. However, the estimate (21) is not available, because the exact solution (u, B) is usually unknown. To overcome this, we use the well-known Z-Z estimate^[6]. The gradient (▽u, ▽B) of the exact solution can be replaced by a recovered gradient (▽^Ru, ▽^RB) by local averaging such as an approximate L²-projection of (▽u_h, ▽B_h) onto V^h × V^h, i.e.,

(26)

Numerical experiments in Ref. [24] showed that the averaging technique was quite more reliable on anisotropic meshes than expected. After this replacement, we introduce the global anisotropic error indicator and the local anisotropic error indicator by

(27)

In conclusion, the relation (18) implies that the adaptive algorithm should be designed such that the local anisotropic error indicators and satisfy, on each triangle K,

(28)

5.2 Adaptive algorithm

We should implement an adaptive algorithm using the software FREEFEM++^[18] with the mesh generator Bamg^[19] based on a right metric field which will be built later. The goal of this adaptive algorithm is to generate a suitable anisotropic mesh so that the relative estimated error is close to a preset tolerance ϕ_TOL, i.e.,

(29)

A sufficient condition to generate such an anisotropic mesh is to make sure that for all K∈T_h, it holds

(22)

where |T_h| denotes the number of elements in T_h, as well as the number |N_h| of mesh vertices which will be used later. In order to get data at the mesh vertices, at each vertex a, we introduce the anisotropic error indicator defined by

(30)

Since

(31)

there holds (22), whenever we generate a mesh satisfying, for all vertices a∈N_h,

(32)

With the mesh generator Bamg, the metric field M is needed to be given at each vertex a as continuous P₁ finite element functions, namely, the two eigenvalues λ_{1, a}, λ_{2, a} (λ_{1, a} ≤ λ_{2, a}) and the corresponding eigenvectors r_{1, a}, r_{2, a} such that M = R^TΛR, where

(33)

Then, in a vicinity of the vertex a, the wanted mesh size can be defined by the metric field M such that the size h is equal to |x|√x^TM(a)x in the direction x∈R², where |x| := √x^Tx.

Then, the adaptive algorithm is designed. First, the values and at all the vertices a of the mesh are defined by

(34)

Thus, the error at each vertex a in the direction of maximum stretching can be represented by , whereas the error at each vertex a in the direction of minimum stretching is denoted by . Second, we compute the value of mesh size h_{1, a} (h_{2, a}) for all vertices a∈N_h by averaging all the values h_{max, K} (h_{min, K}) corresponding to the neighboring triangles K in T_a. The algorithm is carried out as follows. For i = 1, 2, if

(35)

then the value of λ_{i, a} is set to . If

(36)

then the value of λ_{i, a} is set to be . Otherwise, λ_{i, a} is set to be h_{i, a}⁻¹. Finally, for all vertices a∈N_h, let r_{1, a} and r_{2, a} be, respectively, the unit eigenvectors corresponding to the smallest and largest eigenvalues of the Z-Z estimates of the error gradient matrix defined by an average

(37)

With the resulting metric M, we build a new anisotropic mesh using Bamg, which should insure m₁^R (u − u_h, T_h) ~ 1 and m₁^R (B − B_h, T_h) ~ 1, where m₁^R (u − u_h, T_h) and m₁^R (B − B_h, T_h) are the Z-Z estimates of m₁(u − u_h, T_h) and m₁(B − B_h, T_h) obtained by replacing the gradient (▽u, ▽B) of the exact solution by the recovered one (▽^Ru, ▽^RB), i.e.,

(38)

respectively. According to the numerical results of the next subsection, it always holds that both m₁^R (u − u_h, T_h) and m₁^R (B − B_h, T_h) are close to 1, which indicates that the anisotropy of meshes is well adapted to the exact solution of the problem. Therefore, an efficient error estimation is to be expected. In the future, we plan to extend such an efficient anisotropic adaptive FEM to other problems, e.g., nonconforming FEMs^[25-26] for fourth-order elliptic singular perturbation problems with boundary layers^[27].

5.3 Numerical results

In this subsection, for the three test problems, a series of numerical simulations are performed by the anisotropic adaptive FEM developed in this paper. In all the test problems, we consider the domain Ω = (0, 1) × (0, 1) and set the tolerance ϕ_TOL = 0.01.

Example 1 (The Shercliff problem) For this example, there exists an analytic solution, and the numerical data are available in Refs.[1] and [2] such that a comparison can be done between the approximate solution and the exact solution. Here, the walls of the channel are insulated (B = 0 on ∂Ω), and on the solid walls, the velocity is zero (u = 0 on ∂Ω). The external magnetic field B₀ is applied in the x-axis (α = 0) (see Fig. 2).

For the Hartmann numbers Ha = 100 and Ha = 500, the comparisons are done between the approximate solution (u_h, B_h) and the exact solution (u, B) in Tables 1 and 2, respectively. One can see that two results are comparable. In addition, Table 3 presents the corresponding mesh information and Z-Z estimates of the alignment measures after some iterations during the adaptation process. For the sake of comparison, the corresponding meshes and contour plots for approximate solutions are shown in Figs. 3-4 for Ha = 100 and Ha = 500, respectively. It is worth mentioning that the use of anisotropic adaptive FEM can capture the boundary or interior layers in some directions and allow the number of vertices dramatically reduced for a given level of accuracy. In the next two examples, the superior performance will be further demonstrated for much larger Hartmann numbers, such as 10³, 10⁴, and 10⁵.

Example 2 (The two-dimensional square-channel flow with an oblique applied magnetic field) The same MHD problem as Example 1 is considered, except that the externally applied magnetic field B₀ has various positive angles α with the x-axis (see Fig. 6). Here, for the values of α = 0, π/4, and π/3, numerical computations are carried out at different Hartmann numbers. The resulting meshes and contour plots for approximate solutions generated by the anisotropic adaptive FEM are presented at Hartmann numbers Ha = 10³, 10⁴, and 10⁵, respectively, in Figs. 7-15. The numerical results show that when Hartmann number increases, the boundary or interior layers are well dealt with by anisotropic meshes. The corresponding mesh information and approximate alignment measures are presented in Tables 4-6.

Example 3 (The two-dimensional square-channel flow with a partly conducting boundary) Here, the external magnetic field B₀ is perpendicular to the wall at x = −1 (i.e., α = 0). For a length 2l at the center, the wall is electrically insulated (see Fig. 16). In Figs. 17-25, the anisotropic meshes and contour plots for approximate solutions are presented for l = 0.2, 0.5, and 0.7 at different Hartmann numbers (Ha = 10³, 10⁴, and 10⁵), respectively. The corresponding mesh information and approximate alignment measures are presented in Tables 7-9.

6 Conclusions

In this paper, an anisotropic adaptive FEM is proposed to solve the MHD duct flow at high Hartmann numbers. The most distinguish feature of this method is that the layer information from some directions is captured well such that the number of mesh vertices is dramatically reduced for a given level of accuracy. Further, an adaptive algorithm is implemented for several examples on a rectangular domain, where it is shown that even for high values of the Hartmann number, this method gives accurate and stable results. Finally, we note that the anisotropic adaptive FEM can be applied to channels with arbitrary cross-sections.