International Journal of Aerospace Engineering
Volume 2021 (2021), Article ID 6653824, 11 pages
https://doi.org/10.1155/2021/6653824

Research Article

Thermal Structure Strength Analysis of Nozzle of Solid Rocket Motor with the Coupled Algorithm

Chunguang Wang,¹ Weiping Tian,² and Kaining Zhang¹

¹State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, China
²The Fourth Academy of China Aerospace Science and Technology Corporation, Xi’an 710025, China

Correspondence should be addressed to Weiping Tian; tianweiping1964@163.com

Received 22 December 2020; Revised 9 April 2021; Accepted 1 May 2021; Published 25 May 2021

Academic Editor: Jiaqiang E

Copyright © 2021 Chunguang Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In order to achieve the thermal structural integrity analysis of the solid rocket motor nozzle accurately and efficiently, the multifield (flow-thermal-mechanical) coupled numerical investigation was carried out based on the mesh-based parallel code coupled interface. The numerical simulation process and finite element model of the coupled algorithm and engineering algorithm were obtained, while the physical model was simplified appropriately. The coupled interface parameters, internal flow field, temperature field, and stress field of the coupled algorithm were compared with the engineering algorithm results, and the effectiveness and accuracy of the numerical simulation were validated. The numerical investigations shown that both the temperature field and stress field obtained by the coupled algorithm were slightly lower than which obtained by the engineering algorithm. These were considered to be impacted by the Bartz empirical formula and the one-dimensional isentropic flow assumption. Further experimental investigations shown that the exterior surface temperature and strain of the nozzle throat obtained by the coupled algorithm were much closer to the experimental results, which further verified the accuracy of the coupled algorithm.

1. Introduction

After the ignition of the solid rocket motor, the high-temperature and high-pressure gas generated from the propellant burning surface flow out through the nozzle. The severe thermal environment is the main factor causing the structural failure, while the thermal structure analysis is the key to evaluate the reliability of the nozzle. The thermal protection effect and thermal stress value of the nozzle components can be obtained by the numerical simulation. Furthermore, the accurately prediction of the nozzle thermal structure integrity can be obtained by the combination of the numerical and experimental results, which could be used in the structure optimization, as well as the selection of material and performance indicators [1].

At present, the most common calculation methods for the initial boundary condition of the nozzle temperature field and stress field include

(1)Calculated by the analysis software(2)Calculated by the engineering Bartz empirical formula

In reference [2], the finite element method, in the form of the commercial finite element code ADINA, is used to investigate the dynamic thermostructural response of a composite rocket nozzle throat. In reference [3], the coupling simulation of heat transfer and transient temperature of the rocket nozzle wall is carried out. Equations of radiative heat flux on surfaces in an enclosure with inhomogeneous, participating media is developed to compute the radiative heat flux. In reference [4], the finite element model of the complex nozzle was established. The CFD and the thermal and structural analysis part of ANSYS were used to simulate the temperature field and stress field, and the coupled calculation method of the thermal structure was achieved. In order to solve the flow field and obtain the wall temperature distributions, numerical models have been developed incorporating both solid and fluid regions in reference [5]. In reference [6], the variation law of temperature and stress in divergent section at motor working was described by calculating the temperature field and stress field. In order to determine temperature and pressure in the nozzle, the gas was regarded as one-dimension isentropic flow. Based on the axial symmetric finite element model, the transient temperature field in the nozzle was calculated.

In reference [7], numerical simulation was made on the flow field inside the nozzle with the Euler-Lagrange two-phase flow model. Different overloads and particle’s diameter were considered. On the basis of calculation, the temperature field of the 3D thermal structure of the whole nozzle was analyzed with consideration of three heat transfer styles and heat transfer characteristics of each region of the nozzle. To simulate the fluid field and wall temperature distribution, numerical models have been developed considering both solid and fluid regions in reference [8]. The two parts work together. Each one provides a boundary condition for the other, and the solution of the coupled problem has been resolved. In order to estimate the heat protection performance and obtain the temperature load for erosion calculation and structural stress analysis, temperature numerical simulation of the composite nozzle thermostructure was implemented by the Fluent CFD software using the fluid-structure coupled method in reference [9]. In the computation, the RNG k-ω turbulence model and the enhanced wall treatment method were selected. In reference [10], the Reynolds stress viscous model and enhanced wall function are used to solve Navier-Stokes equations. By setting the inner wall temperature distribution as boundary condition, the two-dimensional axisymmetric transient heat transfer equation is solved. By loading the temperature field by step time, a transient static structure analysis is proposed. The flow field, temperature field, and stress field are obtained. The nonlinear transient thermal structure coupling analysis of the nozzle is carried out by using MSC.MARC software in reference [11].

The calculation methods in these references have basically solved the problems in the project, but there is a big gap between the calculation results and the experimental results.

The feedback effect of the solid structure on the temperature field is not taken into account in method (1), which resulting in the decoupled improper [12, 13]. And the calculation of convective heat transfer coefficient by the Bartz empirical formula is inaccurate in method (2), which may result in large calculation errors. Therefore, in order to achieve the thermal structural integrity analysis accurately and efficiently, plenty of researches have been put on the multifield coupled analysis of the solid rocket motor nozzle. In reference [14], a numerical model of fluid-solid integration was established. The three-phase flow-thermal-mechanical coupled analysis of the solid rocket motor nozzle was realized by using the fluid-solid/thermal coupled analysis method of strong coupled of gas flow/solid thermal conductivity and weak coupled of structural stress. The specific method is to embed the thermal structure calculation program of the solid rocket motor nozzle into FLUFNT software based on UDF technology. The coupled heat transfer results were extracted by UDF macrocommand. The automatic identification program of contact boundary was compiled to simulate the nonlinearity contact between nozzle interfaces by the direct constraint method and to realize the integrated calculation of the solid rocket motor nozzle thermal structure. Although the method achieved the three-phase coupled analysis of flow-thermal-mechanical, its inconvenience and low adaptability result in the limited application.

In this paper, the coupled algorithm applied in the investigation of the solid rocket motor nozzle flow-thermal-mechanical field are proposed and validated based on the analysis software of FLUENT, ABAQUS, and MpCCI. Furthermore, the coupled algorithm is compared with the engineering algorithm and the experimental date. The results show that the temperature field and stress field obtained by the coupled algorithm are slightly lower than the engineering algorithm, which is considered to be impacted by the Bartz empirical formula and the one-dimensional isentropic flow assumption. Further investigations show that the nozzle temperature and strain obtained by the coupled algorithm are closer to the experimental data, which further verified the accuracy of the coupled algorithm.

2. Research Model

2.1. Physical Model

The nozzle is composed of metal case, thermal insulation in convergent section, throat insert, thermal insulation in back wall, thermal insulation in divergent section, and so on. The structure of the nozzle in this paper is shown in Figure 1. The diameter of the throat is 210 mm, the expansion ratio of the nozzle is 12, and the total length of the nozzle is 850 mm. The material parameters of each part of nozzle are shown in Table 1.

Figure 1: Two-dimensional axisymmetric model of the nozzle.

Table 1: Performance parameters of nozzle materials.

Component
Throat insert
Thermal insulation in divergent section
Thermal insulation in the back wall
Thermal insulation in convergent section
Case (shell)

Material
Puncture C/C
Carbon cloth winding
High silicon oxygen cloth winding
Carbon Fiber molding
30CrMnSiA

Modulus GPa
9
12.2
17
10.8
210

Poisson ratio
0.1
0.1
0.12
0.1
0.3

1.95
1.45
1.65
1.4
7.8

Specific heat J/(kg·K)
2100
920
850
1000
473

Thermal conductivity W/(m·K)
85
1
0.6
1.5
50

⁶
1
4.2
12
10
13

2.2. Details of the Model

The research model is defined as follows:

(a)Propellants: three-component HTPB propellants(b)Time: 42 s(c)Average pressure: 8.6 MPa (standard state)(d)Combustion temperature: 3600 K(e)Gas specific heat ratio: 1.16(f)Average molecular weight of gas: 29.06(g)Gas specific heat at constant pressure: 3900 J/kg·K

2.3. Basic Assumptions

For the multifield analysis of flow-thermal-mechanical, some minor factors can be neglected. The necessary assumptions for the nozzle are as follows:

(1)The gas is simplified as ideal gas, and its parameters (such as temperature and pressure) will not change with time, and the radiation heat transfer is not considered(2)The assembly clearance between the nozzle components is ignored(3)The contact thermal resistance between the components is ignored(4)The effects of carbonization, ablation, and radiation are ignored. Only the convective heat transfer between the gas and nozzle is taken into consideration(5)The heat transfer between the nozzle exterior surface and the environment is natural convection

3. Multifield Coupled Model Based on the Coupled Algorithm

3.1. Coupled Numerical Simulation Process

In the fluid domain, temperature variations are related to density, velocity, and temperature gradient. In the solid domain, temperature variations are only related to temperature gradient while the velocity is zero and density is constant. In the interface between the fluid domain and solid domain, the temperature should be solved to satisfy both the solution in the fluid domain and the solution in the solid domain. Therefore, in order to acquire the accurate temperature field, the analysis of the flow field in the fluid domain should be coupled with the analysis of the heat transfer in the solid domain. The interface between the fluid domain and the solid domain is set as the coupled boundary. The specific coupled numerical simulation process is shown in Figure 2.

Figure 2: Coupled calculation process.

3.2. Fluid Domain Calculation Model

Since the heat transfer in the solid domain depends on the flow inside the nozzle, accurate simulation of the fluid domain is significant important. The simulations aimed to solve the two-dimensional, coupled, implicit Reynolds Averaged Navier-Stokes (RANS) equations. A two-equation RNG -ε turbulence model was used to approximately evaluate Reynold’s stress and obtain a reasonable estimation of the flow features.

The RNG -ε model is derived using a statistical method called renormalization group theory. It is similar in form to the standard -ε model, but includes the following improvements:

(1)RNG model adds a term to its ε equation, which improves the accuracy of high-speed flow(2)The RNG model takes into account the effect of the eddy current on turbulence and improves the accuracy of eddy flow(3)RNG theory provides an analytical formula for turbulent Prandtl numbers, while the standard -ε uses a user-specified constant value

Although the standard -ε model is a high Reynolds number model, the RNG differential formula theory provides a way to obtain an effective viscosity from the analysis, taking into account the effects of low Reynolds numbers; however, the effective use of this property depends on proper treatment of the near-wall region. These characteristics make the RNG -ε model more accurate and reliable over a wider range of flows than the standard -ε model.

The governing equations were discretized by the second order spatially accurate upwind scheme (SOU). The Courant–Friedrichs–Lewy (CFL) number was maintained at 1.0 with using proper underrelaxation factors to ensure stability. Localized grid refinement was adopted in the near wall domains in order to capture as much flow details as possible and enhance the accuracy of the numerical simulations. The flow-thermal-mechanical couple boundary condition wall condition was applied at the nozzle wall. The entire computational domain was consisted of hexahedral elements.

The mesh height of the first layer is 0.1 mm, which grows to the mainstream area in the proportion of 1.15, and the grid height of the mainstream area is 2 mm.

Figure 3 depicts computational domains and boundary conditions. The time step applied in this paper is set as 0.01 s, and a total of 4200 steps are iterated to complete the calculation of a 42 s working process.

Figure 3: Fluid domain calculation model.

The -ε model is often used with the wall function, and so the requirement for the mesh height in the boundary layer is not very high, and the calculation accuracy can be guaranteed when the Yplus of the wall is around 30. In the calculation results in this paper, the variation trend of Yplus along the nozzle wall is shown in Figure 4, which verifies the rationality of the mesh in this paper and the correctness of the turbulence model selected in this paper.

Figure 4: The trend of the wall Yplus along the axial distance.

3.3. Solid Domain Calculation Model

The simplified axisymmetric model was also applied in the solid domain calculation. Localized grid refinement was adopted in the stress concentration area. The selected calculation element CAX4RT which taken the calculation of temperature and stress into account could obtain more accurate results. The calculation model is shown in Figure 5.The motor and nozzle are treated as binding constraints. The throat insert and insulation are massaged and rubbed while the friction coefficient is 0.25. Besides, the other interfaces are set as binding. The material parameters of each component were shown in Table 1. The coupled tempdisplacement (transient) analysis step was applied, and the calculation time was 42 seconds.

Figure 5: Solid domain computing model.

4. Calculation Model Based on the Engineering Algorithm

4.1. Engineering Algorithm for Wall Boundary Conditions

The convective heat transfer coefficient, wall recovery temperature, and wall pressure distribution are necessary in calculating the nozzle temperature distribution and thermal stress field by the engineering algorithm. The convective heat transfer coefficient is determined by Bartz formula (1):

where is the nozzle throat insert diameter, is the gas dynamic viscous coefficient, is the gas constant pressure specific heat capacity, is the gas Prandtl number, is the gas total pressure at nozzle inlet, is the gas characteristic velocity, is the radius of nozzle throat insert curvature, is the nozzle throat insert area, is the area of the computed section, and is the correction factor of the convective heat transfer coefficient.

The quasione-dimensional isentropic flow in the nozzle can be solved by the following equation.

is the specific heat ratio of gas, is the gas total temperature at nozzle inlet, is the gas total pressure at nozzle inlet, is wall temperature at the computed section, is the temperature at the computed section, is the March number at the computed section, and the is area of the computed section.

4.2. Computational Model

The computational model, material properties, boundary conditions, and analysis steps adopted in this chapter keep the same with the method in chapter 3.3. The heat flux and pressure are loaded in the nozzle interior surface. The subroutines of ABAQUS, namely, DFLUX and DLOAD, are applied to complete the calculation and iteration of (1)~(5) and (6)~(8).

is the heat flux of the section, is the gas recovery temperature, and is the coefficient of recovery.

5. Results and Discussion

5.1. Fluid Domain

The pressure contour, temperature contour, and Mach number contour of the flow field are shown in Figure 6. The calculation results are basically consistent with those in references [15, 16], which verified that the method of the calculating nozzle flow field by fluent in this paper was feasible.

(a)

(b)

(c)

Figure 6: Numerical results of the flow field: (a) pressure, (b) temperature, and (c) Mach number.

5.2. Solid Domain

The nozzle temperature field at 42 s is obtained by numerical simulation and compared with the engineering algorithm, as shown in Figure 7.The temperature distribution of this two methods are basically the same, and no temperature rise on the nozzle exterior surface during the its working time. Besides, the temperature of nozzle interior surface obtained by the coupled algorithm is significantly lower than that obtained by the engineering algorithm. The potential reason is that there are some errors in the calculation of convective heat transfer coefficient by Bartz formula in the engineering algorithm. In addition, some errors are introduced in the hypothesis of the one-dimensional isentropic flow assumption inside the nozzle as well, which resulting in bigger calculation of heat flux density and leading to excessive temperature rise in the nozzle wall.

(a)

(b)

Figure 7: Comparison of nozzle temperature distribution using two algorithms: (a) coupled algorithm and (b) engineering algorithm.

The stress fields of nozzle throat insert at 42 s are showed in Figure 8, which concluding the Mises stress, axial stress, radial stress, and circumferential stress, respectively. It can be easily found that the stress distributions of these two algorithms are almost the same except the circumferential stress. Besides, the stress values obtained by the engineering algorithm are slightly larger than those obtained by the coupled algorithm. The potential reason is that the engineering algorithm results in the higher temperature rising of nozzle throat insert, which finally leading to the higher thermal stress influenced by the circumferential restraint of back wall and metal case.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8: 42 s stress of C/C throat insert in motor operation: (a) coupled algorithm for Mises stress, (b) engineering algorithm for Mises stress, (c) coupled algorithm for radial stress, (d) engineering algorithm for radial stress, (e) coupled algorithm for axial stress, (f) engineering algorithm for axial stress, (g) coupled algorithm for circumferential stress, and (h) engineering algorithm for circumferential stress.

In addition, in these two calculation results, the throat insert stress range is (-86 MPa, 12 MPa), and the maximum stress is the cusp stress concentration. However, the actual stress range should be (-51 MPa, 10 MPa). Both calculation results are within the range of permissible C/C throat insert stress (-90 MPa, 70 MPa). The throat insert structure has high reliability.

6. Comparisons and Analysis of Numerical and Experimental Results

After the SRM thermal test, the C/C throat insert was dissected and analyzed, as shown in Figure 9. As can be seen from the section, the C/C throat insert structure is complete without defects and cracks, which proves that the above judgment results are correct.

(a)

(b)

Figure 9: Solid rocket motor ground test: (a) the nozzle flow field and (b) C/C throat insert cross section after the SRM test.

The nozzle structure in this paper has undergone the examination of the ground test, and the temperature and strain changes of the nozzle exterior surface have been monitored during the experiment. In this paper, the temperature and stress of throat insert are concerned, and the monitoring position is selected as shown in Figure 10. The temperature change at point A is monitored for 120 seconds, and the strain change at that point during the working process of the motor is monitored as well.

Figure 10: Temperature and strain measurement points of the motor ground test.

The measured temperature curve after motor ignition signal to 120 s is extracted and compared with the results calculated by the numerical method as shown in Figure 11.

Figure 11: Temperature variations in exterior case of nozzle throat insert.

It can be found that there is almost no temperature rise in the motor working process (42 s). The results of these three methods are basically the same. The temperature of the nozzle throat case begins to rise at 50 s. The exterior case temperature rise of the nozzle throat obtained by the engineering algorithm is obviously higher than the measured value by the ground test. In contrast, the temperature curve obtained by the coupled algorithm is more consistent with the actual situation. The reason is that the engineering algorithm cannot accurately give the convective heat transfer coefficient after the motor finished working. The Bartz formula and the one-dimensional isentropic flow assumption are difficult to apply to the current nozzle flow state. The coupled algorithm simulates the exhaust process of the nozzle and assigns the boundary of the flow field to the solid boundary through the coupled boundary. This method can obtain more accurate calculation results.

The test results of circumferential and axial strain at point A are extracted and compared with those obtained by the simulation analysis method as shown in Table 2.

Table 2: Comparisons of strain measurements with numerical results.

Method
Circumferential strain
Axial strain

Coupled algorithm
1200
264

Engineering algorithm
822
171

Experimental result
1254
277

It can be found that the strain data obtained by the coupled algorithm are closer to the measured values by the ground test, while the engineering algorithm error is much larger. The reasons can be as follows: there are two main factors causing the nozzle deformation: one is the expansion of the components due to the increase of temperature; the other is that under the internal pressure of the motor, the nozzle wall acts on the pressure varying along the axis, resulting in the overall expansion of the nozzle. The analysis shows that the second is the main factor in overall expansion of the nozzle. In the engineering algorithm, the pressure distribution on the nozzle wall is calculated based on one-dimensional isentropic flow. In the coupled algorithm, the wall pressure is calculated directly by the coupled Algorithm, and the result is more accurate.

7. Conclusions

The main conclusions of this paper are as follows:

(1)Based on the coupled algorithm, the flow-thermal-mechanical multifield coupled analysis of the solid rocket motor nozzle is carried out. Compared with the traditional engineering algorithm, the coupled calculation process is much closer to the actual working conditions. The reason is that the flow field parameters of the coupling algorithm are more accurate than those of the engineering algorithm. This method not only considers the heat transfer effect of the flow field on the solid wall but also considers the effect of the solid wall on the flow field(2)The temperature and stress obtained by the engineering algorithm are slightly higher than those obtained by the coupled algorithm. The analysis shows that the error mainly comes from the Bartz formula and the one-dimensional isentropic flow assumption. The temperature and pressure calculated by the Bartz formula and one-dimensional isentropic flow formula are higher than those calculated by the actual flow field, and so the thermal stress is higher than that obtained by the coupling algorithm(3)Two numerical methods show that the stress of the throat insert is within the required stress range of the material, while the reliable thermal structure of the throat insert is validated by the ground test. According to the analysis of numerical results, the thermal stress is within the allowable stress range of the material (-90 MPa, 70 MPa), and the integrity of the structure can be guaranteed. It means that both methods can be used to analyze the thermal structural integrity of nozzles(4)The temperature and strain of the nozzle throat external case are extracted and compared with the ground test results. The comparison shows that the engineering algorithm has a large error, and the coupled algorithm is closer to the ground test results, which further verifies the accuracy of the coupled algorithm. With the increase of time, the error of the engineering algorithm will be gradually magnified. The engineering algorithm is only suitable for the initial feasibility evaluation of the design scheme, and the fluid-structure coupling algorithm proposed in this paper is more suitable for accurate design(5)As the Bartz equation does not take account the influence of combustion zone distribution in the thrust chamber, thickness variation of the boundary layer, and other actual situations on fuel gas heat flux, it cannot well present the gas heat flux density distribution inside combustion chamber, especially in the combustion zone near the nozzle

Data Availability

As most of the data in this manuscript were related to trade secrets, I cannot provide them completely. In the future, if necessary, I can share some data with reviewers or readers.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This article was funded by the Xi’an Jiaotong University, China.

References

R.-x. Chen, M.-c. Liu, and Z.-m. Li, Design and Research of Solid Rocket Motor, Aerospace Publishing House, Beijing, 2007.

E. V. Morozov and J. F. P. P. de la Beaujardiere, “Numerical simulation of the dynamic thermostructural response of a composite rocket nozzle throat,” Composite Structures, vol. 91, no. 4, pp. 412–420, 2009.

X.-y. Zhang, “Coupled simulation of heat transfer and temperature of the composite rocket nozzle wall,” Aerospace Science and Technology, vol. 15, no. 5, pp. 402–408, 2011.

T. B. Wang, T. S. Xue, C. S. Zhou, X. Chen, and J. Zheng, “Numerical calculation of temperature and stress field of complex nozzle,” Journal of Ballistics, vol. 24, no. 2, pp. 88–105, 2012.

H. Hai-feng, F.-T. Bao, Y.-J. Wang, Q. Cai, and Y. Liu, “Nozzle flow separation fluid-thermal-structure load transfer coupled analysis and test research,” Journal of Astronautics, vol. 32, no. 7, pp. 1534–1542, 2011.

J.-h. Hu, S.-H. Meng, and X.-L. Chang, “Temperature and stress field analysis of conical trapped joint C/C nozzle divergent section,” Journal of Solid Rocket Technology, vol. 35, no. 1, pp. 64–68, 2012.

F. Gao, J.-h. Wang, and H.-y. Xiao, “Numerical simulation for thermal structure of tail-pipe nozzle under high overload,” Journal of Projectiles, Rockets, Missiles and Guidance, vol. 34, no. 4, pp. 111–116, 2014.

H.-y. Huang and Q. Wang, “Coupled numerical study of flow field and heated structure for axisymmetric nozzles,” Journal of Propulsion Technology, vol. 29, no. 2, pp. 223–229, 2008.

F. Xue-jin, H. Deng, and J.-Q. Kou, “Coupled numerical simulation on flow field and thermo-structure of composite nozzle,” Journal of Solid Rocket Technology, vol. 37, no. 6, pp. 781–786, 2014.

S. Lin, F.-t. Bao, and W. Shi, “Thermo-structural analysis of spinning SRM nozzle,” Aeronautical Computing Technique, vol. 44, no. 3, pp. 62–65, 2014.

S.-S. Pei, Y.-c. Hu, and S.-l. Gao, “Nonlinear transient thermal-structure coupling analysis of solid rocket motor nozzle based on MSC.Marc,” Machine Design & Research, vol. 45, no. 4, pp. 6–10, 2018.

X.-y. Z. F.-m. Chen and F.-c. Cai, “Strategy of thermo-structure coupled computation for composite nozzle,” Journal of Aerospace Power, vol. 26, no. 1, pp. 223–227, 2011.

X.-y. Zheng, W.-x. Wang, and F.-c. Cai, “Thermo-structure integrative computation for nozzle in SRM,” Journal of Aerospace Power, vol. 27, no. 9, pp. 2122–2127, 2012.

X.-y. Zheng, W.-x. Wang, and M.-y. Shao, “Study on load transmission line and thermo-structure for composite nozzle,” Chinese Journal of Applied Mechanics, vol. 32, no. 4, pp. 618–624, 2015.

M. Yumuşak, “Analysis and design optimization of solid rocket motors in viscous flows,” Computers & Fluids, vol. 75, no. 75, pp. 22–34, 2013.

R. Dunlap, P. G. Willoughby, and R. W. Hermsen, “Flowfield in the combustion chamber of a solid propellant rocket motor,” AIAA Journal, vol. 12, no. 10, pp. 1440–1442, 1974.