3-D Finite-Element Modelling of a Maglev System using Bulk High-Tc Superconductor and its Application

1. Introduction

Maglev using bulk High-Tc superconductor can realize stable levitation without any active control (Brandt, 1989), and this facsinating property can reduce remarkably the complexity of the levitation system and therefore exhibits promising application in several fields such as maglev bearing (Hull, 2000; Ma, et al., 2003) and maglev transit (Wang, et al., 2002 ; Wang, J. & Wang, S., 2005; Schultz, et al., 2005; Sotelo, et al., 2010). To understand the electromagnetic interaction between the bulk high-Tc superconductor and its applied fields generated by various permanent magnetic devices and to provide a numerical tool to conduct the design for practical application, many methods have been proposed to numerically estimate the characteristics of the magnetic force of the bulk high-Tc superconductor.

The earliest method (Davis, et al., 1988) after the discovery of the high-Tc superconductor was basically established on the traditional mirror-image-model which uses Bean’s critical model (Bean, 1964). However, this model and the later frozen-image model (Kordyuk, 1998) can not reflect the important hysteresis property of the levitation force (Hull & Cansiz, 1999). Though this demerit can be overcome to some extent by introducing additional image dipoles (Yang & Zheng, 2007), this kind of model is essentially a phenomenological one and its applicable scope is also confined to miniature scale systems due to the essential dipole approximation of the levitated body in deducing the model.

Based on the principle of minimum energy, the current distribution in the high-Tc superconductor can be acquired by an iterative process, and then the magnetic force of the high-Tc superconductor can be calculated by Lorentz equation (Sanchez & Navau, 2001; Navau & Sanchez, 2001; Sanchez, et al., 2006). However, up to date, this method is used only to investigate the axisymmetric system with cylindrical high-Tc superconductor and permanent magnet (Navau & Sanchez, 2001) or 2-D translational symmetry system with rectangular high-Tc superconductor and permanent magnetic array (Sanchez, et al., 2006).

According to the state variables used in the governing equations, the methods describing the electromagnetic property of the high-Tc superconductor based on Maxwell’s equations can be classified into three types, i.e., A-V-method (Prigozhin, 1997), T-method (Hashizume, et al., 1991), and H-method (Pecher, et al., 2003). Most previous work (Qin, et al., 2002; Uesaka, et al., 1993; Yoshida, et al., 1994; Luo, et al., 1999; Alonso, et al., 2004; Gou, et al.,2007a; Gou, et al.,2007b) using these methods was to calculate the levitation force of a bulk high-Tc superconductor by aid of the assumption that the induced current only flows in the ab-plane of the high-Tc superconductor. In this case, only the component of the state variable along the c-axis of the high-Tc superconductor is considered in the governing equation and the number of degrees of freedom is therefore reduced and the problem to be solved is actually a 2-D one. This assumption is acceptable when the studied problem has an axisymmetric geometry and the movement of the bulk high-Tc superconductor is restricted to the vertical direction above the center of the magnetic device because in this situation, the induced current in the high-Tc superconductor due to the variation of the applied field will flow along the plane parallel to the ab-plane. As a result, the numerical results of the levitation force compare well with the measured data (Uesaka, et al., 1993; Alonso, et al., 2004; Gou, et al.,2007a), but for other situations, e.g., a bulk high-Tc superconductor is above a rectangular permanent magnet or a magnetic rail, the assumption is no longer satisfied and modelling high-Tc superconductor in full 3-D case is thereby required. Based on the A-V-method, the 3-D model has been proposed to numerically estimate the characteristics of the levitation force as well as the lateral force of a levitating transporter using bulk high-Tc superconductor (Ueda, et al., 2006), and also the dynamic behavior of the levitation system composed of a rectangular high-Tc superconductor and permanent magnet (Alloui, et al., 2009). However, though an anisotropic critical model is employed in the numerical iterative proces (Ueda, et al., 2006), none of them has considered the special electromagnetic anisotropic behavior of the high-Tc superconductor in deducing the governing equations. Those models are hence still an anisotropic 3-D model. For the maglev transit using bulk high-Tc superconductor above a magnetic rail, a 3-D model considering the anisotropic behavior is reported (Lu, et al., 2008). In their model, in order to describe the anisotropic behavior, the bulk high-Tc superconductor is split into two different parts: One is a homogeneous part that is identical to an anisotropic superconductor, and the other is considered as a conductor whose conductivity only along the c-axis is not equal to zero. However, the reasonable explanation of this superposition approach to describe the anisotropic behavior in physics is still a question needed to be answered at present.

In this chapter, we report a 3-D finite-element model using current potential T as state variable. In this model, the anisotropic behavior of the high-Tc superconductor is contained in the 3-D governing equations by considering a tensor resistivity, and the finite-element technique is empolyed to numerically solve the mathmatical formulations on a VC++ software platform. The numerical results of both levitation force and lateral force were validated by the measured data. Lastly, one example using this 3-D finite-element model to optimize the magnetic rail is introduced to present its viable use for practical design of maglev system using bulk high-Tc superconductor.

2. Mathematical formulations

2.1. Formulations to model the anisotropy in high-Tc superconductor

The special microstructure, which consists of the alternating stack of superconductive CuO₂ layers and almost insulating block layers, results in a remarkable anisotropic behavior in the present high-Tc superconductor (Dinger, et al., 1987). Due to this anisotropic behavior arising from the intrinsic pinning and other defects in the high-Tc superconductor, the flux-line curvature will always occur when the high-Tc superconductor is placed in a magnetic field.

Figure 1.

Schematic drawing of the elliptical model (left) and the Cartesian coordinate system with a bulk high-Tc superconductor whose c-axis is parallel to the z-axis (right).

Consequently, the critical current density J_c of the high-Tc superconductor is anisotropic and strongly dependent on the orientation of the applied field (Mikitik & Brandt, 2000). Namely,the value of the critical current density flowing in the ab-plane Jab c is larger than that along the c-axis Jc c (Matsushita, 2007).

Many methods have been proposed to formulate the dependence of the critical current density J_c on the angle φ between the orientation of the local applied field and c-axis. However, a few of them are difficult to be employed in the present calculation due to its complexity in determining the parameters involved in the formulation (Mikitik & Brandt, 2000) or its applicable field range is beyond our focus (Sawamura & Tsuchimoto, 2000). Another model reported in (Yang, et al., 1999) is simple to be employed and has been valiadated in an axisymmetric levitation system. This model is desribed in the following:

When the out-of-plane anisotropic ratio α of the critical current density J_c is defined as α = Jab c/Jc c, the J_c(φ) relation can be rewritten as follows:

$$J_c\left(\phi \right)=J_c^{ab}\left({\cos }^2\phi +{\sin }^2\phi /\alpha \right)$$E1

Here, we will introduce an elliptical model, which has been used to investigate the electromagnetic problem involving the anisotropic ferromagnetic material (Napoli & Paggi, 1983), to describe the angular-dependence property of the critical current density J_c in the high-Tc superconductor. According to the schematic drawing shown in Fig. 1, the elliptic model can be expressed by the following equation:

$${\left(\frac{J_{cx}}{J_c^{ab}}\right)}^2+{\left(\frac{J_{cz}}{J_c^c}\right)}^2=1$$E2

where J_cx and J_cz are the induced current densities in the ab-plane and parallel to the c-axis, respectively. According to (2), the J_c(φ) relation with respect to Jab c and α can be expressed as:

$$J_c\left(\phi \right)=\sqrt{J_{cx}^2+J_{cz}^2}=J_c^{ab}\sqrt{{\cos }^2\phi +{\left(\sin \phi /\alpha \right)}^2}$$E3

The resistivity of the high-Tc superconductor is also anisotropic (Wu, et al., 1991) and can be represented by a tensor while modelling the high-Tc superconductor. The tensor of the resistivity of the high-Tc superconductor can be reduced to a following diagonal matrix when only the out-of-plane anisotropy is considered,

$$\overline{\overline{{\rho }_s}}=\left[\begin{array}{l}{\rho }_{ab}00\\ 0{\rho }_{ab}0\\ 00{\rho }_c\end{array}\right]\text{=}\left[\begin{array}{l}{\rho }_{ab}00\\ 0{\rho }_{ab}0\\ 00\alpha {\rho }_{ab}\end{array}\right]$$E4

where ρ_ab and ρ_c are the resistivity in the ab-plane and along the c-axis, respectively. It is important to remark here that, for the identical electrical field criterion E_c to determine the critical current density J_c, the anisotropic ratio of the resistivity ρ and the critical current density J_c are reciprocal to each other. The above tensor expression of the resistivity illustrates a possible way to establish the governing equations of the high-Tc superconductor including its anisotropic behavior.

2.2. 3-D governing equations

The T-method is adopted in our model because it has the merits that, the number of unknown and space needed to be meshed can be remarkably reduced because the variable T can be defined as zero outside the conductor (Miya & Hashizume, 1988).

As a conducting material, Maxwell’s equations are also valid to describe the electromagnetic phenomena in the high-Tc superconductor. Thus, we have,

$$\nabla \times H=J$$E5

$$\nabla \times E=-\frac{\partial B}{\partial t}$$E6

According to (4), the current density J is a divergence free vector with a quasi-static approximation in the low-frequency problem, i.e., $$\nabla ·J=0,$$and thus, a current vector potential T can be introduced and defined as,

$$J=\nabla \times T$$E7

The Coulomb gauge is applied to vector T to guarantee the uniqueness of the solution, i.e., $$\nabla ·T=0.$$Applying Helmholtz’s theorem to vector T yields the following equation:

$$\begin{array}{l}C\left(P\right)T\left(P\right)=\frac{1}{4\pi }{\displaystyle {\int }_V\left({\nabla }^{\prime }·T\left(P^{\prime }\right)\right){\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{V}^{\prime }}-\frac{1}{4\pi }{\displaystyle {\int }_S\left(n^{\prime }·T\left(P^{\prime }\right)\right){\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}\\ +\frac{1}{4\pi }{\displaystyle {\int }_V\left({\nabla }^{\prime }\times T\left(P^{\prime }\right)\right)\times {\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{V}^{\prime }}-\frac{1}{4\pi }{\displaystyle {\int }_S\left(n^{\prime }\times T\left(P^{\prime }\right)\right)\times {\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}\end{array}$$E8

where R(P, P′ ) is the distance between the source point P′ and the field point P, the superscript ′ refers to the quantity at the source point, n′ is a unit vector out of the surface S′, and the coefficient C(P) takes the following values (Hashizume, et al., 1991)：

$$C\left(P\right)=\{\begin{array}{l}1P\in V^{\prime }(excluding{S}^{\prime })\\ 1/2P\in S^{\prime }\\ 0elsewhere\end{array}$$E9

According to the physical fact that the normal component of J must be zero on the surface of the high-Tc superconductor, i.e., J_n = 0, T has the following boundary condition (Miya & Hashizume, 1988):

$$n^{\prime }\times T=0$$E10

Therefore, only the normal component T_n exists on all the surfaces of the high-Tc superconductor, e.g., for the surface of AA′D′D as shown in Fig.1, vector T is reduced to a scalar T_x. Furthermore, T is zero on all edges of the rectangular high-Tc superconductor bulk, e.g., for edge AA′ as shown in Fig.1, when it is regarded as a part of surface AA′D′D, T_y = T_z = 0, whereas when it is considered as a part of surface AA′B′B, T_x = T_z = 0. Thus, we have T_x = T_y = T_z = 0 on all the edges.

Equation (7)is reduced to the following form when the Coulomb gauge and boundary condition (8) are considered,

$$C\left(P\right)T\left(P\right)=\frac{1}{4\pi }{\displaystyle {\int }_V\left({\nabla }^{\prime }\times T\left(P^{\prime }\right)\right)\times {\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{V}^{\prime }}-\frac{1}{4\pi }{\displaystyle {\int }_S\left(n^{\prime }·T\left(P^{\prime }\right)\right){\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}$$E11

The B-H constitutive law of the high-Tc superconductor can be assumed to be linear as that in vacuum with a good approximation because its applicable conditions (Brandt, 1996) can be easily satisfied in a levitation system using bulk Y-Ba-Cu-O due to its small lower critical field B_c1 (Krusin-Elbaum, et al., 1989) and large applied field as well as geometry. Thus,

$$B={\mu }_{0}H$$E12

The induced field B_s produced by the induced current in the high-Tc superconductor can be expressed in terms of the vector T in the following equation when we combine (9) with Biot-Savart’s law (Miya & Hashizume, 1988),

$$B_s={\mu }_{0}C\left(P\right)T\left(P\right)+\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_S\left(n^{\prime }·T\left(P^{\prime }\right)\right){\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}$$E13

When an equivalent conductivity σ_s, which is nonlinear and dependent on the local electrical field of the high-Tc superconductor, is introduced, the traditional Ohm’s law in the high-Tc superconductor has the following form:

$$J={\sigma }_s\left(\left|E\right|\right)E$$E14

By substituting (6) and (12) into (5) and considering B = B_e + B_s where B_e is the applied field, yields the following equation:

$$\nabla \times \frac{1}{{\sigma }_s}\left(\nabla \times T\right)+\frac{\partial \left(B_e+B_s\right)}{\partial t}=0$$E15

The governing equation of the high-Tc superconductor based on the variable T is finally derived from (11), and (13) as follows,

$$\nabla \times \frac{1}{{\sigma }_s}\left(\nabla \times T\right)+{\mu }_{0}C\left(P\right)\frac{\partial T}{\partial t}+\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_S\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}{\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}+\frac{\partial B_e}{\partial t}=0$$E16

When we replace $${\sigma }_s$$by a tensor resistivity$$\overline{\overline{{\rho }_s}}$$, (14) and consider that,

$$\overline{\overline{{\rho }_s}}\left(\nabla \times T\right)=\left[\begin{array}{l}{\rho }_{ab}00\\ 0{\rho }_{ab}0\\ 00\alpha {\rho }_{ab}\end{array}\right]\left[\begin{array}{l}\left(\frac{\partial T_z}{\partial y}-\frac{\partial T_y}{\partial z}\right)\widehat{x}\\ \left(\frac{\partial T_x}{\partial z}-\frac{\partial T_z}{\partial x}\right)\widehat{y}\\ \left(\frac{\partial T_y}{\partial x}-\frac{\partial T_x}{\partial y}\right)\widehat{z}\end{array}\right]={\rho }_{ab}\left[\left(\frac{\partial T_z}{\partial y}-\frac{\partial T_y}{\partial z}\right)\widehat{x}+\left(\frac{\partial T_x}{\partial z}-\frac{\partial T_z}{\partial x}\right)\widehat{y}+\alpha \left(\frac{\partial T_y}{\partial x}-\frac{\partial T_x}{\partial y}\right)\widehat{z}\right]$$E17

The following equality for the first term in the left side of (14) is satisfied,

$$\begin{array}{l}\nabla \times \frac{1}{{\sigma }_s}\left(\nabla \times T\right)=\nabla \times \overline{\overline{{\rho }_s}}\left(\nabla \times T\right)={\rho }_{ab}\left[\alpha \frac{{\partial }^2{T}_y}{\partial x\partial y}-\alpha \frac{{\partial }^2{T}_x}{\partial y^2}-\frac{{\partial }^2{T}_x}{\partial z^2}+\frac{{\partial }^2{T}_z}{\partial x\partial z}\right]\widehat{x}\\ +{\rho }_{ab}\left[\frac{{\partial }^2{T}_z}{\partial y\partial z}-\frac{{\partial }^2{T}_y}{\partial z^2}-\alpha \frac{{\partial }^2{T}_y}{\partial x^2}+\alpha \frac{{\partial }^2{T}_x}{\partial y\partial x}\right]\widehat{y}+{\rho }_{ab}\left[\frac{{\partial }^2{T}_x}{\partial z\partial x}-\frac{{\partial }^2{T}_z}{\partial x^2}-\frac{{\partial }^2{T}_z}{\partial y^2}+\frac{{\partial }^2{T}_y}{\partial z\partial y}\right]\widehat{z}\end{array}$$E18

According to the Coulomb gauge, we have$$\nabla \left(\nabla ·T\right)=0$$ i.e.,

$$\left(\frac{{\partial }^2{T}_x}{\partial x^2}+\frac{{\partial }^2{T}_y}{\partial y\partial x}+\frac{{\partial }^2{T}_z}{\partial z\partial x}\right)\widehat{x}+\left(\frac{{\partial }^2{T}_x}{\partial x\partial y}+\frac{{\partial }^2{T}_y}{\partial y^2}+\frac{{\partial }^2{T}_z}{\partial z\partial y}\right)\widehat{y}+\left(\frac{{\partial }^2{T}_x}{\partial x\partial z}+\frac{{\partial }^2{T}_y}{\partial y\partial z}+\frac{{\partial }^2{T}_z}{\partial z^2}\right)\widehat{z}=0$$E19

The following identities can be derived from (17),

$$\frac{{\partial }^2{T}_z}{\partial z\partial x}=-\frac{{\partial }^2{T}_x}{\partial x^2}-\frac{{\partial }^2{T}_y}{\partial y\partial x},$$$$\frac{{\partial }^2{T}_z}{\partial z\partial y}=-\frac{{\partial }^2{T}_x}{\partial x\partial y}-\frac{{\partial }^2{T}_y}{\partial y^2},$$

$$\frac{{\partial }^2{T}_x}{\partial x\partial z}+\frac{{\partial }^2{T}_y}{\partial y\partial z}=-\frac{{\partial }^2{T}_z}{\partial z^2}$$E20

Equality (16) can be further written in the following way when (18) is taken into account,

$$\begin{array}{l}\nabla \times \overline{\overline{{\rho }_s}}\left(\nabla \times T\right)={\rho }_{ab}\left[-\frac{{\partial }^2{T}_x}{\partial x^2}-\alpha \frac{{\partial }^2{T}_x}{\partial y^2}-\frac{{\partial }^2{T}_x}{\partial z^2}+\left(\alpha -1\right)\frac{{\partial }^2{T}_y}{\partial x\partial y}\right]\widehat{x}\\ +{\rho }_{ab}\left[-\alpha \frac{{\partial }^2{T}_y}{\partial x^2}-\frac{{\partial }^2{T}_y}{\partial y^2}-\frac{{\partial }^2{T}_y}{\partial z^2}+\left(\alpha -1\right)\frac{{\partial }^2{T}_x}{\partial x\partial y}\right]\widehat{y}+{\rho }_{ab}\left(-\frac{{\partial }^2{T}_z}{\partial x^2}-\frac{{\partial }^2{T}_z}{\partial y^2}-\frac{{\partial }^2{T}_z}{\partial z^2}\right)\widehat{z}.\end{array}$$E21

Besides,

$${\nabla }^{\prime }\frac{1}{R\left(P,P^{\prime }\right)}=\frac{\partial }{\partial x^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)\widehat{x}+\frac{\partial }{\partial y^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)\widehat{y}+\frac{\partial }{\partial z^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)\widehat{z},$$E22

and

$$B_e\text{=}{B}_{ex}\widehat{x}+B_{ey}\widehat{y}+B_{ez}\widehat{z}.$$E23

Finally, the following 3-D governing equations are obtained once (19) to (21) are substituted into (14) and ρ_ab is replaced by its reciprocal σ_ab.

$$\frac{1}{{\sigma }_{ab}}\left(-\frac{{\partial }^2{T}_x}{\partial x^2}-\alpha \frac{{\partial }^2{T}_x}{\partial y^2}-\frac{{\partial }^2{T}_x}{\partial z^2}+\left(\alpha -1\right)\frac{{\partial }^2{T}_y}{\partial x\partial y}\right)+{\mu }_{0}C\left(P\right)\frac{\partial T_x}{\partial t}+\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_S\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}\frac{\partial }{\partial x^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)d{S}^{\prime }}+\frac{\partial B_{ex}}{\partial t}=0$$E24

$$\frac{1}{{\sigma }_{ab}}\left(-\alpha \frac{{\partial }^2{T}_y}{\partial x^2}-\frac{{\partial }^2{T}_y}{\partial y^2}-\frac{{\partial }^2{T}_y}{\partial z^2}+\left(\alpha -1\right)\frac{{\partial }^2{T}_x}{\partial x\partial y}\right)+{\mu }_{0}C\left(P\right)\frac{\partial T_y}{\partial t}+\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_S\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}\frac{\partial }{\partial y^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)d{S}^{\prime }}+\frac{\partial B_{ey}}{\partial t}=0$$E25

$$\frac{1}{{\sigma }_{ab}}\left(-\frac{{\partial }^2{T}_z}{\partial x^2}-\frac{{\partial }^2{T}_z}{\partial y^2}-\frac{{\partial }^2{T}_z}{\partial z^2}\right)+{\mu }_{0}C\left(P\right)\frac{\partial T_z}{\partial t}+\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_S\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}\frac{\partial }{\partial z^{\prime }}\left(\frac{1}{R\left(P,P^{\prime }\right)}\right)d{S}^{\prime }}+\frac{\partial B_{ez}}{\partial t}=0$$E26

where σ_ab is the conductivity in the ab-plane. It is worth noting that, compared with the traditional T-Ω method, the complexity of the governing equations is reduced due to the omission of the variable Ω (Miya & Hashizume, 1988), and thus the number of varibles in the problem is three with unknown T_x, T_y, and T_z.

2.3. Nonlinear E-J relations

Typically, there are three different models to address the highly nonlinear relationship between E and J of the high-Tc superconductor, i.e., Bean’s critical current model (Bean, 1964), power law model (Rhyner, 1993) and flux flow and creep model (Yamafuji & Mawatari, 1992). For Bean’s critical current model, it fails to investigate problems such as force relaxation (Luo, et al., 1999) and drift under vibration (Gou, et al., 2007) due to the assumption that the current flowing in the superconductor is constant with time and the lack of material related parameters in its model. In addition, Bean’s critical current model can be considered as an infinite case of the power law model (Rhyner, 1993). This model is thereby not employed in the following calculation.

When an index n is introduced and defined as n=U₀/kΘ where U₀ is the pinning potential of the superconductor at an absolute temperature Θ and k is the Boltzmann constant, the power law model is expressed as:

$$E=E_c{\left(\left|J\right|/J_c\right)}^{n}J/\left|J\right|$$E27

According to Anderson theory (Anderson, 1962), the responding behavior between E and J due to flux flow and creep model phenomena in the superconductor can be described as:

$$E=\{\begin{array}{l}2{\rho }_c{J}_c\sinh (\frac{U_0}{\kappa \Theta }\frac{\left|J\right|}{J_c})\exp (-\frac{U_0}{\kappa \Theta })\frac{J}{\left|J\right|}0\le \left|J\right|\le J_c\\ E_c+{\rho }_f{J}_c(\frac{\left|J\right|}{J_c}-1)\frac{J}{\left|J\right|}\left|J\right|>J_c\end{array}$$E28

where ρ_c and ρ_f are the creep and flow resistivities respectively, and $$E_c={\rho }_c{J}_c\left[1-\exp (-2U_0/\kappa \Theta )\right]\approx {\rho }_c{J}_c.$$

3. Numerical solution using finite-element technique

The previous work in handling linear eddy current problem has proven that Boundary Element Method is an effective method to numerically solve the governing equations including both differential and integral terms (Miya & Hashizume, 1988). Unfortunately, for the case studied here, the high-Tc superconductor is a highly nonlinear media, and this special property leads to difficulty when Boundary Element Method is employed to handle the 3-D governing equations. Consequently, the finite element method is a frequent choice to discretize the governing equation of the high-Tc superconductor in space (Pecher, et al., 2003; Uesaka, et al., 1993; Yoshida, et al., 1994; Luo, et al., 1999; Alonso, et al., 2004; Gou, et al.,2007a; Gou, et al.,2007b; Lu, et al., 2008). In particular, one format of the finite difference method, named Crank-Nicolson-θ method, is employed to numerically perform the time domain of the 3-D governing equations.

3.1. Finite-element matrices

The governing equations (22) to (24) incorporated with the boundary condition are numerically solved by finite-element technique via Galerkin’s method. The tetrahedral element is chosen to mesh the domain of the high-Tc superconductor. The final algebraic equations corresponding to (22) to (24) can be compactly expressed by the following matrix equation,

$$\left[K_i\left({\sigma }_{ab}\right)\right]\left\{T_i\right\}+\left[K_{12}\left({\sigma }_{ab}\right)\right]\left\{T_{3-i}\right\}+\left[Q_{0i}\right]\left\{\frac{\partial T_i}{\partial t}\right\}+\left[Q_{1i}\right]\left\{\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}\right\}=\left\{L_i\right\}(i=1,2,3)$$E29

where

i	Ki	Ti	T3-i	Q0i	Q1i	Li	αi	α3-i	u′i	Bei
1	Kx	Tx	Ty	Q0x	Q1x	Lx	1	α	x′	Bex
2	Ky	Ty	Tx	Q0y	Q1y	Ly	α	1	y′	Bey
3	Kz	Tz	0	Q0z	Q1z	Lz	1	1	z′	Bez

Table 1.

Specification of the variable parameters in (27)

$$\left[Q_{0i}\right]={\displaystyle \sum _e{\left[Q_{0i}\right]}_e}={\mu }_0{\displaystyle \sum _e{\displaystyle {\int }_{V_e}{\left[N\right]}_e^T{\left[C\left(P\right)N\right]}_{e}dV}}$$$$\left[Q_{1i}\right]={\displaystyle \sum _e{\left[Q_{1i}\right]}_e}=\frac{{\mu }_0}{4\pi }{\displaystyle \sum _e{\displaystyle {\int }_{V_e}{\left[N\right]}_e^T\left({\displaystyle {\int }_S{\left[N\right]}_e\frac{\partial }{\partial {u^{\prime }}_i}\frac{1}{R\left(P,P^{\prime }\right)}d{S}^{\prime }}\right)dV}}$$

$$\left\{L_i\right\}={\displaystyle \sum _e{\left\{L_i\right\}}_e}=-{\displaystyle \sum _e{\displaystyle {\int }_{V_e}{\left[N\right]}_e^T\frac{\partial B_{ei}}{\partial t}dV}}$$E30

N is the shape function of the linear tetrahedral nodal element. The subscript e represents its detailed formula in each element. The conductivity is different for each element and it is therefore represented by σe ab here. All the other parameters with a subscript including i can be determined by Table 1.

The finite-element matrices (27) for three components are integrated into one matrix in the numerical program, i.e.,

$$\left[K\left({\sigma }_{ab}\right)\right]\left\{T\right\}+\left[Q_0\right]\left\{\frac{\partial T}{\partial t}\right\}+\left[Q_1\right]\left\{\frac{\partial \left(n^{\prime }·T^{\prime }\right)}{\partial t}\right\}=\left\{L\right\}$$E31

where

$$\left[K\left({\sigma }_{ab}\right)\right]={\displaystyle \sum _{i=1}^2\left[\left[K_i\left({\sigma }_{ab}\right)\right]+\left[K_{12}\left({\sigma }_{ab}\right)\right]\right]}+\left[K_3\left({\sigma }_{ab}\right)\right]$$$$\left[Q_0\right]={\displaystyle \sum _{i=1}^3\left[Q_{0i}\right]}$$, $$\left[Q_1\right]={\displaystyle \sum _{i=1}^3\left[Q_{1i}\right]}$$,

$$\left\{L\right\}={\displaystyle \sum _{i=1}^3\left\{L_i\right\}}$$E32

After applying Crank-Nicolson-θ method to (28) for time discretization, the matrix at nth time step can be written as:

$$\left(\frac{\left[Q_{\text{0}}\right]+\left[Q_1\right]}{\Delta t}+\theta \left[K\left({\sigma }_{ab}\right)\right]\right)\left\{T^n\right\}=\theta \left\{L^n\right\}+\left(1-\theta \right)\left\{L^{n-1}\right\}+\left(\frac{\left[Q_0\right]+\left[Q_1\right]}{\Delta t}-\left(1-\theta \right)\left[K\left({\sigma }_{ab}\right)\right]\right)\left\{T^{n-1}\right\}$$E33

The disadvantage of the T-method is that the coefficient of (29) is a dense matrix. In the previous work, an over-relaxation iterative solution approach has been proposed to handle this dense matrix on a linear eddy current problem (Takagi, et al., 1988). However, to the situation we are facing, the nonlinearity of the E–J characteristic brings an additional iterative procedure in determining the conductivity of the high-Tc superconductor. Therefore, if the over-relaxation approach is employed in our computation, there would be two iterative procedures in the numerical program, and this would give rise to numerical instability in the calculation as well as the complexity of the numerical program. The dense feature of the coefficient matrix of (29) arises from the dense matrix [Q₁], which is related to the integral term of (27). Therefore, the coefficient matrix will be a sparse one if [Q₁] is transferred to the right side of (29), and this operation gives the following form,

$$\left(\frac{\left[Q_0\right]}{\Delta t}+\theta \left[K\left({\sigma }_{ab}\right)\right]\right)\left\{T^n\right\}=\theta \left\{L^n\right\}+\left(1-\theta \right)\left\{L^{n-1}\right\}+\left(\frac{\left[Q_0\right]}{\Delta t}-\left(1-\theta \right)\left[K\left({\sigma }_{ab}\right)\right]\right)\left\{T^{n-1}\right\}-\left[Q_1\right]\left\{T^{n-1}-T^{n-2}\right\}$$E34

In (30), the dense matrix [Q₁] is related to the difference of unknown variables between the last and the previous two time steps. The value of the unknown at the adjacent time step will approach to each other as the continuous decrease of the time step size △t. Therefore, in order to obtain a higher precision, we should assign a sufficient small time step size.

3.2. Nonlinear equation solution approach

The common Newton-Raphson method is employed to solve the nonlinear equations obtained from (30). Basically, the nonlinearity of the equations is eliminated by introducing a linear residual, and the corresponding linear equations can be integrated after calculating the Jacobi matrix for each element. To improve the stability of the calculation, a relaxation coefficient is also introduced as suggested in (Grilli, et al., 2005). The relaxation coefficient is assigned with an initial value before calculation. During the calculation, it will be reduced to a smaller one once the convergence can not be achieved within a threshold of the total iterative steps, and in this case the current time step is recalculated with the new value.

Because of the significant increase of the order of the coefficient matrix when the problem is extended from 2-D to 3-D and also the element of the matrix with the anisotropy consideration involved (an additional second derivative in term of x and y appears in (22) and (23) when the anisotropic behavior is taken into account), the Incomplete Cholesky-Conjugate Gradient method (Kershaw, 1978), which is regarded as an effective approach to solve linear algebraic equations with large symmetric positive definite matrix, is employed to solve the linear equations at each iterative step in an accelerated form. In this accelerated form, an accelerated factor is introduced in the incomplete decomposition process (Fujiwara, et al., 1993) after the coefficient matrix and column vector are preconditioned by a method proposed in (Cui, 1989).

3.3. Numerical procedure

The above-discussed numerical method is implemented by a self-written numerical program based on a VC++ platform with the following typical steps:

Step 1: Initial σ_ab in all elements is assumed to be the same and sufficiently large value before the first time step.

Step 2: Calculating [K(σ_ab)] with the present σ_ab in each element at the current time step, and vector T is obtained by solving (30).

Step 3: The σ_ab at the kth iterative step (σk ab) is calculated in each element according to the following method using either power law model or flux flow and creep model.

If J_c(φ) is expressed by (1),

$${\sigma }_{ab}^k=\left|J_k\right|/\left|E_k\right|/\left({\cos }^2\phi +{\sin }^2\phi /\alpha \right)$$E35

If J_c(φ) is expressed by (3),

$${\sigma }_{ab}^k=\left|J_k\right|/\left|E_k\right|/\sqrt{{\cos }^2\phi +{\left(\sin \phi /\alpha \right)}^2}$$E36

At the end, σ_ab in each element is replaced by its new value σk ab.

Step 4: Repeating steps 2 and 3 until the residual becomes less than a prescribed tolerance ε, i.e.,

$${\left[b-f\left(T^{\left(k\right)}\right)\right]}^2<\epsilon$$E37

where b is a column vector corresponding to the right side of (30), and f (T^(k)) is the result of the left side of (30) at the kth iterative step.

Step 5: The current density J is obtained by (6), and then the magnetic force is calculated by Lorentz equation,

$$F={\displaystyle {\int }_{V}J\times B_{e}dV}$$E38

Step 6: t = t + Δt until the maximum number of time step is achieved, and the steps 2-5 are repeated.

3.4. Numerical precision

Based on a levitation system composed of a bulk high-Tc superconductor（single-domain with a cubical shape）and magnetic rail（assembled by the permanent magnets with opposite magnetization direction as shown in the inset of Fig. 2）, the dependence of the computational precision of the levitation force on both mesh density and time step is discussed in this section. In the calculation, the bulk high-Tc superconductor was downward in a speed of 1 mm/s from the filed cooling position to the nearest gap, and then upward to its original position. On the basis of the geometrical and material parameters listed in Table 2 and power law model, the levitation force of the bulk high-Tc superconductor was calculated in this vertical down-and-up movement with different mesh densities and time steps.

wSC (mm)	lSC (mm)	tSC (mm)	Jab c(A/m2)	Ec (V/m)	U0 (eV)	α
10	10	10	2.5×108	1×10-4	0.1	3

Table 2.

Parameters of high-Tc superconductor used in the calculation for investigating the numerical precision.

Fig. 2 shows the dependence of the levitation force versus gap curve and the maximum levitation force F_max at the nearest gap on the total finite element node. This figure clearly illustrates that the levitation force and its hysteresis behavior will approach to a saturated state with the continuous increment of the total finite element node, or in other words, the continuous increment of the mesh density. It can be seen from the inset that the F_max continuously increases with the total finite element node but the increased rate is gradually

Figure 2.

Levitation force versus gap with different total nodes and the maximum levitation force F_max as a function of the total finite element node (inset).

Figure 3.

Levitation force versus gap with different number of time step and the maximum levitation force F_max as a function of the number of time step (inset).

reduced, e.g., the rate of the increment is about 57.1% when the total finite element node is raised from 27 to 64, but it is only 2.1% when the node is raised from 343 to 512. This indicates that the computed result is consistent with the finite-element thoery, i.e., the computed value will approach to the real value with the continuous increment of mesh density. Fig. 3 plots the levitation force versus gap curve with different number of time steps and the dependence of the F_max on the number of time step. It is obvious that the levitation force and its hysteresis behavior will also approach to a saturated state with the increment of the number of time step, e.g., the reduced rate of the F_max is less than 0.024% when the number of time step is raised from 1000 to 2000 and the value of the two cases is almost identical as shown in the inset of Fig. 3. This indicates that, the approach used to handle the coupling problem of the governing equations is applicable when the number of time step adopted is sufficiently large.

4. Experimental validation

The experimental validation is of great importance to support the development of any theoretical model and to confirm the practical application of the theoretical mode to the real world. Compared with the previous validations of the model, we will validate the 3-D finite-element model in a more comprehensive way, i.e., in which two different types of motion are considered, i.e., vertical movement (perpendicular to the surface of the magnetic rail) and transverse movement (parallel to the surface of the magnetic rail), are considered, and the associated magnetic force computed with our model are compared with the experimental data.

wSC (mm)	lSC (mm)	tSC (mm)	Jab c(A/m2)	Ec (V/m)	U0 (eV)	ρf (Ωm)	α
30	30	15	1.6×108	1×10-4	0.1	5×10-10	3

Table 3.

Parameters of the high-Tc superconductor used in the calculation for verifying the 3-D finite-element model.

4.1. Brief introduction of the experiment

The self-developed maglev measurement system was used to measure the magnetic force of the bulk high-Tc superconductor (Wang, 2010). The Y-Ba-Cu-O sample used here is of single-domain with a rectangular shape, and its geometrical parameters and photo as well as its schematic drawing are presented in Table 3 and Fig. 4 respectively. The geometric parameters and photos of the two magnetic rail demonstrators (Rail_1 and Rail_2) used to generate the applied field for high-Tc superconductor are also shown in Fig. 4.

In the experiment, the sample was firstly bath-cooled in a liquid nitrogen vessel for several minutes to to insure the superconductive state was established at a certain position above the center of the magnetic rail, and then, for the vertical movement, the sample was put downward to an expected nearest gap (distance from the bottom of the high-Tc superconductor to the upper surface of the magnetic rail) and then upward to its original position in the levitation case, whereas in the suspension case, the movement was opposite; For the transverse movement, the sample was put downward (levitation case) or upward (suspension case) to its levitation height, and then it was forced to move along the transverse direction parallel to the surface of the magnetic rail (y-axis shown in Fig. 1) after a 30 seconds’ relaxation at the levitation height to reduce the influence of the force relaxation on the subsequent measurement or calculation. The experimental and computed speed was chosen to be 1 mm/s in all cases and the maximum lateral displacement for traverse movement was 5 mm.

Figure 4.

Photos of the Y-Ba-Cu-O sample and magnetic rail demonstrators and associated schematic drawings of magnetic rail’s cross sections.

4.2. Experimental validation of the computed results

In the calculation, the magnetic rail were calculated by a 3-D analytical model in which the finite geometry of each magnet used in the magnetic rail is taken into account (Ma, et al., 2009). With the restriction of the available experimental tools, the material parameters involved in the E–J relation and the angular-dependence of the critical current density formulation can not be directly measured. In the following calculation, the necessary material parameters were determined according to the published literatures. The reported values of the pinning potential U₀ and flow resistivity ρ_f are very scattered. Here, pinning potential and flow resistivity were chosen to be 0.1 eV and 5×10^-10 Ωm respectively. Both of them are the frequently used values in the calculation Gou, et al.,2007a; Gou, et al.,2007b; Yoshida, et al., 1994). The anisotropic ratio for the melt-processed-single-domain Y-Ba-Cu-O has been measured to be ～3 at 77 K (Murakami, et al., 1991). Critical current density in the ab-plane Jab c, was determined by fitting one of the levitation force versus gap curve (the first curve in the following Fig. 7). In addition, Kim’s model (Kim, et al., 1962) was also employed to describe the field amplitude dependence characteristic of the critical current density in the calculation. All the parameters of the Y-Ba-Cu-O were summarized in Table 3.

4.2.1. Numerical results with different E-J constitutive relations and angular- dependence of critical current density formulas

Firstly, we evaluated the levitation force versus gap behavior with different angular-dependence of the critical current density formulations and E–J constitutive relations. In the 3-D finite-element model, the three different components of the magnetic force, i.e., vertical force perpendicular to the surface of the magnetic rail (z-axis), transverse force along the magnetic rail’s width (y-axis) and longitudinal force along the magnetic rail’s length (x-axis) can be calculated. From the computed results presented in Fig. 5, we can see that, there is no noticeable discrepancy among the computed levitation force for different test cases. This indicates that there is almost no difference between the power law model and flux flow and creep model when the pining potential U₀ and operating temperature are identical (the index n in power law model is ～ 15 at liquid nitrogen temperature with the pinning potential presented in Table 3). Moreover, the elliptical model is also viable to describe the anisotropic behavior of the high-Tc superconductor. Furthermore, as expected, the transverse force and longitudinal force shown in the inset of Fig. 5 are almost zero all the time with the variation of the gap because the applied field is symmetrical along the longitudinal and transverse directions of the magnetic rail.

In the following calculation, we choose the elliptical model for the purpose of introducing a new way to describe the anisotropic behavior of the critical current density and power law model to describe the E–J constitutive relation of the high-Tc superconductor because it seems to be more numerically stable to converge when compared to the flux flow and creep model.

Figure 5.

Comparison of levitation force versus gap curve between different E–J constitutive relations and angular-dependent of the critical current density formulations as well as the transverse and longitudinal force versus gap behavior with power law and elliptical model (inset). The applied field was provided by Rail_1 shown in Fig. 4.

Figure 6.

Comparison of the levitation force versus gap curve between the computed results and measured results under vertical movement in the applied field generated by Rail_1 shown in Fig. 4. The Y-Ba-Cu-O was bath-cooled at a position which was 60 mm above the rail where the applied field is weak and so this case can be approximately considered as zero field-cooled condition.

4.2.2. Comparison of the levitation force under vertical movement

The vertical force of the Y-Ba-Cu-O above the center of the magnetic rail1 with different field-cooled positions was calculated and compared with the measured data. From the compared results presented in Figs. 6-8, we can see that, the computed results agree well with the measured results in both hysteresis loop as well as detailed values at a certain gap.

For the levitation case shown in Figs. 6-7, we can find from both the computed and measured results that, on the downward branch of the levitation force versus gap curve, the levitation force and slope of the curve are continuously enhanced with the decrease of the gap. For the same gap, the levitation force of the downward branch is always larger than that of the upward branch, and this reveals an obvious hysteresis behavior of the levitation force. Furthermore, the levitation force versus gap behavior was calculated and measured another two times after the first one for the field-cooled test case at 25 mm above the magnetic rail1. It can be seen clearly that the computed results also compare well with the measured data for the second time and third time. This further confirms the validation of the 3-D finite-element method.

For the suspension test case shown in Fig. 8, the suspension force versus gap curve of both computed and measured cases indicates that, according to the stiffness defined in (Hull, 2000), the suspension force increases with the gap and its stiffness is always positive before the maximum absolute value of the suspension force is achieved at a certain gap. Therefore, the suspension system is always stable within this gap. The suspension force will decrease when the gap is further increased and the suspension system becomes unstable. Also, there is a clear hysteresis behavior of the suspension force by comparing the upward and downward branch of the curve. In addition, we also presented the computed results with flux flow and creep model, and the good agreement between the results of power law model and flux flow and creep model demonstrates that, the two E–J constitutive relations are also identical for the suspension case.

Figure 7.

Comparison of the levitation force versus gap curve between the computed results and measured results under vertical movement in the applied field generated by Rail_1 shown in Fig. 4. The Y-Ba-Cu-O was bath-cooled at a position of 25 mm above the rail which was a typical field-cooled condition. The levitation force versus gap curve was calculated and measured with continous three times here.

Figure 8.

Comparison of the suspension force versus gap curve between the computed results and measured results under vertical movement in the applied field generated by Rail_1 shown in Fig. 4. The Y-Ba-Cu-O is bath-cooled at a position of 8 mm above the rail.

Figure 9.

Comparison of the magnetic force as a function of lateral displacement between the computed and measured results under transverse movement in the applied field generated by Rail_2 shown in Fig. 4. The Y-Ba-Cu-O sample was field-cooled at a position of 30 mm, and the levitation gap where the magnetic force was calculated or measured was 18 mm above the rail.

4.2.3. Comparisons of the magnetic force under transverse movement

In order to verify the robustness of the performance of the 3-D finite-element model, the Rail_2 shown in Fig. 4, which is a Halbach array (Halbach, 1985), was employed to produce the applied field in this section.

During the transverse movement, both the vertical force and transverse force are a function of the lateral displacement. The compared results between the numerical results and measured results for different cases, i.e., levitation case with field-cooled above the levitation position, field-cooled at the levitation position, and suspension case with field-cooled below the suspension position, are presented in Figs. 9-11, respectively.

Both the numerical results and measured results shown in Figs. 9-11 indicate that, the absolute value of the transverse force increases with the lateral displacement and the

Figure 10.

Comparison of the magnetic force as a function of lateral displacement between the computed and measured results under transverse movement in the applied field generated by Rail_2 shown in Fig. 4. The Y-Ba-Cu-O sample was field-cooled at the identical height (18 mm) with the levitation height where the magnetic force was calculated or measured.

Figure 11.

Comparison of the magnetic force as a function of lateral displacement between the computed and measured results under transverse movement in the applied field generated by Rail_2 shown in Fig. 4. The Y-Ba-Cu-O sample is field-cooled at a position of 13 mm, which was below the levitation gap (18 mm) where the magnetic force was calculated or measured.

direction of the transverse force is opposite to the transverse movement, which indicates an inheret stable levitation can be acquired. Also, the magnetic force (i.e., the vertical force and the transverse force) are hysteretic as a function of the lateral displacement. Despite that the compared results of the transverse case is not as good as that of vertical case, as a whole, the numerical results are well comparable to the measured results in quality especially for the suspension case shown in Fig. 11. The fitting value of the critical current density in the ab-plane Jab c was derived from the case above the Rail_1 and the shift of the applied field are likely to be responsible for the discrepancy between the computed and measured results.

5. Optimization of the magnetic rail using the 3-D finite-element model

Magnetic rail is a key component to privide the applied magnetic field for the present magnet levitation system using bulk high-Tc superconductor, and the cost spend in building the magnetic rail occupies most part of the entire investment because the magnetic rail is required along the whole line. It is thereby meaningful to optimize the structure and geometric parameters of the magnetic rail in the purpose of getting a magnetic rail that holds the required levitation capability with a reduced cost. In this aspect, a lot of previous work has been reported and the dependence of the levitation force/guidance force on the parametes of the magnetic rail has been invesigated. However, all of these results were concluded from the numerical data calculated by a 2-D model that can just fit the experimental results in quality but fails in quantity with a reasonable value of critical current density (Song, et al., 2006; Dias, et al., 2010).

wSC (mm)	lSC (mm)	tSC (mm)	Jab c(A/m2)	Ec (V/m)	U0 (eV)	α
42	21	9	2.5×108	1×10-4	0.1	3

Table 4.

Parameters of high-Tc superconductor used in the calculation for obtaining an optimized magnetic rail.

Figure 12.

Two different magnetic rails with five permanent magnets derived from the Halbach array

Figure 13.

Chart of the the maglev system with three bulk Y-Ba-Cu-O undergoing the vertical or transverse movement above the rail

In this section, basing on the 3-D finite-element model introduced in this chapter, we calculate the levitation force and guidance force of three Y-Ba-Cu-O bulks samples above two differnt magnetic rails deriving from the traditional Halbach array. The geometric parameters of magnetic rail such as height, width are variable in the calculation, and then the depedence of levitation capability of the bulk Y-Ba-Cu-O on those parameters was studied. Compared with the previous work, the merit of the present calculation is that, the computed levitation force/guidance force is comparable to the real system with reasonable vaule of the critical current density and thus, the computed results can be used to conduct the practical design directly.

The geometric and material parameters of the Y-Ba-Cu-O sample were shown in Table 4 The speed of the samples in both vertical and transverse direction is 1 mm/s and the magnetization M₀ of the permanent magnet employed to assemble the magnetic rail is 8.9×10⁵A/m in all cases.

As for the structure of the magnetic rail, Halbach array is a better choice than the original type used in the high-Tc superconducting maglev train demonstrator because this structure can concentrate the magnetic field to its upper space where the high-Tc superconductors are placed, and thus can improve the utilization of the magnetic field (Jing, 2007). From the basic structure of the Halbach array shown in Fig. 12, we can derive two different types of magnetic rail, i.e., one has three permanent magnets magnetized in the horizontal direction and two permanent magnets magnetized in the vertical direction, the other has three permanent magnets magnetized in the vertical direction and two permanent magnets magnetized in the horizontal direction. These two magnetic rails derived from the Halbach array are presented in Fig. 12.

In the calculation, we ignore the possible interaction among the Y-Ba-Cu-O samples for simplicity, and calculate only two samples beacuse of the symmetry of the levitation system shown in Fig. 13. In the default case, the high-Tc superconductors were field-cooled at a position of 30 mm above the surface of the magnetic rail for the calculation of the levitation force, and for the calculation of the guidance force, the tranverse movement occurs at the same height as the field-cooled position, that was 12 mm above the surface of the magnetic rail. The main parameters that should be optimized in the two structure, i.e., Rail_A and Rail_B shown in Fig. 12, are the ratio between the width of two different magnetized magnets, the width of and the height of the magnetic rail. The following section shows the computed results of the optimization by varying those parameters.

5.1. Width ratio of the two different magnetized magnets

In this section, the levitation force and guidance force on the high-Tc superconductors with the variation of the width ratio, i.e., w_A1/w_A2 for Rail_A or w_B1/w_B2 for Rail_B, were calculated. In this calculation, the total width and the height of both rails were assumed respectively to be 130 mm and 30 mm and invariant with the change of the width ratio. For simplicity in drawing the following figures, the width ratio was replaced by an order and the corresponding relationship between the width ratio and the order was given in Table. 5.

Order	1	2	3	4	5	6	7	8	9	10	11
12	13	14	15	16	17	18	19	20	21
Width; ratio	0	0.01	0.1	0.16667	0.25	0.333333	0.4	0.5	0.66667	0.83	1
1.22	1.5	2	2.5	3	4	6	10	100	∞

Table 5.

Corresponding relationship between the width ratio and the order in optimazing the width ratio of the magnetic rails.

For the case of vertical movement, twenty-one different width ratios from zero to infinity were considered and the changing curve of the levitation force along with the width ratio at position of 5, 10, 15 and 20 mm above the Rail_A and the Rail_B was shown respectively in Fig. 14(a) and Fig. 14(b). Note that the extreme case with a width ratio of zero or infinity denotes that Halbach struture disappears and permanent magnets employed in the Rail only have one magnitized direction (vertical or horizontal direction). Both Fig. 14(a) and Fig. 14(b) display that, the levitation force, at all positions we presented, increases with the growth of the width ratio and reaches to a maximum value when the width ratio was ～0.83 for Rail_A and was ～1 for Rail_B and then drops with continuous growth of the width ratio. This finding indicates, the halbach array has a better performance because the levitation

Figure 14.

The changing curve of the levitation force on the high-Tc superconductors at position of 5, 10, 15 and 20 mm above the Rail_A (a) and Rail_B (b) with the growth of the width ratio. The width ratio is replaced by an order and the corresponding relationship between the width ratio and the order (abscissa) is given in Table. 5.

Figure 15.

The changing curve of the guidance force on high-Tc superconductors at a lateral displacement of 6 mm with the growth of the width ratio. The corresponding relationship between the width ratio and the order (abscissa) can be found in Table. 5. The total number of width ratio considered in this calculation was reduced and thus the order starts at 6 and ends at 16.

force at the extreme case is always the smallest one as has been verified by experiment (Jing, et al., 2007; Sotelo, et al., 2010), on the other hand, the optimized struture for Rail_A is that the width of magnet magnetized in horizontal direction is slightly larger than that magnetized in vertical direction, but for Rail_B, the two type magnets have an identical width ratio.

For the case of the transverse movement, only eleven different width ratios around 1 were considered because larger levitation force can be obtained in that area according to the above computed results of levitation force. The high-Tc superconductors have a field-cooled height of 12 mm and then were drived transversely at the same height with a maximum lateral displacement of 6 mm. Fig. 15 presents the guidance force with respect to the width ratio for the Rail_A and the Rail_B. We can find from this figure that, the guidance force exhibits a same tendency with the growth of the width ratio as that found in the case of levitation force, but the width ratio where the largest guidance force occurs is found to be 1.22～1.5 for Rail_A and is still ～1 for Rail_B.

As a whole, we suggest that, for the Halbach array as the magnetic rail used in the maglev train, the optimized width ratio between the vertical and transverse magnetized magnet is ～1.

5.2. Total width of the magnetic rail

In this calculation, the total width of both rails was a variable parameter with a fixed height of 30 mm and a fixed width ratio of 1.

Fig. 16 shows the changing curve of the levitation force on the high-Tc superconductors at positions of 5, 10, 15 and 20 mm above the Rail_A and Rail_B with the growth of the total width from 60 mm to 170 mm. These figures show that, with the increase of the position, the total width where the maximum value of the levitation force appears is shift from small to large value, e.g., the levitation force increases almost linearly with the total width at the position of 20 mm for both rails and the maximum value of the levitaion force even does not appear within our calculational range. This figure also indicates that, the increase of the total width does not always bring an ehencement of the levitation force especially for the case at a

Figure 16.

The changing curve of the levitation force on the high-Tc superconductors at positions of 5, 10, 15 and 20 mm above the Rail_A (a) and Rail_B (b) with the growth of the total width. The total width was varied from 60 to 170 mm.

Figure 17.

The changing curve of the guidance force on the high-Tc superconductors at a lateral placement of 6 mm with the growth of the total width. The total width was varied from 60 to 170 mm.

low position. Of course, if a larger levitaion capability is required at a high position, it is a viable way to lengthen the total width of the rail. For the typical position of 15 mm in the maglev train, it seems that, the optimized value of the total width appears around 130 mm for both rails because the slope of the curve begins to dcrease at this value.

Fig. 17 shows the changing curve of the guidance force on the high-Tc superconductors at a lateral displacement of 6 mm with the growth of the total width from 60 to 170 mm. From this figure, we can find that, the guidance force also increases firstly and then decreases with the growth of the total width, and the maximum value of the guidance force occurs around 120 mm for Rail_A and around 135 mm for Rail_B. Although the total width where the maximum value of the guidance force occurs is scattered above Rail_A and Rail_B, the variation of the guidance force above both rails in the range of 120～140 mm is not evident. It is therefore this range can be considered as the optimized range which is comparable to the total width of the high-Tc superconductors (126 mm) in determining the total width of the rail for a practical design.

5.3. Height of the magnetic rail

The height of the magnetic rail is another factor influencing the levitation capability of the high-Tc superconductor. Here, we fixed the width ratio to be 1 and total width to be 130 mm, and varied the height of the magnetic rail from 15 mm to 100 mm. The levitation forces above the Rail_A and Rail_B were plotted as a function of the height and the associate curves were given in Fig. 18. This figure displays clearly that, no matter the position where the high-Tc superconductors are placed, the levitation force increases drastically at first and then goes into gradually a saturated state with heightening the rail. This phenomenon indicates that, considering the cost of the rail, it is not reasonable to get a better levitation performance by increasing the height of the rail when the height is already high enough.

Figure 18.

The levitation force versus height of the Rail_A (a) and of the Rail_B (b). The height of both rails were varied from 15 to 100 mm.

6. Conclusion remarks

On the basis of Maxwell’s equations and Helmholtz’s theorem, the 3-D governing equations containing the anisotropic behavior of the HTS were deduced at length by introducing a current vector potential T and a simplified expression of the resistivity tensor for the high-Tc superconductor. The two common models to handle the highly nonlinear E–J characteristic, i.e., power law model and flux flow and creep model, were presented and, the critical current density J_c as a function of the angle φ between the orientation of the applied field and the c-axis of the high-Tc superconductor, was expressed by an elliptical model. The discrete matrix of the 3-D governing equations was derived by utilizing the Galerkin’s finite-element method and the Crank-Nicolson-θ method to numerically treat the governing equations one by one in space and time domain respectively. The corresponding nonlinear algebraic equations of the discrete matrix were effectively resolved by incorporating the Newton-Raphson method with an extended format of the Incomplete Cholesky-Conjugate Gradientmethod.

The computed results of the levitation force of a bulk high-Tc superconductor above a magnetic rail indicate that, the levitation force and associated hysteresis behavior will approach to a saturated state with the increment of the mesh density or the number of the time steps. This result is consistent with the finite-element theory, and also, supports the approach used to simplify the coupling problem between the differential and integral term presented in the governing equations.

Subsequentely, the 3-D finite-element model was validated by comparing the computed magnetic forces of a bulk Y-Ba-Cu-O sample with the measured data under vertical and transverse movements. Basing on the 3-D finite-element model, we also found that, the computed results of magnetic force of the high-Tc superconductor are almost identical for different angular-dependence of critical current density formulations (i.e., Yang’s model and elliptical model) and E–J constitutive relations (i.e., power law model and flux flow and creep model).

Lastly, using this 3-D finite-element model, we conduct a preliminary numerical work with aim of geting an optimized geometric parameter of the magnetic rails derived from the Halbach array. The results show that, when the width ratio of the two different magnetized magnets used in the magnetic rail is approximately identical and the total width of the rail is roughly indentical to the width of the high-Tc superconductors, the rail exhibits a better performance in considering both the levitation capability and the cost (be proportional to the voulme of the rail). Moreover, the levitation force has a saturated tendency with increasing continuely the height of the rail, so it is not reasonable to get a better levitation performance by increasing the height of the rail when the height is already high enough.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant 51007076 and 50777053, and the Innovation Foundation of Southwest Jiaotong University for Ph.D. Candidate. The authors would like to thank Dr. Frank N. Werfel of the Adelwitz Technologiezentrum GmbH (ATZ) for providing the Y-Ba-Cu-O sample employed in the experiment and to thank Dr. Hong-Hai Song, Dr. Zi-Gang Deng, Dr. Yi-Yun Lu, Dr. Min-Xian Liu, Dr. Jun Zheng, Dr. Fei Yen, Mr. Wei Liu, Mr. Qun-Xu Lin, Mr. Dong-Hui Jiang for their many fruitful discussions