Perturbation solutions for a micropolar fluid flow in a semi-infinite expanding or contracting pipe with large injection or suction through porous wall

1 Introduction

Since Eringen [1, 2] proposed a mathematical model to describe the non-Newtonian behaviour of liquids such as polymers, colloidal suspensions, animal blood and liquid crystals, there has been interest in micropolar fluids. In particular, micropolar fluids flowing in porous channels or pipes have received more attention due to their relevance to a number of practical biological problems. For example, Mekheimer and Elkot [3] presented a micropolar model for axisymmetric blood flow through an axially nonsymmetric but radially symmetric mild-stenosis-tapered artery. Mekheimer et al. [4] investigated the effects of an induced magnetic field on peristaltic transport of an incompressible micropolar fluid in a symmetric channel. Furthermore, Subhardra Ramachandran et al. [5] used Van Dyke's singular perturbation technique to study the heat transfer of a micropolar fluid past a curved surface with suction and injection. Anwar Kamal and Hussain [6] examined the steady, incompressible and laminar flow of micropolar fluids inside an infinite channel where the flow was driven due to a surface velocity proportional to the streamwise coordinates. Joneidi et al. [7] obtained similarity equations for a micropolar fluid in a porous channel and used the homotopy analysis method (HAM) to discuss the velocity distribution. In addition, Ariman et al. [8] and Lukaszewicz [9] gave reviews of micropolar fluid mechanics and its applications.

The purpose of this paper is to extend previous investigations by presenting analytical solutions for the flow inside a deforming porous pipe with large injection or suction. Equations describing the unsteady flow of an incompressible Newtonian fluid in a porous expanding channel are presented by White [10] as one of the new exact Navier-Stokes solutions attributed to Dauenhauer and Majdalani [11] . In their work [11], they numerically discussed the influence of the expansion ratio and Reynolds number on the velocity and pressure distribution. Furthermore, Majdalani, Zhou and Dawson [12] also obtained an asymptotic solution for the flow in a porous channel, with slowly expanding or contracting walls, by considering the permeation Reynolds number and expansion ratio as two small parameters. Boutros et al. [13, 14] also discussed the flow through an expanding porous channel or pipe using the Lie group method and obtained the analytical solution with the perturbation method. Recently Si et al. [15] also investigated the micropolar fluid in a porous deforming channel and discussed the effects of the micropolar parameter and the expansion ratio on the velocity and microrotation distribution. As further research, Li, Lin and Si [16] numerically analyzed the flow of a micropolar fluid through a porous pipe with an expanding or contracting wall.

In this paper, asymptotic solutions are constructed for the flow of a micropolar fluid through an expanding or contracting porous pipe. For large injection, analytical solutions are constructed using the Lighthill method, which eliminates singularity of the solution in the high order derivative [17-19];a series expansion matching method is used for large suction. The accuracy of the analytical solutions for each case is compared with its numerical results.

2 Preliminaries

Consider a micropolar fluid flowing through a pipe with a vertical moving porous wall. Here we assume that one end of the pipe is closed by a complicated solid membrane and the wall of the pipe moves in the radial direction and expands or contracts uniformly at a time-dependent rate a˙(t). In order to neglect the influence of the opening at the end, the length of the pipe is assumed to be semi-infinite [20] . Under the porous wall stipulation, the fluid is injected or aspirated uniformly and vertically through the pipe wall with an absolute velocity ν_w, which is proportional to the moving velocity at the wall surface. u = (u, v, 0) and ω = (0, 0, N) are the velocity vector and microrotation vector, respectively. N is the component of microrotation in the direction vertical to (r, z) plane, and u, v are the components of velocity in the direction of z and r, respectively. The flow configuration and the coordinate system are shown in Fig. 1.

Fig. 1

A model of a micropolar fluid through a porous expanding pipe.

Under these assumptions, the governing equations of the incompressible and homogeneous micropolar fluid flowing with no body force are expressed as follows:

where ρ and μ are the density and the dynamic viscosity, and j, y and κ are the micro-inertial coefficient, spin gradient viscosity and vortex viscosity, respectively. Here y is assumed to be

The corresponding boundary conditions are [11, 12, 15]

where A is the measure of the permeability. Here we also assume that there is a strong concentration of microelements, and the microelements close to the wall are unable to rotate [21, 22].

Introduce the stream function Ψ(r, z, t) such that

In this paper, the stream function Ψ and the microrotation velocity N are assumed as follows:

where η=(ra)2 and ν=μρ.

Similar to Dauenhauer and Majdalani [11], Uchida and Aoki [20], and Boutros et al. [13, 14], we substitute Eqs. (6) and (7) into governing equations and consider the similarity solutions with respect to space and time, then the following ordinary equations can be obtained:

where (f,g)=(FRe,GRe), Re=aνwν is the permeation Reynolds number, α=aa˙ν is the expansion ratio, and K=κμ and ζ=a2j are the micropolar parameters. In physical meaning, α is positive for expansion and negative for contraction.

The corresponding boundary conditions can be written as

3 Perturbation analysis for this problem

3.1 Solution for the large injection Reynolds number

For large injection Reynolds numbers, the asymptotic solution of (8) and (9), subject to the boundary conditions (10), is obtained by the Lighthill method. One treats ε=2Re as the perturbation parameter, the equations (8) and (9) then become

ε(1+K)(ηf‴+f″)+εα2(ηf″+f′)+εK4(ηg′+g)−f′2+ff″=λ,(11)

ε(1+K2)(η2g″+2ηg′)−εKζ2(ηg+2ηf″)+εα2(η2g′+2ηg)+12(fg+2ηfg′−2ηf′g)=0,(12)

where λ is an integral constant. Firstly, we introduce a variable transformation of η

η=ξ+εX1(ξ)+ε2X2(ξ)+O(ε3),(13)

where the functions X₁, X₂ are unknown and will be determined in the following process. One assumes that the functions f, g and the constant λ are expanded as

f(η)=∑i=0∞εifi(ξ),  g(η)=∑i=0∞εigi(ξ),  λ=∑i=0∞εiλi.(14)

Substituting (13)-(14) into (11)-(12) and collecting the same powers of ε, one can obtain the leading solution

f0f¨0−f˙02=λ0,  ξf˙0g0−ξf0g˙0−12f0g0=0,(15)

and the first order solution

f0f¨1−2f˙0f˙1+f¨0f1=−(1+K)(ξf¨0+f¨0)−α2(ξf¨0+f˙0)−K4(ξg˙0+g0)+λ1+2λ0X˙1,(16)

2ξf0g˙1−2ξf˙0g1+f0g1=−(2+K)(ξ2g¨0+2ξg˙0)+Kζ(ξg0+2ξf¨0)−α(ξ2g˙0+2ξg0)−f1g0−f0g0X˙1−2X1f0g˙0−2ξf1g˙0+2X1f˙0g0+2ξf˙1g0.(17)

Here ˙ denotes the derivative with respect to ξ.

3.1.1 A. the transformed boundary conditions at the wall of the pipe

We assume ξ˜ is the root of (13) at η = 1, then

ξ˜=1−εX1(ξ˜)−ε2X2(ξ˜)+O(ε3)=1−ε{X1(1)+X1˙(1)[−εX1(ξ˜)−ε2X2(ξ˜)]+⋯}−ε2{X2(1)+X1˙(1)[−εX1(ξ˜)−ε2X2(ξ˜)]+⋯}+⋯=1−εX1(1)−ε2[X2(1)−X1˙(1)X1(1)]+O(ε3),(18)

thus the conditions at the wall can be obtained

f|η=1=1⇒1=f|ξ=ξ˜=f|ξ=1+f˙|ξ=1{−εX1(1)−ε2[X2(1)−X˙1(1)X1(1)]+⋯}+⋯=f0|ξ=1+ε(f1−X1f˙0)|ξ=1+O(ε2),(19)

f′|η=1=0⇒0=f˙|ξ=ξ˜=f˙|ξ=1+f¨|ξ=1{−εX1(1)−ε2[X2(1)−X˙1(1)X1(1)]+⋯}+⋯=f˙0|ξ=1+ε(f˙1−X1f¨0)|ξ=1+O(ε2),(20)

g|η=1=0⇒0=g|ξ=ξ˜=g|ξ=1+g˙|ξ=1{−εX1(1)−ε2[X2(1)−X˙1(1)X1(1)]+⋯}+⋯=g0|ξ=1+ε(g1−X1g˙0)|ξ=1+O(ε2).(21)

Hence, the boundary conditions of f_i and g_i at η = 1 are

f0|ξ=1=1,  f1−X1f˙0|ξ=1=0,  ⋯,(22)

f˙0|ξ=1=0,  f˙1−X1f¨0|ξ=1=0,  ⋯,(23)

g0|ξ=1=0,  g1−X1g˙0|ξ=1=0,  ⋯.(24)

3.1.2 B. the transformed boundary conditions at the center of the pipe

One supposes that ξ^ is the root of (13) at η = 0, then

ξ^=−εX1(ξ^)−ε2X2(ξ^)+O(ε3)=−εX1(0)−ε2[X2(0)−X˙1(0)X1(0)]+O(ε3),(25)

thus we can induce

f|η=0=0⇒0=G|ξ=ξ^=f|ξ=0+G˙|ξ=0{−εX1(0)−ε2[X2(0)−X˙1(0)X1(0)]+⋯}+⋯=f0|ξ=0+ε(f1−X1f˙0)|ξ=0+O(ε2).(26)

Hence, the boundary conditions of f_i at η = 0 are

f0|ξ=0=0,  (f1−X1f˙0)|ξ=0=0,  ⋯.(27)

Using (22)-(24) and (27), the solution for (15) can be obtained

f0=sin(π2ξ),  g0=0,(28)

and then λ0=−π24 can be obtained. Substituting above results into (16), (17) yields the equations for f₁, g₁

sin(π2ξ)f¨1−πcos(π2ξ)f˙1−π24sin(π2ξ)f1=(1+K)[π38ξcos(π2ξ)+π24sin(π2ξ)]−α2[−π24ξsin(π2ξ)+π2cos(π2ξ)]−π22X˙1+λ1,(29)

and

ξsin(π2ξ)g˙1−π2ξcos(π2ξ)g1+12sin(π2ξ)g1+Kζπ24ξsin(π2ξ)=0.(30)

Here it should be noted that direct use of the method of variation of parameters will cause a singularity in the third-order derivative of f₁ [17-19]. In order to eliminate the singularity and to simplify the equation of f₁, we can set

(1+K)[π38ξcos(π2ξ)+π24sin(π2ξ)]−α2[−π24ξsin(π2ξ)+π2cos(π2ξ)]−π22X˙1+λ1=0.(31)

Then we have

X1(ξ)=1+K2ξsin(π2ξ)−α2πξcos(π2ξ)+2λ1π2ξ+C1^,(32)

where C1^ is a constant. Thus, the equation for f₁ becomes

sin(π2ξ)f¨1−πcos(π2ξ)f˙1−π24sin(π2ξ)f1=0.(33)

The solution of Eq. (33) is

f1(ξ)=C2^cos(π2ξ)+C3^[2πsin(π2ξ)−ξcos(π2ξ)],(34)

where C2^,  C3^ are still integral constants. According to the boundary conditions (22)-(24) and (27), C1^=0,  C2^=0,C3^=0, λ1=−(1+K)π24 can be determined. Thus, f₁, g₁ can be achieved as follows:

f1=0,  g1=Kζπ24ξ−12sin(π2ξ)∫ξ1t12csc(π2t)dt.(35)

Finally, one obtains the asymptotic solutions of f, g in terms of ξ

f(η)=sin(π2ξ),g(η)=εKζπ24ξ−12sin(π2ξ)∫ξ1t12csc(π2t)dt,(36)

where

η=ξ+ε[1+K2ξsin(π2ξ)−α2πξcos(π2ξ)−1+K2ξ].(37)

Fig. 2 shows the profiles of f′(η) and g(η) against η for the asymptotic and numerical results. Tables 1 and 2 give asymptotic and numerical values of f″(1) and g′(1) for some values of large injection Reynolds number and expansion ratio α, respectively. The results agree well.

Table 1

Values of f″(1) and g′(1) for large injection Reynolds number (α =-5, K = 0.2, ζ = 10).

	f″(1)		g′(1)
Re	Numerical	Asymptotic	Numerical	Asymptotic
50	-2.6641619556581	-2.7339624379749	-0.2128340193109	-0.2077811452861
100	-2.5700001017056	-2.5955567129755	-0.1036184958045	-0.1012267118060
150	-2.5366158007420	-2.5517506351567	-0.0681380895087	-0.0669125722108
200	-2.5195961015861	-2.5302623043972	-0.0507066428440	-0.0499726805118
500	-2.4884675811385	-2.4922614078153	-0.0199682862609	-0.0198384008062

Table 2

A comparison of f″(1) and g′(1) for different α (Re = 100, K = 0.3, ζ =10).

	f″(1)		g′(1)
α	Numerical	Asymptotic	Numerical	Asymptotic
-5	-2.569310686399770	-2.595556712975505	-0.155025023543228	-0.151840067709067
-2	-2.507901592906879	-2.517499337080236	-0.147379010820407	-0.149539460622566
2	-2.430977747070986	-2.418783550899264	-0.138027490934619	-0.146578283184495
5	-2.376702021744101	-2.348507888421025	-0.131584147099143	-0.144433235137893

Fig. 2

Variation of f′(η) and g(η) as Re = 100, α = -5, K = 0.2, ζ =10.

3.2 Solution for large suction Reynolds number

The boundary layer happens near the wall not only for the velocity but also for the microrotation when there is large suction. The solutions of (8) and (9), subject to the boundary conditions (10), can be obtained for large suction by using the method of matched asymptotic expansion. One treats ε=−2Re as the perturbation parameter, the equations (8) and (9) then become

ε(1+K)(ηf‴+f″)+εα2(ηf″+f′)+εK4(ηg′+g)+f′2−ff″=k,(38)

ε(2+K)(η2g″+2ηg′)−εKζ(ηg+2ηf″)+εα(η2g′+2ηg)−fg−2ηfg′+2ηf′g=0,(39)

where k is a constant of integration,

k=ε(1+K)σ+εαδ2+εKω4+δ2,(40)

and

δ=f′(0)=δ0+εδ1+ε2δ2+O(ε3),σ=f″(0)=σ0+εσ1+ε2σ2+O(ε3),ω=g(0)=ω0+εω1+ε2ω2+O(ε3),(41)

where the coefficients δ_i, σ_i and ω_i (i = 0,1,2, …) are constants determined by matching with the inner solution.

We assume that the forms of the outer solutions are

f=f0+εf1+ε2f2+O(ε3),   g=g0+εg1+ε2g2+O(ε3).(42)

Substituting (41)-(42) into (38)-(39), and equating the same power of the coefficient ε, one obtains

ε0:f′02−f0f″0=δ02,(43)

f0g0+2ηf0g′0−2ηf′0g0=0.(44)

ε1:  2f′0f′1−f″0f1−f0f″1=−(1+K)(ηf‴0+f″0)−α2(ηf″0+f′0)−K4(ηg′0+g0)+2δ0δ1+(1+K)σ0+αδ02+Kω04,(45)

2ηf′0g1−2ηf0g′1−f0g1=−(2+K)(η2g″0+2ηg′0)+Kζ(ηg0+2ηf″0)−2ηf′1g0−α(η2g′0+2ηg0)+2ηf1g′0+f1g0.(46)

ε2:  2f′0f′2−f″0f2−f0f″2=−(1+K)(ηf‴1+f″1)−α2(ηf″1+f′1)−K4(ηg′1+g1)+f1f″1−f′12+2δ0δ2+δ12+(1+K)σ1+αδ12+Kω14,(47)

2ηf′0g2−2ηf0g′2−f0g2=−(2+K)(η2g″1+2ηg′1)+Kζ(ηg1+2ηf″1)+f1g1−α(η2g′1+2ηg1)+f2g0+2ηf1g′1−2ηf′1g1−2ηf′2g0+2ηf2g′0.(48)

ε3:  2f′0f′3−f″0f3−f0f″3=−(1+K)(ηf‴2+f″2)−α2(ηf″2+f′2)−K4(ηg′2+g2)+f″1f2+f1f″2−2f′1f′2+2δ0δ3+2δ1δ2+(1+K)σ2+αδ22+Kω24,(49)

2ηf′0g3−2ηf0g′3−f0g3=−(2+K)(η2g″2+2ηg′2)+Kζ(ηg2+2ηf″2)−α(η2g′2+2ηg2)+f1g2+f2g1+f3g0+2ηf1g′2+2ηf2g′1−2ηf′1g2−2ηf′2g1−2ηf′3g0+2ηf3g′0.(50)

…….

The corresponding boundary conditions for the outer solutions are

fi(0)=0,  limη→0η12f″i=0, limη→0η12gi=0,  i=0,1,2,⋯.(51)

According to the boundary conditions (51), the solutions for (43) and (44) can be obtained

f0=δ0η,  g0=C0˜η12,(52)

where C0˜ is an integral constant. Because g₀ means the microrotation velocity of particles, its first-order derivative should be bounded, thus g′₀(0) → ∞ leads to C0˜=0. Then ω₀ = g₀(0) = 0, σ0=f″0(0)=0. Substituting (52) into (45), one obtains

2f′1−ηf″1=2δ1.(53)

Then the solution of (53) that satisfies the boundary conditions (51) is

f1=δ1η+C1˜η3,(54)

where C1˜ is an integral constant, which will be determined next. Substituting (52), (54) into (46), one obtains

g1−2ηg′1=0,(55)

whose solution is

g1=C2˜η12,(56)

where C2˜ is an integration constant. Similarly, g′₁(0) → ∞ leads to C2˜=0. Then

f1=δ1η+C1˜η3,  g1=0.(57)

Thus one obtains σ1=f″1(0)=0,   ω1=g1(0)=0. Substituting Eqs. (51) and (57) into (47), one obtains

f2=−6(1+K)C1˜δ0η2+αC1˜2δ0(3η3lnη−η3)+δ2η+3C1˜210δ0η5+C3˜η3,(58)

where C3˜ is an integral constant. Similarly, f‴2(0)→∞ leads to C1˜=0, thus

f1=δ1η,  f2=δ2η+C3˜η3.(59)

Then one obtains σ2=f″2(0)=0. Substituting (51), (57) and (59) into (48), one obtains

g2=C4˜η12.(60)

Similarly, g′2(0)→∞ leads to C4˜=0. Then one obtains ω₂ = g₂(0) = 0. Substituting (51), (57), (59) and (60) into (49), one obtains

f3=−6(1+K)C3˜δ0η2+αC3˜2δ0η3(3lnη−1)+δ3η+C5˜η3,(61)

where C5˜ is an integration constant. Similarly, f‴3(0)→∞ leads to C3˜=0, thus

f2=δ2η,  g2=0.(62)

In order to obtain the inner solution in the viscous layer, we introduce a stretching transformation τ = (1 - η)/ε. Substituting into Eqs. (11) and (12) yields

(1+K)[−f⃛+ε(τf⃛+f¨)]+α2[εf¨−ε2(τf¨+f˙)]+K4[−ε2g˙+ε3(τg˙+g)]+f˙2−ff¨=ε2δ2+ε3[(1+K)σ+αδ2+Kω4],(63)

(2+K)[ε2(τ2g¨+2τg˙)−ε(2τg¨+2g˙)+g¨]−Kζ[−ε3τg+ε2g−2ετf¨+2f¨]+ α[−ε3(τ2g˙+2τg)+ε2(2τg˙+2g)−εg˙]+2fg˙− 2f˙g+ε(2τf˙g−fg−2τfg˙)=0(64)

Here ˙ denotes the derivative with respect to τ. According to the boundary conditions (10), we assume the inner solutions near the wall to be

f(τ)=1+∑i=1∞εiϕi(τ),   g(τ)=∑i=1∞εiψi(τ).(65)

Substituting (65) into (63) and (64) yields the following equations

ε1:(1+K)ϕ⃛1+ϕ¨1=0,(66)

(2+K)ψ¨1−2Kζϕ¨1+2ψ˙1=0.(67)

ε2:(1+K)(ϕ⃛2−τϕ⃛1−ϕ¨1)−α2ϕ¨1−ϕ˙12+ϕ¨2+ϕ1ϕ¨1=−δ02,(68)

(2+K)(ψ¨2−τψ¨1−2ψ˙1)−2Kζϕ¨2−αψ˙1+2ψ˙2+ 2ϕ1ψ˙1−2ϕ˙1ψ1−ψ1=0.(69)

…….

The boundary conditions corresponding to the inner solution are

ϕi(0)=0,  ϕ˙i(0)=0,  ψi(0)=0,  i=1,2,⋯.(70)

The solution of (66) satisfying the boundary conditions (70) is

ϕ1=D1(e−A0τ+A0τ−1),(71)

where A0=11+K, and D₁ is an integral constant. Then the first two terms of the inner solution of f can be expressed as

f(τ)=1+D1(e−A0τ+A0τ−1)ε.(72)

The outer solution of f, expressed in terms of the inner variable τ, is

f(η)=δ0+(δ1−δ0τ)ε+(δ2−δ1τ)ε2+⋯.(73)

As τ → ∞, matching the inner solution (72) with (73) gives

δ0=1,  δ1=1A0,  D1=−1A0.(74)

The solution of (67) satisfying the boundary conditions (70) is

ψ1=KB0ζB0−A0(e−A0τ−1)+D2(e−B0τ−1),(75)

where B0=22+K, and D₂ is an integral constant. Similarly, as τ → ∞ in (75), D₂ can be determined:

D2=−KB0ζB0−A0.(76)

The solution of (68) satisfying the boundary conditions (70) is

ϕ2=−(α2+2A0)τe−A0τ−(αA0+4A02−D3A0)e−A0τ+(D3−α2−2A0)τ+4A02+αA0−D3A0,(77)

where D₃ is an integral constant. Then the inner solution of f can be expressed as

f=1−εA0(e−A0τ+A0τ−1)+ε2[−(α2+2A0)τe−A0τ−(αA0+4A02−D3A0)e−A0τ+(D3−α2−2A0)τ+4A02+αA0−D3A0]+⋯.(78)

As τ → ∞, matching the inner solution with outer solution, one obtains

D3=1A0+α2,   δ2=3A02+α2A0.(79)

similarly, one can obtain

ψ2=ζK(2+K)(4+3K)[(48K+108K2+78K3+18K4+ 8Kα+10K2α+3K3α)τe−A0τ−(24K+46K2+25K3+ 3K4+8Kα+14K2α+6K3α)τe−B0τ+(8Kα+18K2α+ 13K3α+3K4α−48−156K−186K2−96K3− 18K4)e−A0τ+(48+156K+182K2+84K3+6K4− 4K5−8Kα−18K2α−13K3α−3K4α)e−B0τ+ (4K2+12K3+12K4+4K5)e−(A0+B0)τ].(80)

Hence, the complete solutions of (8) and (9) satisfying the boundary conditions (10) for large suction can be obtained as follows:

f=η+ε{(1+K)η−eη−1ε(1+K)[α2+3(1+K)−(2+2K+α2)η]}+ε2{3(1+K)2+α(1+K)2−eη−1ε(1+K)[α(1+K)2+3(1+K)2]},(81)

and

g=εζ{eη−1ε(1+K)[8(1+K)+α−(α+6K+6)η]+ 1+K2+Ke2(η−1)ε(2+K)[−10−3K−2α+(6+K+2α)η]}+ ε2ζ{eη−1ε(1+K)(1+K)(Kα−6−6K)K+ e2(η−1)ε(2+K)[6+12K+2K2K+4K(7+7K−K3)(2+K)(4+3K)−(1+K)α]+ 4K(1+K)3(2+K)(4+3K)e(4+3K)(η−1)ε(1+K)(2+K)}.(82)

Fig. 3 shows the profiles of f′(η) and g(η) against η for the asymptotic and numerical results. Tables 3 and 4 give asymptotic and numerical values of f″(1) and g′(1) for some values of large suction Reynolds number and expansion ratio, respectively.

Fig. 3

variation of f′(η) and g(η) as Re = -100, α - 5, K = 0.2, ζ = 10.

Table 3

Values of f″(1) and g′(1) for large suction Reynolds number (α = -5, K = 0.2, ζ =10)

	f″(1)		g′(1)
Re	Numerical	Asymptotic	Numerical	Asymptotic
-60	-26.029180587850	-26.083333333333	-1.9697087541999	-1.9620553359684
-90	-38.553222227512	-38.583333333333	-1.9168665563568	-1.9140974967062
-100	-42.723767379927	-42.750000000000	-1.9066699865053	-1.9045059288538
-125	-53.146883680302	-53.166666666667	-1.8885407560495	-1.8872411067194
-150	-63.567422760150	-63.583333333333	-1.8766007832311	-1.8757312252964

Table 4

The comparison of f″(1) and g′(1) for different α (Re = -100, K = 0.2, ζ =10)

	f″(1)		g′(1)
α	Numerical	Asymptotic	Numerical	Asymptotic
-5	-42.7237673799268	-42.7500000000000	-1.90666998651	-1.90450592885
-2	-41.4054335101753	-41.5000000000000	-1.91280714642	-1.90450592885
2	-39.6395810279868	-39.8333333333333	-1.92210393045	-1.90450592885
5	-38.3082775719552	-38.5833333333333	-1.93007245250	-1.90450592885