Convenient Debonding Strength Evaluation for Spray Coating Based on Intensity of Singular Stress

Abstract

Hearth rolls are used in continuous annealing furnace to produce thin steel sheet. The roll surface is usually coated by using thermal spraying, which has high adhesive strength and wears resistance. However, in the ceramics coating, thermal stress caused during heating and cooling process in the furnace may lead to debonding due to the low toughness of ceramics. In order to improve the heat resistance of the thermally sprayed coating, it is essential to evaluate the debonding strength. Generally, heat resistance of thermal spray coating is evaluated by thermal shock test prescribed by JIS H8304 although few research is available in terms of singular stress at the end of the interface for JIS specimen under thermal shock. This paper focuses on the intensity of the singular stress at the end of interface. Then, the most suitable conditions are discussed with varying the coating material and the coating thickness.

1. Introduction

In continuous annealing furnace for producing steel sheet, hearth rolls as shown in Fig. 1(a) are used to produce thin sheet. To prevent against steel adhesive and to improve wear resistance, ceramics spraying coating is usually applied to the roll surface. Due to its excellent chemical stability and high hardness at high temperatures, the wide use of ceramics spraying coating is expected. However, the peeling of the coating layer should be considered after long-term use, which is mainly caused by heating and cooling process in the furnace. Therefore, the ceramic coating is not usually used as a single layer, while multi-layer coating with adhesive layer is usually employed.

Fig. 1.

Specimen for thermal shock test and the test conditions. (a) Roll used in annealing furnace (b) JIS 8304 specimen for thermal shock test 50×50×(H+h+l) (c) Thermal cycle of thermal shock testing. (Online version in color.)

In order to improve the thermal shock resistance of the ceramic coating, it is necessary to evaluate the peel strength of the sprayed coating accurately. The thermal shock resistance of ceramic coating is usually evaluated by the thermal shock test as shown in Figs. 1(b), 1(c), which is specified by JIS8304 (2007).¹⁾ Mutoh et al have reported the thermal shock damage characteristics based on experimental results and finite element analysis.²⁾ In order to evaluate the interfacial strength of dissimilar adhesive structures more accurately, it is necessary to consider the singularity of thermal stress at the corner of the interface. So far, few studies have been conducted to the research of intensity of singular stress on this multi-layer structure under thermal shock. Therefore, in this study, the yttrium (Y₂O₃) stabilized zirconia (ZrO₂) sprayed coating shown in Table 1 is considered mainly focusing on the intensity of singular stress at the edge of coating layer under the thermal shock. Then, the effect of the coating thickness on the peeling strength of the sprayed coating is investigated. In our previous study, a useful method for analyzing the intensity of singular stress has been proposed for an adhesive bonding plate under bending and tension.^3,4) Thus in this paper, this method will be applied to analyzing ceramics coating under thermal shock.

Table 1. Material Properties depending on temperature.

Material		Young’s modulus E [GPa]	Poisson’s ratio ν	Thermal expansion αe[10−6/K]	Dundurs parameter α	Dundurs parameter β	Singular index λ
Top Coating	20YZ (ZrO2-20wt%Y2O3)	20	0.25	7.5–9.7	0.969	0.199	0.8399– 0.8899 Fig. 4(c)
	12YZ (ZrO2-12wt%Y2O3)			8.3–9.93
	8YZ (ZrO2-8wt%Y2O3)			8.7–10.05
Bonded Coating (CoNiCrAlY)		103–180 Fig. 4(a)	0.33	2.8–16.8 Fig. 4(b)
					0.978	0.188	0.9987– 0.9988 Fig. 4(c)
Substrate (SUS304)		85–200 Fig. 4(a)	0.3	16.8–19.2 Fig. 4(b)

2. Thermal Shock Test for Evaluating the Peeling Strength of Sprayed Coating

JIS 8304 specimen is cubic shaped with length×width×thickness=50 mm×50 mm×10 mm. The substrate SUS304 and the top coating yttrium Y₂O₃with 8wt% – 20wt% of ZrO₂ (In the following % will be omitted) are bonded with CoNiCrAlY (See Fig. 1). In the experiments, the bond coating thickness is fixed as 100 μm, and the bond coating surface roughness is adjusted to Ra=4.5 μm to ensure the adhesive strength with the top coating. After the top coating is sprayed to the bond coating surface, the top coating is polished to the thickness of 150 μm and surface roughness of Ra=1.0 μm. Thermal shock test was carried out by using an electric heating furnace by adjusting the ambient temperature. As shown in Fig. 1(c), (1) the ambient temperature was first raised from 20°C to 1000°C during the first 3000 seconds, (2) and then during the following 1800 seconds, the temperature was kept as 1000°C to ensure that the specimen is uniformly heated, (3) after that, the specimen was taken out from the furnace and immersed into cold water of 25°C. Repeat this process until the delamination is observed. Then, number of cycle Nc is recorded by visual inspection when cracking, peeling or blistering of specimen occurs. The thermal shock resistance is evaluated by numbers of cycle Nc.

Table 1 shows material properties depending on temperature, that is, Young’s modulus, Poisson’s ratio, thermal expansion coefficient, Dundurs parameter defined later by Eq. (1) and singular index.^5,6) Three kinds of top coatings having different content of yttrium Y₂O₃ are considered. Notations 20YZ, 12YZ, 8YZ stand for 20%, 12%, 8% yttrium content, respectively. As shown in Fig. 2, with increasing yttrium content number of cycle until debonding Nc decreases. In Fig. 2 the y-axis Nc normalized by the standard number of cycle No.

Fig. 2.

Results of thermal shock test.

3. Analysis Method for the Singularity of Thermal Stress

In this study, finite element method analysis is applied to evaluating the singular stress at the end of interface assuming the delamination appearing at this point. A two-dimensional models shown in Fig. 3 are use to simulate the specimen. In Eq. (1) the subscript 1, 2 and 3 represent the top coating, bond coating and substrate respectively, whose material properties are shown in Table 1. Thermo-structural elastic analysis is performed by employing FEM code MARC MENTAT (MSC Software Corporation, 2012) with the 4 nodes quad element. Multi-frontal method is used in the solution of simultaneous equations. Notations α and β are known as Dundurs’ parameters defined by the following equations,^5,6) here v is Possion’s ratio and G = E / 2(1+ v) is shear modulus.   

$$\begin{array}{l}\alpha =\frac{G_1\left({\kappa }_2+1\right)-G_2\left({\kappa }_1+1\right)}{G_1\left({\kappa }_2+1\right)+G_2\left({\kappa }_1+1\right)},\\ \beta =\frac{G_1\left({\kappa }_2-1\right)-G_2\left({\kappa }_1-1\right)}{G_1\left({\kappa }_2+1\right)+G_2\left({\kappa }_1+1\right)},\\ {\kappa }_j=\{\begin{array}{c}\frac{3-{\nu }_j}{1+{\nu }_j}\left(\text{plane}\hspace{0.17em}\text{stress}\right)\\ 3-4{\nu }_j\left(\text{plane}\hspace{0.17em}\text{strain}\right)\end{array}\hspace{1em}\left(j=1,2\right)\end{array}$$(1)

Fig. 3.

Analytical model and FEM mesh. (a) Overview (b) End of interface (c) FEM mesh. (Online version in color.)

Furthermore, the singular index λ, that is, the exponent of r as shown in Fig. 3(b), can be obtained by solving the following equation. Here r is the distance from the end of interface in a dissimilar materials bonded plate.   

$$\begin{array}{l}\left[si{n}^2\left(\frac{\pi }{2}\lambda \right)-{\lambda }^2\right]{\beta }^2+2{\lambda }^2\left[si{n}^2\left(\frac{\pi }{2}\lambda \right)-{\lambda }^2\right]\alpha \beta \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}+{\lambda }^2\left({\lambda }^2-1\right){\alpha }^2+\frac{si{n}^2\left(\pi \lambda \right)}{4}=0\end{array}$$(2)

Figures 4(a) and 4(b) show the material properties of three layers: Young’s modulus E and coefficient of expansion α_e. It’s found that all properties of these three materials vary with the temperature except the Young’s modulus of top coating. Therefore, the singular stress field is also changed by the variety of material properties.

Fig. 4.

Relationship between temperature. (a) Relationship between temperature and modulus of longitudinal elasticity (b) Relationship between temperature and the coefficient of liner expansion (c) Relationship between temperature and singular index λ. (Online version in color.)

Figure 4(c) shows singular index λ₁ and λ₂ as functions of temperature. Singular index λ₂ of the interface between the bond coating and substrate is close to 1, which means the singularity is very small. Therefore, this study will mainly focus on the singular stress with λ₁ between the top coating and bond coating as shown in Fig. 3(a). During the thermal cycles as shown in Fig. 1, the water-cooled condition is simulated by giving temperature 20°C to the entire surface of the specimen at t=4800 s. Then the temperature distribution is obtained by thermal analysis and after that thermal stress is determined by thermo-elastic analysis.

The thermal stress σ_y is considered under the temperature change T→T+ΔT in the bonded structure as shown in Fig. 5(a). It should be noted that the non-singular term σ₀ should be eliminated.⁷⁾ It is known that the remaining singular term (σ_y−σ₀) has a singularity of r^1−λ expressed as Eq. (3). The intensity of singular stress field caused by the thermal stress is equivalent to the one subjected to the tension of σ₀ determined by Eq. (4).⁸⁾   

$${\sigma }_y-{\sigma }_0=\underset{r\to 0}{lim}\frac{K_{\sigma }}{r^{1-\lambda }},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{K}_{\sigma }=-{\sigma }_0{F}_{\sigma }{W}^{1-\lambda }$$(3)

  

$${\sigma }_0=-\mathrm{\Delta }\alpha \mathrm{\Delta }E\mathrm{\Delta }T,$$(4)

where,   

$$\mathrm{\Delta }\alpha =\{\begin{array}{c}{\alpha }_1-{\alpha }_2 \left(\text{plane}\mathrm{\hspace{0.17em} }\text{stress}\right)\\ \left(1+{\nu }_1\right){\alpha }_1-\left(1+{\nu }_2\right){\alpha }_2\left(\text{plane}\mathrm{\hspace{0.17em} }\text{strain}\right)\end{array}$$

  

$$\mathrm{\Delta }E=\frac{8}{\{\frac{\left({\kappa }_1-3\right)}{G_1}-\frac{\left({\kappa }_2-3\right)}{G_2}\}}$$

  

$${\kappa }_j=\{\begin{array}{c}\frac{3-{\nu }_j}{1+{\nu }_j}\left(\text{plane}\mathrm{\hspace{0.17em} }\text{stress}\right)\\ 3-4{\nu }_j\left(\text{plane}\mathrm{\hspace{0.17em} }\text{strain}\right)\end{array}$$

  

$$G_i=\frac{E_i}{2\left(1+{\nu }_i\right)}(\begin{array}{c}i=1:\text{Top}\mathrm{\hspace{0.17em} }\text{coat},\\ i=2:\text{Middle}\mathrm{\hspace{0.17em} }\text{coat}\end{array}$$

Fig. 5.

Comparing known problem and unknown problem. (a) Thermal singular stress field for ceramic coating as a unknown problem for the analytical (b) The bonded strip model as the known problem whose the intensity of singular stress have been analyzed. (Online version in color.)

The bonded plate subjected to tension shown in Fig. 5 has been accurately analyzed by the body force method.^3,4) Since the singular stress fields are similar if the local geometries around the singular points are same, the stress intensity factor of unknown problem shown in Fig. 5(a) can be determined by apply this known solution in Fig. 5(b) to Eq. (5).^3,4,9)   

$$\frac{K_{\sigma }^{}}{K_{\sigma }^{*}}=\frac{F_{\sigma }{\sigma }_0^{}{W}^{1-\lambda }}{F_{\sigma }{\sigma }_0^{*}{W}^{1-\lambda }}=\frac{{\sigma }_y^{FEM}-{\sigma }_0^{}}{{\sigma }_y^{FEM*}}$$(5)

Here the superscript* means known reference problem. Table 2 shows FEM results ${\sigma }_y^{FEM}$ with minimum element size 1/3⁸ mm and 1/3¹⁰ mm. The stress distributions along the interface in Figs. 5(a), 5(b) at the time t=4800 s are indicated. The value of ${\sigma }_y^{FEM}$ in Table 2 are depending on the mesh size. However, the stress ratio $\left({\sigma }_y-{\sigma }_0\right)/{\sigma }_y^{FEM*}$ is independent of the mesh size and does not depend on the distance r from the end of interface. This is because FEM error due to the finiteness of the mesh can be eliminated by applying the same mesh pattern to the unknown problem and reference problem.^3,4) Since the reference problem has already been solved in the previous study,^3,4) the intensity of singular stress for problem shown in Fig. 5(a) can be obtained by using Eq. (5).

Table 2. Stress distribution on the interface and the ratio of the stress.

Mesh size=1/38 mm				Mesh size=1/310 mm
r/W	${\sigma }_y^{FEM*}$ [MPa]	${\sigma }_y-{\sigma }_0$ [MPa]	$\frac{{\sigma }_y^{FEM}-{\sigma }_0}{{\sigma }_y^{FEM*}}$	r/W	${\sigma }_y^{FEM*}$ [MPa]	${\sigma }_y-{\sigma }_0$ [MPa]	$\frac{{\sigma }_y^{FEM}-{\sigma }_0}{{\sigma }_y^{FEM*}}$
0/6 561	2300.0	1389.8	0.604	0/59 049	2946.4	1778.4	0.604
1/6 561	1741.1	1052.3	0.604	1/59 049	2230.4	1346.2	0.604
2/6 561	1639.3	990.9	0.604	2/59 049	2099.9	1267.4	0.604
3/6 561	1576.8	953.2	0.604	3/59 049	2019.8	1219.1	0.604
4/6 561	1535.7	928.2	0.604	4/59 049	1967.1	1187.3	0.604
5/6 561	1503.4	908.6	0.604	5/59 049	1925.8	1162.4	0.604

4. Results and Discussion for the Intensity of Singular Stress Field

4.1. Occurrence Time of the Largest Singular Stress

In the first place, the thermal load as shown in Fig. 1(c) is applied to the standard specimen with coating thicknesses H=0.15 mm and h=0.1 mm (See the model in the lower left corner of Fig. 6(b)). Although residual stress may appear in the ceramic spraying coating, but a number of micro-cracks are observed in the coating studied in this paper. Therefore, the residual stress is not considered because most of the residual stress is released after the crack formation. The thermal stress is not very large during the heating process in the thermal cycling, thus we mainly consider the thermal stress during the cooling process. Figure 6(a) shows the time variation of the interface stress ${\sigma }_y^{FEM}$ obtained by FEM around the start of cooling t=4800 s. The results for three types of top coatings are indicated. As shown in Fig. 6, the maximum value of ${\sigma }_y^{FEM}$ always occurs immediately after the start of cooling. The increment percentage of 8YZ is the smallest among all. Figure 6(b) shows the time variation of the intensity of singular stress field for 8YZ around the start of cooling. Different from ${\sigma }_y^{FEM}$ in Fig. 6(a), the intensity K_σ in Fig. 6(b) takes a maximum value at the constant temperature of 1000°C and decreases immediately after the start of cooling. It should be noted that singular index varies depending on temperature, but from Fig. 4(c), singular index is almost constant as λ≅0.88 during the first moment of cooling. Therefore, the value of K_σ = −F_σσ₀W^1-λwith λ≅0.88 can be compared from Fig. 6(b).

Fig. 6.

Stress and intensity of the singular stress at the end of interface. (a) Stress at the end of interface by FEM (mesh size 1/3¹⁰ mm) (b) Intensity of the singular stress at the interface for top coating 8YZ. (Online version in color.)

To figure out the reason that the maximum value of ${\sigma }_y^{FEM}$ occurs immediately after the start of cooling, while K_σ decreases at the same time monotonically, all variables in Eq. (3) are demonstrated in Table 3 and Fig. 7, so that a comprehensive consideration can be taken. The range of time is from the start of cooling t=4800 s to 4800.12 s. From Table 3 and Fig. 7, it is found that F_σ and W^1-λ are almost constant within the variation of 1.9% and 5.2% respectively, while −σ₀ decreases 18.7% over time. For this reason, the tensile interfacial stress ${\sigma }_y^{FEM}$ increases at the early stage of cooling while the intensity of singular stress field decreases by 16.5% as is shown in Fig. 7. Therefore, in the following analysis of K_σ, we mainly pay attention to the results of temperature 1000°C.

Table 3. This table shows that K_σ value of the influence factor changes with time for top coating 8YZ.

T [s]	4800	4800.01	4800.02	4800.03	4800.04	4800.06	4800.08	4800.12	Rate of change
T [°C]	1000	978.4	966.8	957.1	948.4	933.5	921.1	901.0
Kσ[MPa·m0.1127]	421.2	397.0	391.3	386.1	381.0	371.4	363.2	351.6	16.5%↓
Fσ	0.423	0.425	0.427	0.427	0.427	0.425	0.423	0.419	1.9%↓
−σ0 [MPa]	629.6	592.2	578.8	568.0	558.7	543.2	530.8	511.8	18.7%↓
W1−λ	1.554	1.576	1.584	1.592	1.598	1.608	1.617	1.640	5.2%↑
${\sigma }_y^{FEM}$ [MPa]	761.2	791.9	806.7	815.9	821.7	827.2	828.4	825.5	7.8%↑

Fig. 7.

K_σ value of the influence factor changes with time.

4.2. Intensity of Singular Stress for 2-layer Coating

Here, three materials of the top coating with different content of yttrium (20YZ, 12YZ, 8YZ) are analyzed to determine the maximum intensity of singular stress field under thermal shock test. Figure 8 shows the intensities of singular stress for different materials with varying the top coating thickness H while keeping the bond coating layer thickness h=0.1 mm. From Fig. 8, it is found that the material 8YZ with the lowest content of Y₂O₃ exhibits the best thermal shock resistance under any thickness, which means thermal shock resistance has a negative correlation to the amount of added yttrium. It is also found that thinner top coating has higher thermal shock resistance because of smaller value of K_σ.

Fig. 8.

For 20YZ, 12YZ, 8YZ at 1000°C with varying ceramic coat thickness H when h=0.1 mm. (Online version in color.)

Table 4 shows the results of K_σ when both thicknesses H and h are changed. From Table 4 it is found that the singular stress intensity reaches the minimum value when the thickness ratio H/h≅2. In engineering application the most suitable ratio H/h≅2 is often used in the coating thickness design empirically. The present study has verified the effectiveness of H/h≅2 in terms of the intensity of singular stress.

Table 4. K_σ for 8YZ at 1000°C of two layers model. (Red figure shows minimum value when H=const.) [MPa·m^0.1127] (Online version in color.)

4.3. Intensity of Singular Stress for 3-layer Coating

In the previous section, the intensity of singular stress K_σ in a two-layer coating with one layer of top coating and one layer of bond coating was discussed. Then, it was found that yttrium-rich 20YZ coating has the weakest thermal shock resistance because of the largest stress singularity. While, in a different perspective, it is known that yttrium-rich coating is useful for preventing against steel plate adhesive. Therefore, to improve thermal shock resistance for high yttrium coating is expected. To accomplish this, two layers of top coating whose thickness is (H1+H2) as shown in Fig. 9 is employed. The outer layer of top coating 20YZ has higher adhesive preventing ability and higher roughing resistance. The second layer coating 8YZ has higher peeling resistance with the bond coat.

Fig. 9.

Three layers ceramic coating model used to analysis. (Online version in color.)

Table 5 shows the results of K_σ when the thickness of the first layer of top coating is fixed as H₁=0.075, the second layer of top coating H₂ varies from 0.015 mm to 0.225 mm and bond coating h varies from 0.005 mm to 0.15 mm. Similarly to the two-layer coating, when the ratio of top coating thickness to the bond coating thickness (H₁+H₂)/h≅2, the intensity of singular stress K_σ reaches minimum value. In addition, with decreasing thickness (H₁+H₂), K_σ decreases and thermal shock resistance increases.

Table 5. K_σ at 1000°C of three layers model when fixed H₁=0.075 mm. [MPa·m^0.1127] (Online version in color.)

Table 6 shows the results of K_σ when the thickness of the second layer of top coating is fixed as H₂=0.075, and the first layer of top coating H₂ varies from 0.015 mm to 0.225 mm and bond coating h varies from 0.005 mm to 0.15 mm. As shown in Table 6, when the ratio of top coating thickness and the bond coating thickness (H₁+H₂)/h≅2, the intensity of singular stress K_σ reaches minimum value which is similar as Tables 4 and 5. Figure 10 shows the distribution of the intensity of singular stress K_σ when the thickness of the top coating and the thickness of bond coating are fixed as H₁+H₂=0.15 mm and h=0.075 mm. From Fig. 10, it is found that the intensity of singular stress increases with increasing the top coating thickness H₁.

Table 6. K_σ at 1000°C of three layers model when fixed H₂=0.075 mm. [MPa·m^0.1127] (Online version in color.)

Fig. 10.

K_σ at 1000°C of three layers model when fixed (H₁+H₂)=0.15 mm, h=0.075 mm. (Online version in color.)

Fig. 11.

Comparison of the results for the three-layer and two-layer. (Online version in color.)

Next, the strength improvement of the 3-layer coating is considered in contrast with the two layers coating. Assume the first layer is 20YZ with excellent higher adhesive preventing ability with the bond coat thickness h=0.075 mm. The minimum value of singular stress field intensity K_σ=424 MPa·m^0.1127 is obtained when H₁=H₂=0.075 mm. Here, H₁=H₂=0.075 mm is the minimum allowable thickness for thermal spraying coating. On the other hand, the minimum value of K_σ=461 MPa·m^0.1127 for a two-layer coating when H=0.15 mm and h=0.075 mm is about 8% higher than that of a three-layer coating. As shown in this example, when the top coating thickness H(=H₁+H₂) and bond coating thickness h are fixed, it is found that the three-layer coating is better that two-layer coating in strength because of smaller K_σ. Figure 12 shows the relationship between K_σ and dimensionless cycle number N_c/N_o. Here, vertical axis is dimensionless in the same manner as Fig. 2. It is seen that the number of cycles until peeling increases about three times by using three-layer coating.

Fig. 12.

Relationship between N_c and K_σ.

5. Conclusion

This study focuses on the improvement of peeling resistance for the thermal sprayed coating roll using in the continuous annealing furnace. Considering the two-dimensional model based on the thermal shock test specimen, the effect of material properties and thickness of the sprayed coating on the intensity of singular stress field are discussed by using FEM analysis. As a result the following conclusions can be obtained.

(1) In the two-layer coating consisting of top coating and bond coating, experimental results shows that higher thermal shock peeling resistance can be obtained by smaller content of Y₂O₃. This can be explained from the viewpoint of singular stress intensity K_σ near the end of interface for the top coating and bond coating. In other words, with decreasing Y₂O₃ content, K_σ decreases and therefore thermal shock resistance increases.

(2) Compared with the 2-layer structure, the 3-layer structure has higher thermal shock resistance because the stress intensity factor is 8% less than the 2-layer structure without losing a good roughening resistance. The first layer 20YZ (higher content of Y₂O₃) has excellent roughening resistance while the second layer 8YZ (lower content of Y₂O₃) can keep the thermal shock resistance.

(3) The optimal thickness ratio of top coating bond coating H/h is discussed from the view of intensity of singular stress near the end of interface. The theoretical analysis and engineering application all indicated that the optimal thickness ratio H/h≅2.0 both in 2-layer structure and 3-layer structure.

References

• 1)  JIS H 8304 2007: 2007, Ceramics Sprayed Coatings.
• 2)   Y.  Mutoh,  I.  Sakamoto,  O.  Waranabe and  T.  Nishimura: Trans. Jpn. Weld. Soc., 7 (1989), No. 2, 42.
• 3)   Y.  Zhang,  N. A.  Noda,  K.  Takaishi and  X.  Lan: Trans. Jpn. Soc. Mech. Eng. A, 77 (2011), No. 774, 128.
• 4)   Y.  Zhang,  K.  Takaishi,  N. A.  Noda and  X.  Lan: Trans. Jpn. Soc. Mech. Eng. A, 77 (2011), No. 784, 2076.
• 5)   D. B.  Bogy: Trans. ASME, J. Appl. Mech., 38 (1971), 377.
• 6)   D. B.  Bogy: Trans. ASME, J. Appl. Mech., 35 (1968), 460.
• 7)   S.  Ioka,  S.  Kubo,  K.  Ohji and  J.  Kishimoto: Trans. Jpn. Soc. Mech. Eng. A, 60 (1994), No. 572, 141.
• 8)   D. H.  Chen,  K.  Nonomura and  K.  Ushijima: Trans. Jpn. Soc. Mech. Eng. A, 66 (2000), No. 646, 130.
• 9)   T.  Teranishi and  H.  Nishitani: Trans. Jpn. Soc. Mech. Eng. A, 65 (1999), No. 638, 16.