Table of Contents Author Guidelines Submit a Manuscript
Shock and Vibration
Volume 2015 (2015), Article ID 736256, 9 pages
Research Article

Dynamic Analysis of Three-Layer Sandwich Beams with Thick Viscoelastic Damping Core for Finite Element Applications

1Deusto Institute of Technology (DeustoTech), Faculty of Engineering, University of Deusto, Avenida de las Universidades 24, 48007 Bilbao, Spain
2Faculty of Engineering, University of Deusto, Avenida de las Universidades 24, 48007 Bilbao, Spain

Received 12 October 2014; Revised 15 January 2015; Accepted 23 February 2015

Academic Editor: Ahmet S. Yigit

Copyright © 2015 Fernando Cortés and Imanol Sarría. 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.


This paper presents an analysis of the dynamic behaviour of constrained layer damping (CLD) beams with thick viscoelastic layer. A homogenised model for the flexural stiffness is formulated using Reddy-Bickford’s quadratic shear in each layer, and it is compared with Ross-Kerwin-Ungar (RKU) classical model, which considers a uniform shear deformation for the viscoelastic core. In order to analyse the efficiency of both models, a numerical application is accomplished and the provided results are compared with those of a 2D model using finite elements, which considers extensional and shear stress and longitudinal, transverse, and rotational inertias. The intermediate viscoelastic material is characterised by a fractional derivative model, with a frequency dependent complex modulus. Eigenvalues and eigenvectors are obtained from an iterative method avoiding the computational problems derived from the frequency dependence of the stiffness matrices. Also, frequency response functions are calculated. The results show that the new model provides better accuracy than the RKU one as the thickness of the core layer increases. In conclusion, a new model has been developed, being able to reproduce the mechanical behaviour of thick CLD beams, reducing storage needs and computational time compared with a 2D model, and improving the results from the RKU model.

1. Introduction

In the last years many studies have been presented concerning the structural vibration reduction making use of passive damping control techniques by means of surface treatments with viscoelastic materials. A survey of different subjects on viscoelastic treatments can be found in [1]. This kind of vibration control technique is largely used nowadays for several industrial applications, such as aeronautical and automotive components. The free layer damping (FLD) and constrained layer damping (CLD) technologies are two of these viscoelastic surface treatments, consisting of adding a damping viscoelastic layer to the structural system. Specifically, FLD consists of adding that viscoelastic layer on a vibrating metallic base, and the configuration can be analyzed as the flexural behavior of a two-layer beam. In this context, Oberst and Frankenfeld’s model [2] is traditionally used. This model assumes that the layers are thin and the shear effect is negligible. In one of the authors’ works [3] a model was presented in which shear effects are taken into account, improving the results given by Oberst and Frankenfeld’s model for thick layers. On the other hand, CLD technologies are based on a viscoelastic core working under shear deformation. In this configuration, this damping viscoelastic material is in the core of a three-layer sandwich structure. The other two layers are vibrating structural elements, usually made of metallic materials subjected to flexural moment and axial and shear loads. The effect of this shear load is more important as the thickness of the layers increases. Compared with other damping techniques, this viscoelastic treatment presents the disadvantage of adding mass into the original system but offers a better damping-weight ratio than the free layer damping configuration [4, 5].

In Figure 1 the cross section of a three-layer composed beam is represented, the thickness of the base, viscoelastic, and constraining layers being , , and , respectively, and being the width of the beam. Although in CLD beams , in this work three-layer sandwich beams with any value for are considered. The properties for these elastic materials are Young’s moduli and , Poisson’s coefficients and , and densities and , whereas those of the viscoelastic core are , , and , respectively. If the materials are considered isotropic, the shear moduli , , and are related to the extensional ones by , , and for the three layers.

Figure 1: Cross section of a three-layer sandwich beam with a viscoelastic damping layer.

The complex nature of the modulus of viscoelastic materials, represented by the asterisk as , is due to their ability to dissipate mechanical energy. A complex relationship between stress and strain allows representing the hysteretic behaviour of this kind of materials [6]. The dissipative property may be described by the loss factor :where and are the real and imaginary components of the complex modulus , which are known as the storage and loss modulus, respectively, the former being also frequently represented by . For polymers, both storage and loss moduli and consequently loss factor vary with temperature and frequency in terms of three different behaviours: rubbery, vitreous, and transition (see, e.g., [7]). In this sense, Jones [8] reviews experimental data for complex modulus of some typical viscoelastic materials, including elastomers, adhesives, and specially compounded materials. The experimental characterisation of low stiffness damping materials, which are not appropriate to prepare self-supporting test specimens, may be achieved by means of the ASTM E 756-05 standard [9], wherefrom the nomenclature has been taken.

The classical RKU [10] model is one of the first models developed in order to calculate the flexural stiffness for sandwich beams with viscoelastic core and is one of the most used models nowadays, although it may lose accuracy in several cases, for example, when rotational or extensional inertias are not negligible and specifically when a thick viscoelastic core in a CLD beam is studied. This is due to the assumptions which this model takes.(i)The shear deformation in the elastic base and constraining layers is considered to be negligible.(ii)The longitudinal direct stress in the viscoelastic layer is negligible.(iii)The shear stress in the viscoelastic layer is assumed to be uniform.(iv)The plane cross section in each layer remains plane after deformation.There are some other models beyond these restrictions such as [1118] or more recently [19, 20] but they are not so computationally efficient for finite elements applications. In fact, RKU is used in most engineering applications basically due to its very easy computational implementation. Also, RKU strictly applies to simply supported beams, although it is much more generally used, for other boundary conditions, such as the fixed support. We will refer to the results obtained as RKU results, although RKU theory for simply supported beams has been used in this work for a cantilever beam, in the context of a finite element analysis with displacement fields restricted to those of the RKU theory.

In short, following a similar approach as for the study of the FLD configuration for thick beams previously mentioned [3], this work is aimed at developing a new model improving the accuracy of the RKU model for three-layer sandwich beams with a thick viscoelastic core but maintaining its computational benefit. Thus, an equivalent complex flexural stiffness is derived considering quadratic shear stress based on Reddy-Bickford’s theory. In order to prove the improvement achieved by the new model, a numerical application for a cantilever beam is presented comparing the solutions provided by three different finite element models: a 2D model (whose results are considered to be the reference ones) and two 1D models, based on the RKU theory and the new one, respectively. The damping material is characterised by means of a fractional derivative model involving the frequency dependence of the complex modulus, which implies important disadvantages for the dynamic analysis. Thus, the extraction of the eigenvalues and eigenvectors is carried out by a simple and effective iterative algorithm [21]. Finally, the dynamic response of the three models is compared in terms of the frequency response function.

As a result, the new model is able to reproduce the mechanical behaviour of three-layer sandwich beams, reducing storage needs and computational time compared with a 2D model and improving the results provided by the classical RKU model.

2. Homogenised Model for a CLD Beam

Next the theoretical study of a three-layer sandwich beam is presented, in which quadratic shear stress is taken into account. An equivalent flexural stiffness will be deduced for pinned-pinned beams. In order to simplify the notation, it is assumed that the behaviour of all materials is linear and elastic (i.e., all the magnitudes are real) and the complex character of the modulus in the viscoelastic core will be taken into account in the numerical examples of the subsequent sections. This substitution of a complex modulus into an elastic solution is usually known as the “correspondence principle” (see any standard text on viscoelasticity, e.g., [22]).

The equivalent flexural stiffness may be obtained by the addition of the individual contribution of each layer:where , , and are the complex flexural stiffness of the three layers and , , and are the complex cross-sectional second order moments, computed with respect to the neutral axis, given by respectively, in which the complex position of the neutral axis is represented by Following the same methodology as in [3] and decoupling the transverse displacement of any cross section in a term due to the flexural moment and in another term derived by the shearing force (see any book of strength of materials for details, e.g., [23]), the equivalent shear stiffness of the cross section can be found to be decomposed as For the geometry represented in Figure 1, the stiffness of , , and of the individual layers satisfiesrespectively, whererespectively. In these equations, is the flexural stiffness given by (2).

By solving these integrals, it yieldsrespectively, where , , , , , and . For base and constraining layers made of metallic materials, which are much thinner and more rigid than the viscoelastic layer, the terms and may be neglected with respect to the coefficient of the polymeric layer . It can be pointed out that if the three layers were composed of the same material, (14) provides the well-known result for homogeneous rectangular sections:where is the total cross-sectional area.

By considering the shear coefficient , the flexural field equation in free vibration is given by Timoshenko’s formula [24]:where is the mass per unit length. In (18) the rotational inertia has been neglected. When the shear stiffness tends to infinity, (18) degenerates on the well-known Euler-Bernoulli field equation. Following the same methodology as in [3], the equivalent flexural stiffness , considering shear effects and satisfying the homogenised Euler-Bernoulli field equation, can be found to beThe function takes into account the shear effects and is given byIn order to compare RKU and the new model, it can be noted that, for the static values, , both models give the same result, the equivalent . However, as a difference, when frequency tends to infinity , the stiffness of the present model tends to zero whereas that of the RKU one tends to a finite value.

These two models have another common property when shear stiffness tends to infinity, that is, when shear deformations are negligible. In this case, the function tends to zero and the equivalent flexural stiffness tends to the classic .

3. Dynamic Analysis of a Three-Layer Sandwich Beam Using Finite Element Procedures

3.1. Problem Definition

In this section the harmonic analysis of a three-layer sandwich beam in a CLD configuration will be completed using finite element procedure techniques. Three different thicknesses of the viscoelastic core layer will be studied, so as to evaluate the accuracy of the homogenized models compared to a 2D model, whose solution will be considered as exact, considering the nonexistence of experimental results.

The length of the beam is , the width is , the thickness of the base and constraining metallic layers is , and for the viscoelastic layer , and is chosen.

The properties of the materials are taken from the experimental characterisation effectuated by Cortés and Elejabarrieta [25] on AISI T 316 L stainless steel laminated sheet and on Soundown Vibration Damping Tile material [26]. Indeed, Young’s modulus and density of the elastic material are and , respectively, and density of the damping material is . Poisson’s coefficient is chosen for the elastic materials, and is chosen for the viscoelastic material. The experimental data of the storage modulus and loss factor for the viscoelastic core layer were fitted to a four-parameter fractional model [27, 28], given bywhere and represent the relaxed and unrelaxed modulus, respectively; is the relaxation time and is the fractional parameter. The parameter values are summarised in Table 1.

Table 1: Parameters of the fractional derivative model.

The dynamic behaviour of the three-layer sandwich beam in a CLD configuration is studied on the basis of three different finite element models. The first of them is a 2D model discretised in bilinear quadrilateral elements with four nodes under plane-stress assumption (see, e.g., [2931] for details about finite element formulations). All the three layers are modelled with 4 elements along thickness to assure the continuous evolution of the shear stress and with 60 elements along the length (see Figure 2) in order to obtain the first three eigenvalues accurately enough: it has been checked that this number of finite elements is enough to get convergence for any of the results shown in the tables.

Figure 2: Finite element model for the 2D dynamical analysis of a sandwich beam.

The consistent mass matrix and the stiffness matrix obtained using reduced integration with Kosloff and Frazier [32] hourglass control are summarised as follows.

Mass and Stiffness Matrices for the 2D Model Quadrilateral Finite Elements. See Figure 3.

Figure 3: Rectangular Finite Element (Plane-Stress Assumption).

Consistent Mass Matrix . Consider

Stiffness Matrix Obtained by Reduced Integration with Hourglass Control. Consider

Extensional and shear strain and transverse, extensional, and rotational inertias are considered in this 2D finite element model which is then able to reproduce the dynamical behaviour of the three-layer sandwich beam with any thickness of the viscoelastic core.

The two other models are 1D beam models whose th mass and stiffness matrices arerespectively, where is the length of the th finite element. The complex flexural stiffness of (25) is given by (2) for the RKU model and by (19) for the new thick beam model. The discretisation is also made with 60 finite elements along span.

3.2. Extraction of Eigenvalues

The equation from which the complex eigenvalues of the system under study can be obtained iswhere and are the complex eigenvalue and eigenvector of the th mode, respectively, is the mass matrix, and is the complex stiffness matrix, which is dependent on frequency. This is the real part of the square root of the complex eigenvalue :which induces a nonlinearity into the eigenproblem. There are several methods such as those of Lanczos [33] or Arnoldi [34], which use iterative procedures involving important computational time. In order to decrease this computational effort, Cortés and Elejabarrieta [21] developed an iterative procedure that approximates in a simple and accurate way the complex eigenpair. This method begins by considering the static stiffness matrix in (26) yieldingand then the undamped eigensolutions and can be obtained. Then, Nelson’s method [35] is used in order to calculate the eigenvector derivative . From this eigenvector derivative and taking into account the variation of the complex stiffness for the obtained eigenfrequency,and by means of Taylor’s series approach, a complex finite increment of the eigenvector is obtained. The complex eigenvector can be approximated withwith which the complex eigenvalue is estimated according towhere denotes the Hermitian transpose operator, that is, the complex conjugate transposition. Equations (29)–(31) can be iterated making use of the new eigenfrequency given by (27), in order to obtain the desired convergence tolerance. As a main difference with other iterative methods, this one presents the advantage of solving only once the undamped eigenproblem, and the iterations are carried out on the derivatives, reducing computational resources.

If damping in the system is very large, the accuracy of the method can be improved by means of the incremental approach of the method, as seen in [21].

Making use of this new method, with the corresponding incremental approach, the first three modal natural frequencies derived from (27) and loss factor derived fromcan be computed. The corresponding results for the three thicknesses and for the three models under study are shown in Tables 24.

Table 2: Modal properties of the sandwich cantilever beam with  mm.
Table 3: Modal properties of the sandwich cantilever beam with  mm.
Table 4: Modal properties of the sandwich cantilever beam with  mm.

It can be pointed out that the results for the natural frequency are practically the same for the thinnest beam (see Table 2), and the differences between the present model and the RKU one are more important as the thickness of the viscoelastic layer increases, which is an expected behaviour; the reason is that the shear contribution is more important for larger thickness, and the present model considers a more realistic shear stress distribution. Also, the most important differences take place at higher order modes. This is because, at higher frequencies, the shear effects acquire more importance, and, as previously mentioned, the present model takes into account shear effects in a more effective way.

Specifically, for the third mode of the beam with a viscoelastic layer thickness equal to 5 mm (see Table 3), the present model improves the RKU result in a 6.8% (from 10.5% down to 3.7%). This improvement is even more important for the beam with 10 mm of viscoelastic layer (see Table 4), where the difference between the present model and the RKU one with respect to the 2D model goes up to 16.5% (from 21.5% to 5.0%).

As for the results of the modal loss factor , an erratic behaviour in both RKU and thick beam models can be noted. It should be highlighted that this parameter cannot be directly compared, because the damping of the viscoelastic material represented by the loss factor depends on frequency according to (21), and the natural frequencies for the models are not the same. Instead of modal loss factor , the amplitudes of the resonance peaks will be compared in the next section.

3.3. Analysis of the Frequency Response Function

The matrix system for a cantilever beam is given bywhere and are the nodal vectors of the amplitude of the force and the displacement, respectively. The frequency dependence of the complex stiffness matrix involves the fact that the modal superposition cannot be applied, and therefore (33) must be solved for at each desired frequency . Figure 4 represents the frequency response up to 3 kHz of the free-end of the three-layer sandwich beam in all the three models, when a harmonic unitary force is applied on the right-hand side of our system.

Figure 4: Frequency response function for the free-end for the cantilever beam: (a) , (b) , and (c) .

Figure 4(a) indicates that the responses provided by the three models are practically the same, as expected for such a thin beam in which the shear is not so important. On the contrary, Figures 4(b) and 4(c) illustrate a displacement of the curves to the right with respect to that given by the 2D model. This is due to the overestimation of their natural frequencies, as it was mentioned in the previous section. The differences are more patent as the frequency is larger and as the thickness is bigger, mostly for the RKU model.

Regarding the amplitudes of the resonance peaks, Tables 57 show how the model presented in this paper improves the results of the RKU model, the improvement being better as the thickness of the viscoelastic layer increases. For instance, for the first mode of the thickest beam, the difference between the errors in RKU and present models goes down from 7.3% to 3.4%, which is a significant difference taking into account the logarithmic scale.

Table 5: Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  mm.
Table 6: Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  mm.
Table 7: Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  mm.

Since a quadratic distribution of the shear stress is representative of the state of stress in beams, the differences between the new model and the 2D model are mainly due to the inertias.

4. Conclusions

In this paper a formulation for three-layer sandwich beams has been presented, which takes into account the deflection due to shearing forces. For this, quadratic distribution of shear stress was considered along thickness. In order to test the performance of the presented model, a dynamical analysis has been carried out on a CLD beam, in which the viscoelastic material has been modelled by means of a fractional derivative model. Consequently, the natural frequencies and resonance peak amplitudes provided by the new model and by the classical RKU model have been compared with those given by a 2D finite element model. As a conclusion, both models present similar results for a thin core layer, but the results are significantly improved as the thickness of the core layer increases. This formulation is very simple to be implemented in finite element analysis, presenting the same complexity level as the RKU one. This represents an important advantage with respect to other formulations.

In short, a new model for finite element calculations in three-layer sandwich beams with viscoelastic core layer has been developed, improving the results of the RKU one, because the former considers shearing force deflection, without enlarging computational efforts.

Conflict of Interests

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


  1. F. Cortés, M. Martinez, and M. J. Elejabarrieta, Viscoelastic Surface Treatments for Passive Control of Structural Vibration, Nova Publishers, New York, NY, USA, 2012.
  2. H. Oberst and K. Frankenfeld, “Über die dämpfung der biegeschwingungen düner bleche durch fest haftende beläge,” Acustica, vol. 2, pp. 181–194, 1952. View at Google Scholar
  3. F. Cortés and M. J. Elejabarrieta, “Structural vibration of flexural beams with thick unconstrained layer damping,” International Journal of Solids and Structures, vol. 45, no. 22-23, pp. 5805–5813, 2008. View at Publisher · View at Google Scholar · View at Scopus
  4. A. D. Nashif, D. I. G. Jones, and J. P. Henderson, Vibration Damping, John Wiley & Sons, New York, NY, USA, 1995.
  5. C. T. Sun and Y. P. Lu, Vibration Damping of Structural Elements, Prentice Hall, Upper Saddle River, NJ, USA, 1995.
  6. N. O. Myklestad, “The concept of complex damping,” Journal of Applied Mechanics, vol. 19, pp. 284–286, 1952. View at Google Scholar
  7. I. M. Ward and D. W. Hadley, An Introduction to the Mechanical Properties of Polymers, John Wiley & Sons, Chichester, UK, 1993.
  8. D. I. G. Jones, Handbook of Viscoelastic Vibration Damping, John Wiley & Sons, Chichester, UK, 2001.
  9. American Society for Testing and Material, “Standard test method for measuring vibration-damping properties of materials,” ASTM E 756-05, 2005. View at Google Scholar
  10. D. Ross, E. M. Kerwin, and E. E. Ungar, “Damping of plate flexural vibration by means of viscoelastic laminae,” in Structural Damping, Section II, ASME, New York, NY, USA, 1959. View at Google Scholar
  11. R. A. DiTaranto, “Theory of vibratory bending for elastic and viscoelastic layered finite-length beams,” Journal of Applied Mechanics, vol. 32, no. 4, p. 881, 1965. View at Publisher · View at Google Scholar
  12. D. J. Mead and S. Markus, “The forced vibration of a three-layer, damped sandwich beam with arbitrary boundary conditions,” Journal of Sound and Vibration, vol. 10, no. 2, pp. 163–175, 1969. View at Publisher · View at Google Scholar · View at Scopus
  13. M. J. Yan and E. H. Dowell, “Governing equations for vibrating constrained–layer damping of sandwich beams and plates,” Journal of Applied Mechanics, vol. 39, no. 4, pp. 1041–1046, 1972. View at Publisher · View at Google Scholar · View at Scopus
  14. Y. V. K. Sadasiva Rao and B. C. Nakra, “Vibrations of unsymmetrical sandwich beams and plates with viscoelastic cores,” Journal of Sound and Vibration, vol. 34, no. 3, pp. 309–326, 1974. View at Publisher · View at Google Scholar · View at Scopus
  15. M. Macé, “Damping of beam vibrations by means of a thin constrained viscoelastic layer: evaluation of a new theory,” Journal of Sound and Vibration, vol. 172, no. 5, pp. 577–591, 1994. View at Publisher · View at Google Scholar · View at Scopus
  16. Z. Liu, S. A. Trogdon, and J. Yong, “Modeling and analysis of a laminated beam,” Mathematical and Computer Modelling, vol. 30, no. 1-2, pp. 149–167, 1999. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. J. M. Bai and C. T. Sun, “Effect of viscoelastic adhesive layers on structural damping of sandwich beams,” Mechanics of Structures and Machines, vol. 23, no. 1, pp. 1–16, 1995. View at Publisher · View at Google Scholar · View at Scopus
  18. Y. Frostig and M. Baruch, “Free vibrations of sandwich beams with a transversely flexible core: a high order approach,” Journal of Sound and Vibration, vol. 176, no. 2, pp. 195–208, 1994. View at Publisher · View at Google Scholar · View at Scopus
  19. A. C. Galucio, J.-F. Deüi, and R. Ohayon, “Finite element formulation of viscoelastic sandwich beams using fractional derivative operators,” Computational Mechanics, vol. 33, no. 4, pp. 282–291, 2004. View at Publisher · View at Google Scholar · View at Scopus
  20. Y. Frostig and O. T. Thomsen, “High-order free vibration of sandwich panels with a flexible core,” International Journal of Solids and Structures, vol. 41, no. 5-6, pp. 1697–1724, 2004. View at Publisher · View at Google Scholar · View at Scopus
  21. F. Cortés and M. J. Elejabarrieta, “An approximate numerical method for the complex eigenproblem in systems characterised by a structural damping matrix,” Journal of Sound and Vibration, vol. 296, no. 1-2, pp. 166–182, 2006. View at Publisher · View at Google Scholar · View at Scopus
  22. R. M. Christensen, Theory of Viscoelasticity, Academic Press, London, UK, 1982.
  23. T. H. G. Megson, Structural and Stress Analysis, Butterworth-Heinemann, Oxford, UK, 2005.
  24. S. P. Timoshenko, “LXVI. On the correction for shear of the differential equation for transverse vibrations of prismatic bars,” Philosophical Magazine Series 6, vol. 41, no. 245, pp. 744–746, 1921. View at Publisher · View at Google Scholar
  25. F. Cortés and M. J. Elejabarrieta, “Viscoelastic materials characterisation using the seismic response,” Materials & Design, vol. 28, no. 7, pp. 2054–2062, 2007. View at Publisher · View at Google Scholar · View at Scopus
  26. Soundown Corporation, November 2012,
  27. R. L. Bagley and P. J. Torvik, “On the fractional calculus model of viscoelastic behaviour,” Journal of Rheology, vol. 30, no. 1, pp. 133–155, 1986. View at Publisher · View at Google Scholar · View at Scopus
  28. T. Pritz, “Analysis of four-parameter fractional derivative model of real solid materials,” Journal of Sound and Vibration, vol. 195, no. 1, pp. 103–115, 1996. View at Publisher · View at Google Scholar · View at Scopus
  29. O. C. Zienkiewicz and R. L. Taylor, Finite Element Method, Butterworth-Heinemann, Oxford, UK, 2000.
  30. T. J. R. Hughes, The Finite Element Method, Dover Publications, New York, NY, USA, 2000.
  31. K. J. Bathe, Finite Element Procedures, Prentice Hall, Upper Saddle River, NJ, USA, 1996.
  32. D. Kosloff and G. A. Frazier, “Treatment of hourglass patterns in low order finite element codes,” Numerical and Analytical Methods in Geomechanics, vol. 2, no. 1, pp. 57–72, 1978. View at Publisher · View at Google Scholar
  33. C. Lanczos, “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators,” Journal of Research of the National Bureau of Standards, vol. 45, pp. 255–282, 1950. View at Publisher · View at Google Scholar · View at MathSciNet
  34. W. E. Arnoldi, “The principle of minimized iteration in the solution of the matrix eigenvalue problem,” Quarterly of Applied Mathematics, vol. 9, pp. 17–29, 1951. View at Google Scholar · View at MathSciNet
  35. R. B. Nelson, “Simplified calculation of eigenvector derivatives,” AIAA Journal, vol. 14, no. 9, pp. 1201–1205, 1976. View at Publisher · View at Google Scholar · View at MathSciNet