Abstract
This paper alternatively derives the exact element stiffness equation for a beam on Kerr-type foundation. The shear coupling between the individual Winkler-spring components and the peripheral discontinuity at the boundaries between the loaded and the unloaded soil surfaces are taken into account in this proposed model. The element flexibility matrix is derived based on the virtual force principle and forms the core of the exact element stiffness matrix. The sixth-order governing differential compatibility of the problem is revealed using the virtual force principle and solved analytically to obtain the exact force interpolation functions. The matrix virtual force equation is employed to obtain the exact element flexibility matrix based on the exact force interpolation functions. The so-called “natural” element stiffness matrix is obtained by inverting the exact element flexibility matrix. One numerical example is utilized to confirm the accuracy and the efficiency of the proposed beam element on Kerr-type foundation and to show a more realistic distribution of interactive foundation force.
1. Introduction
As a numerical counterpart of the continuous medium model, the continuum finite element model has been widely used by geotechnical researchers in studying several complex soil-structure interaction (SSI) problems due to drastic advances in computer technology. The problem of beams on deformable foundation is the most commonly encountered SSI problem and has many applications in engineering and science [1–3]. Even though the continuum finite element model yields the most comprehensive data on the stress and deformation variations within the beam-foundation system, there is still a substantial need in routine engineering practice to use the mechanical subgrade model to analyze and design the beam-foundation system. This lies in the fact that considerable experience and skill of practicing geotechnical engineers are required in constructing a suitable continuum element mesh, interpreting the analysis results, and implementing the numerical model. These could limit the model access by practicing geotechnical engineers. Furthermore, only small beam-foundation systems can be realistically investigated using the continuum finite element model due to high computational costs, and the beam response along the beam-foundation interface, not the stresses or strains inside the foundation medium, is of high interest by the designing engineers. Therefore, most structural-analysis platforms available in the industry still employ the mechanical subgrade model to represent the supporting foundation with a reasonable degree of accuracy.
The Winkler foundation model [4] is the most rudimentary mechanical subgrade model and has been widely adopted in studying the problem of beams on elastic foundation. In the Winkler foundation model, a set of 1D independent springs is attached along the beam to form the beam-foundation system. This type of foundation model is often referred to as the “one-parameter” foundation model since it is characterized only by the vertical stiffness of the Winkler-foundation springs. Though simple, the Winkler foundation model can lead to some peculiar responses of many practical beam-foundation systems due to omission of the shear stress inside the foundation medium [5, 6]. This omission results in the uncoupling of the individual Winkler foundation springs and the neglect of the existence of the foundation medium beyond either end of the loaded beam and leads to an unrealistic abrupt change in the foundation surface displacement between the loaded and the unloaded regions. To bridge the gap between the continuum finite element model and the rather crude Winkler foundation model, several researchers [7–9] had improved the Winkler foundation model by introducing the second foundation parameter to account for the existence of shear stress inside the foundation medium, resulting in the so-called “two-parameter” foundation model. Even though each researcher group has its own particular way to visualize the second foundation parameter, its proposed expressions for the interactive foundation forces can simply be written in the same mathematical form. For example, Filonenko-Borodich [7] regarded the second foundation parameter as the magnitude of the pretensioned force in the elastic membrane inserted between the beam and Winkler-foundation springs. A detailed discussion of several two-parameter foundation models can be found in Kerr [6].
To further improve the two-parameter foundation model, Hetényi [10] and Kerr [11] had added the third foundation parameter, leading to the so-called “three-parameter” foundation model. The major role of the third foundation parameter is to provide more flexibility in controlling the degree of foundation-surface continuity between the loaded and the unloaded regions of the beam-foundation system. This is in compliance with the observation made by Foppl [12] that the foundation-surface displacement outside the loaded region predicted by the continuous medium model decreased too gradually as opposed to what happened in reality, and hence a certain degree of discontinuity at the loaded-unloaded boundary existed. Furthermore, Kerr [6] concluded that for several types of foundation materials (e.g., soil, soft filament, foam, etc.), neither the Winkler-foundation model nor the continuous medium model can realistically represent the interaction mechanisms between the beams and the contacting media. Among several three-parameter foundation models, the Kerr-type foundation model [11] is of particular interest since it stems from the famous Winkler-Pasternak two-parameter foundation model [9] for which several applications and solutions have been available. In the Kerr-type foundation model, the foundation medium is visualized as consisted of lower and upper spring beds sandwiching an incompressible shear layer. Three parameters characterizing the Kerr-type foundation model are the lower and upper spring moduli and the shear-layer section modulus. It is noted that the interactive foundation force of the Kerr-type foundation model can be written in the same mathematical form as obtained with the simplified continuum models of Reissner [13], Horvath [14, 15], and Worku [16]. Synthetic and hierarchical correlations between several mechanical subgrade models and simplified continuum models are comprehensively presented in Horvath [17] and Horvath and Colasanti [18]. The pros and cons of each model are summarized in Horvath and Colasanti [18].
Even though the Kerr-type foundation model was developed since the mid-sixties, there have been only a limited number of researchers studying the problem of beams resting on Kerr-type foundation. Avramidis and Morfidis [19] used the principle of stationary potential energy to derive the governing differential equilibrium equations of the beam-foundation system and its essential boundary conditions. Subsequently, Morfidis [20, 21] derived the exact beam-foundation stiffness matrix based on the exact solution of the governing differential equilibrium equations for static and dynamic analyses, respectively, and calibrated the foundation parameters with the analysis results obtained with high fidelity 2D finite element models. The problem of beams resting on tensionless Kerr-type foundation was also investigated by Zhang [22] and Sapountzakis and Kampitsis [23]. Wang and Zeng [24] used the Kerr-type foundation model to study the interface stress between piezoelectric patches and host structures.
It is worth mentioning that a series of research papers on the so-called “modified Kerr-Reissner hybrid” foundation model have been presented by Horvath and Colasanti [18] and Colasanti and Horvath [25]. This foundation model is also regarded as the three-parameter foundation model and is formulated based on the combination of the modified Kerr-type foundation model with the Reissner simplified continuum subgrade model. In the modified Kerr-type foundation model, a shear layer is replaced by a tensioned membrane for the sake of modeling ease. The modified Kerr-Reissner hybrid foundation model is attractive particularly to practicing geotechnical engineers since it combines the advantages of both mechanical subgrade model and simplified continuum model as comprehensively discussed in Colasanti and Horvath [25]. Horvath and Colasanti [18] discuss the detailed derivation of this foundation model; Colasanti and Horvath [25] illustrate the modeling approach of this foundation model using commercially available structural analysis software. Subsequently, the modified Kerr-Reissner hybrid foundation model is applied to the planar geosynthetics used for tensile earth reinforcement under vertical loads [26].
In this paper, the virtual force principle is employed to reveal the governing differential compatibility equations of the beam-Kerr foundation system, as well as its natural boundary conditions. Thus, this paper can naturally be considered as a companion paper to the earlier work on the beam-Kerr foundation system by Avramidis and Morfidis [19] and Morfidis [20]. Unlike the structural component-based approach used by Horvath and Colasanti [18] and Colasanti and Horvath [25], all system components can be combined effectively into a single element, thus rendering the proposed model more attractive and unique from the theoretical and modeling point of view. The exact beam-Kerr foundation stiffness matrix is alternatively derived based on the exact beam-Kerr foundation flexibility matrix. The exact force interpolation functions of the beam-Kerr foundation system are at the core of the derivation of the exact element flexibility matrix. The governing differential equilibrium equations and constitutive relations of the beam on Kerr-type foundation are first presented. Next, the sixth-order governing differential compatibility equations, as well as the associated end-boundary compatibility conditions, are derived based on the virtual force principle. The exact force interpolation functions of the beam-foundation system are derived from the analytical solution of the governing differential compatibility equations of the problem. The matrix virtual force equation is employed to obtain the exact element flexibility matrix using the exact force interpolation functions. It is worth mentioning that the element flexibility matrix presented in this paper is different from that presented in Limkatanyu and Spacone [27] in that the foundation force distribution in Limkatanyu and Spacone [27] has to be assumed, thus resulting in the approximate moment interpolation functions and the approximate element flexibility matrix. The exact element stiffness matrix can be obtained directly from the exact element flexibility matrix following the natural approach [28]. It is noted that the natural approach had been used with successes in deriving the exact element stiffness matrices for beams on Winkler foundation [29] as well as beams on Winkler-Pasternak foundation [30]. It is also imperative to emphasize that, in the proposed model, the applied distributed load does not influence the exact force interpolation functions as long as it varies uniformly along the whole length of the beam. This finding renders the proposed flexibility-based model attractive since the analytical solution to the governing differential compatibility equation requires only the homogeneous part. Unfortunately, this beneficial effect is not available in the exact stiffness-based model presented by Avramidis and Morfidis [19] and Morfidis [20] since the analytical solution to the governing differential equilibrium equation requires both homogeneous and particular part with the presence of the applied distributed load. Therefore, the derivation of the exact displacement interpolation functions becomes more involved. A brief discussion on the efficient way to account for the extended-foundation effect is also introduced. All symbolic calculations throughout this paper are performed using the computer software Mathematica [31], and the resulting beam-foundation model is implemented in the general-purpose finite element platform FEAP [32]. A numerical example is used to verify the accuracy and the efficiency of the natural beam element on Kerr-type foundation and to show a more realistic distribution of interactive foundation force. A 2D finite element package VisualFEA [33] is also used to analyze this numerical example for comparison purpose.
2. Governing Equations of Beams on Kerr-Type Foundation
2.1. Differential Equilibrium Equations: Direct Approach
A beam-Kerr foundation system is shown in Figure 1(a) and comprises a beam, the upper and lower springs, and an intermediate shear layer. The governing differential equilibrium equations of the system are derived in a direct manner as follows.
(a)
(b)
(c)
A differential segment taken from the beam on Kerr-type foundation is shown in Figure 1(b). The vertical equilibrium of the infinitesimal beam segment is written as where is the beam-section shear force; is the transverse distributed load; and is the interactive force in the upper spring and acts at the bottom face of the beam. Considering the moment equilibrium, this yields where is the beam-section bending moment. Following the Euler-Bernoulli beam theory, only flexural contributions are considered in the paper. Enforcing the beam shear equilibrium of (2), (1) and (2) can be combined into a single relation; thus,
A differential segment taken from the shear layer resting on the lower foundation springs is shown in Figure 1(c). The vertical equilibrium of the infinitesimal shear-layer segment can be written as where is the shear-layer section shear force and is the interactive force in a lower spring. Equations (3) and (4) form a set of governing differential equilibrium equations of the system and are coupled through the upper-spring interactive force .
It is noteworthy to remark that this system is internally statically indeterminate and the internal forces cannot be determined simply by equilibrium conditions since there are 4 internal force unknown fields, , , , and , at any system section while only two equilibrium equations are available.
2.2. Deformation-Force Relations
The system sectional deformations can be related to their conjugate-work forces as follows: where is the beam-section curvature; is the shear-layer section shear strain; is the lower-spring deformation; is the upper-spring deformation; IE is the flexural rigidity; is the shear-layer section modulus; is the lower-spring modulus; and is the upper-spring modulus. Following the comprehensive work by Worku [16], the three foundation parameters (, , and ) can be related to the elastic modulus, Poisson ratio, and depth of the soil continuum underneath the beam.
2.3. Differential Compatibility Equations and End Compatibility Conditions: The Virtual Force Principle
The virtual force equation is an integral expression of the system compatibility equations and can be expressed in the general form as where is the system total complementary virtual work; is the system internal complementary virtual work; and is the system external complementary virtual work.
In the case of the beam-Kerr foundation system, and can be expressed as where is the beam vertical displacement; the vector contains shear forces and moments acting at beam ends and shear forces acting at the shear-layer ends; and the vector contains their conjugate-work displacements and rotations at the beam ends and displacements at the shear-layer ends. At the moment, external force quantity, , is arbitrarily chosen to be zero. Thus, (6) becomes
Eliminating the internal deformation fields through the deformation-force relations, (8) can be written as
The lower and upper spring forces and their virtual counterparts can be eliminated in (9) through the governing differential equilibrium equations of (3) and (4). Therefore, the system virtual force equation can be written as
In order to move all differential operators to the bending moment and the shear-layer section shear force , integration by parts is applied once to the forth term and twice to the fifth and sixth terms of (10), respectively, hence resulting in the following expression:
Considering the shear-force definition of (2) and following the Cartesian sign convention, (11) can be written as
Due to the arbitrariness of and , the governing differential compatibility equations of the beam and shear-layer components are obtained; thus
Equations (13) and (14) form a set of governing differential compatibility equations of the system. It is noted that the compatibility equations of the lower and upper spring deformations are not involved in the virtual force equation since their conjugate-work forces are eliminated in (10). However, they can be obtained simply by considering the geometrical deformations of the lower and the upper springs in Figure 1(a) as where is the shear-layer vertical displacement. Considering the deformation-force relations of (5) and enforcing the equilibrium equations of (3) and (4) and compatibility equations of (15), (13), and (14) are reduced to
It now becomes clear that (13) and (14) simply state the definitions of the beam section curvature and shear-layer section shear strain, respectively.
To make use of (13) and (14), there is a need to establish the relation between and . This could be accomplished by recalling the governing differential equation of the foundation surface subjected to a continuously distributed load as given by Kerr [11]:
Enforcing the equilibrium equations of (3) and (4) as well as compatibility equations of (15), recalling the curvature definition of (16), and considering the deformation-force relations of (5), (18) relates the first derivative of the shear-layer section shear force to the bending moment and its fourth-order derivative as
Differentiating (19) twice and substituting into (13) yield the following sixth-order differential equation: where ; and .
It is noted that when the upper-spring modulus approaches infinite, (20) is reduced to a fourth-order governing differential compatibility equation of the beam on Winkler-Pasternak foundation as given by Limkatanyu et al. [30] and when the shear-layer section modulus is equal to zero, (20) becomes a fourth-order governing differential compatibility equation of the beam on Winkler foundation as given by Limkatanyu et al. [29]. Furthermore, when compared to the governing differential equation derived by Avramidis and Morfidis [19] using the principle of stationary potential energy (equivalent to the virtual displacement principle), it becomes clear that (20) and the one derived by Avramidis and Morfidis [19] are dual. This illustrates the dualism of the virtual displacement and virtual force principles.
The end-boundary compatibility conditions are obtained by accounting for the arbitrariness of in (12) as
It is observed that the homogeneous and particular contributions to the end displacements are clearly separated in (21). This observation is unique to the proposed formulation and very useful in determining the equivalent fixed-end force vector due to the element load .
3. “Exact” Element Stiffness Matrix: Natural Approach
In this paper, the “exact” element stiffness matrix is derived simply by inverting the “exact” element flexibility equation. This is feasible since the system does not experience any rigid-body motion (neither rigid-body translation nor rigid-body rotation) due to the presence of supporting foundation. Thus, the exact element flexibility matrix is at the core of the element formulation and requires the “exact” moment interpolation functions. The analytical solution to the sixth-order governing differential compatibility equation of (20) is central to obtain the exact moment interpolation functions. For the sake of simplicity, the applied distributed load is assumed to be uniform along the whole length of the beam. Thus, only homogeneous solution is required to derive the exact moment interpolation functions. This merit comes from the fact that the force terms on the right-hand side of (20) disappear as long as varies uniformly along the whole length of the beam, thus rendering the proposed flexibility-based model attractive and unique.
Thanks to the comprehensive investigation performed by Morfidis [34] and Avramidis and Morfidis [19] on all possible solutions to the similar sixth-order differential equation, the general solution of (20) can be written as where , , , , , and are real functions and their forms depend on the values of the system parameters (, , and ) as given in Appendix A; and , , , , , and are constants of integration to be determined by imposing force boundary conditions. These six force boundary conditions are
The first four boundary conditions on the beam ends can be imposed directly while the last two on the shear-layer ends cannot be enforced at the first stage since the beam-section bending moment is the only variable field in the governing differential compatibility equation of (20). This difficulty can be overcome by establishing the relation between and .
Recalling the shear-layer compatibility condition of (14) and the relation of (19), the shear-layer shear force can be expressed in terms of the beam-section bending moment and its derivatives as where , , and .
By imposing force boundary conditions of (23), the moment interpolation relation can be expressed as where is an array containing the moment interpolation functions. Imposing the relation of (24) and differential equilibrium equations of (3) and (4), the shear-layer shear force , the lower-spring force , and the upper-spring force can be expressed in terms of as where is an array containing the shear-layer shear-force interpolation functions; is an array containing the lower-spring force interpolation functions; and is an array containing the upper-spring force interpolation functions.
Applying the virtual force expression of (9), substituting (25)-(26), and accounting for the arbitrariness of yield the following element flexibility equation: where is the element flexibility matrix, defined as where , , , and are the beam, the shear-layer, the lower-spring, and the upper-spring contributions to the element flexibility matrix, respectively:
It is noted that element end displacements due to the applied load is supplemented into (27). In the case of linear variation of , can be written in a simple expression as given in Appendix B.
Based on the element flexibility expression of (27), the element stiffness equation can be written as where the complete element stiffness matrix is and the fixed-end force vector due to is simply computed as . It is worthwhile to note that the subscript stands for “natural.” This is due to the fact that the approach employed herein to obtain the element stiffness matrix is known as the natural approach [28]. The configuration of the natural beam element on Kerr-type foundation is shown in Figure 2.
Unlike the stiffness-based formulation, the displacement fields cannot be computed directly since no displacement interpolation function is available in the element formulation. However, the following compatibilities can be used to retrieve the vertical displacement and rotational fields of the beam component and the vertical displacement field of the shear-layer component once the internal force distributions are obtained:
4. Restrained Effects of Extended Kerr-Type Foundation on the Beam End
When the foundation on either end of the beam is infinitely extended, appropriate modeling of the beam-end condition is deemed essential to account for the foundation continuity [35]. One efficient way to consider this end effect is to place a vertical spring with a stiffness of at the associated beam end as suggested by Eisenberger and Bielak [36]. A detailed derivation of this stiffness value can be found in Alemdar and Gülkan [35] and Colasanti and Horvath [25]. For the case of finitely extended foundation, a virtual beam-foundation element with a small value of the flexural rigidity and large value of the upper-spring modulus can be assumed beyond its physical end to account for the existence of the extended foundation.
5. Numerical Example
A free-free beam on an infinitely long Kerr-type foundation subjected to various loads along its length is shown in Figure 3. This beam-foundation system was also studied by Morfidis [20] and is used in this study to verify the accuracy and to show the efficiency of the natural beam element on Kerr-type foundation. The flexural stiffness and width of the beam are ?kN-m2 and 1?m, respectively. The elastic soil mass underneath the beam is 10?m depth and is assumed to be loose sand with elastic modulus and Poisson ratio . Following the modified Kerr-Reissner model [18, 19], the lower-spring and shear-layer section moduli are found to be and , respectively. As suggested by Avramidis and Morfidis [19], the upper-spring modulus is related to the lower-spring modulus as where is a factor expressing the relative stiffness of the upper and the lower springs. Following comprehensive correlation studies between the three-parameter foundation and the high fidelity 2D finite element models by Avramidis and Morfidis [19], the optimal values of are suggested depending on the system parameters. In this example, the value of is equal to 7 which is the optimal value for soft soils [19]. Thus, the value of is equal to . To account for the effect of infinitely extended foundation beyond both beam ends, a vertical spring with stiffness is placed at each end as shown in Figure 3. Seven natural beam-foundation elements (elements AB, BC, CD, DE, EF, FG, and GH) are used to discretize the system, thus resulting in twenty-four nodal unknowns.
The beam-foundation system of Figure 3 is also analyzed by the 2D finite element model [33]. Figure 4 shows the 2D finite element mesh of the beam-foundation system of Figure 3. The virtual soil mass of the length ?m is assumed beyond each beam end to account for the existence of the infinitely long foundation. The soil mass is discretized into 950 rectangular plane-strain elements while the beam is modeled with 35 conventional beam elements. In order to ensure the sufficiency of the model discretization, a finer finite element mesh was used but yielded the same analysis results.
Figure 5 shows the obtained beam vertical displacement, beam rotation and shear-layer vertical displacement diagrams while Figure 6 shows the obtained beam shear force, beam moment, and shear-layer force diagrams. The exact displacement-based responses given by Avramidis and Morfidis [19] and Morfidis [20] as well as the responses obtained with 2D finite element analysis (2D FEM model) are also superimposed for comparison on the respective diagrams. Clearly, the natural beam-foundation model is capable of representing the exact displacement and force responses using only one element for the beam span. Winkler-foundation responses obtained with the model by Limkatanyu et al. [29] are also presented in the same respective diagrams. The results presented in Figure 5 indicate that the Kerr-type foundation model plays a role in reducing the vertical displacement and rotation of the beam, thanks to the coupling between the Winkler-foundation springs. This coupling effect renders the Kerr-type foundation model with the ability to resemble the displacement and rotation diagrams obtained with the 2D FEM model. When compared to the Winkler-foundation model, Figure 6 shows that the Kerr-type foundation model affects the bending moment response more than the shear-force response along the beam. Furthermore, the bending moment response obtained with the Kerr-type foundation model is closer to that obtained with the 2D FEM model when compared to the Winkler-foundation model. It should be kept in mind that a complete comparison between the proposed model and the more sophisticated finite element model is not to be expected. This is due to the fact that a full compatibility at the beam-soil interface is assumed in the finite element model while only the vertical displacement compatibility is enforced in the proposed model [37]. In this example, introducing the more refined foundation model generally results in reducing the negative moment (concave) but slightly increasing the positive moment (convex).
(a)
(b)
(c)
(a)
(b)
(c)
Figure 7 shows the upper and the lower spring force diagrams. Obviously, the proposed beam-foundation element is capable of representing the exact foundation-spring force distributions along the beam length. Figure 7(a) compares the interactive foundation force acting at the bottom face of the beam obtained with the Winkler and Kerr-type foundation models. Evidently, the distribution characteristics of these two foundation interactive forces are distinctively different. The interactive foundation force distribution obtained with the Kerr-type foundation model corresponds well to the observation made by Foppl [12] that there exists a certain degree of discontinuities in system responses between the loaded and the unloaded regions of the actual beam-foundation system. This feature is unique to the Kerr-type foundation model and clearly indicated by incompatibilities between the upper and the lower foundation spring forces at both beam ends ( and ), thus resulting in accurately representing the peripheral reactions of the beam ends [6].
(a)
(b)
6. Summary and Conclusions
The “natural” element stiffness matrix and the fixed-end force vector for a beam on elastic foundation subjected to a uniformly distributed load are derived in this paper. The Kerr-type foundation model is employed to model the underlying foundation continua, thus taking into account the shear coupling between the individual Winkler-spring components through the shear-layer component and determining the level of vertical-displacement continuity at the boundaries between the loaded and the unloaded soil surfaces. This feature is unique to the Kerr-type foundation model. The element flexibility matrix forms the core of the natural element stiffness matrix and is derived based on the virtual force principle using the “exact” force interpolation functions. The exact force interpolation functions are obtained by solving analytically the sixth-order governing differential compatibility equation. Compared to the stiffness-based models published in the literatures, the effect of the applied element load can readily be included in the proposed formulation. One numerical example is employed to verify the accuracy and efficiency of the natural beam-foundation model. This numerical example shows that the natural beam-foundation element is capable of giving exact system responses. Therefore, the exactness of the proposed element obviates the requirement for discretizing the beam into several elements between loading points. The number of elements needed in the analysis of a beam-foundation system is largely dictated by the convenient way of representing loadings (concentrated or distributed loads). Furthermore, the Kerr-type foundation model results in more realistic interactive foundation forces as compared to the Winkler foundation model. A 2D finite element model is also used to confirm the superiority of the proposed model. The next step forward in this research direction is to account for system nonlinearities and to apply the resulting model to practical soil-structure interaction problems. This could be accomplished by first following the evolution of beam and foundation mechanical parameters and then updating the force interpolation functions accordingly. Another interesting topic worth investigating in future works is the derivation of consistent mass and geometric stiffness matrices based on the force interpolation functions.
Appendices
A. Homogenous Solution to the Sixth-Order Governing Differential Compatibility Equation (20)
The homogeneous form of (20) can be written as
For simplicity, the following auxiliary variables are introduced instead of terms of system parameters , , and :
There are many solution types to (A.1) depending on the sign of the auxiliary parameter . Thanks to the thorough investigations of all possible solution cases, Avramidis and Morfidis [19] and Morfidis [20] suggest the most common solution case corresponding to the positive sign of the auxiliary parameter as follows.
Solution Case (when ). The homogeneous solution of (20) can be written as where
B. Nodal Displacements due to
The nodal displacements due to the uniformly distributed load may be written as
Acknowledgments
This study was partially supported by the Thai Ministry of University Affairs (MUA), by the Thailand Research Fund (TRF) under Grant MRG4680109 and Grant RSA5480001, by STREAM Research Group under Grant ENG-51-2-7-11-022-S, Faculty of Engineering, Prince of Songkla University, and by the National Research Foundation of Korea (NRF) Grant funded by the Korean government (MEST) (NRF-2011-0028531). Any opinions expressed in this paper are those of the authors and do not reflect the views of the sponsoring agencies. Special thanks are due to the senior lecturer Mr. Wiwat Sutiwipakorn for reviewing and correcting the English of this paper and Mr. Paitoon Ponbunyanon and Mr. Attpon Sangkeaw for preparing the numerical example. In addition, the authors would also like to thank three anonymous reviewers for their valuable and constructive comments.