Journal of Applied Mathematics
Volume 2015 (2015), Article ID 420649, 12 pages
http://dx.doi.org/10.1155/2015/420649

Research Article

Solution of Contact Problems for Nonlinear Gao Beam and Obstacle

J. Machalová and H. Netuka

Faculty of Science, Palacký University in Olomouc, 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic

Received 21 May 2015; Revised 24 July 2015; Accepted 14 August 2015

Academic Editor: Hui-Shen Shen

Copyright © 2015 J. Machalová and H. Netuka. 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.

Contact problem for a large deformed beam with an elastic obstacle is formulated, analyzed, and numerically solved. The beam model is governed by a nonlinear fourth-order differential equation developed by Gao, while the obstacle is considered as the elastic foundation of Winkler’s type in some distance under the beam. The problem is static without a friction and modeled either using Signorini conditions or by means of normal compliance contact conditions. The problems are then reformulated as optimal control problems which is useful both for theoretical aspects and for solution methods. Discretization is based on using the mixed finite element method with independent discretization and interpolations for foundation and beam elements. Numerical examples demonstrate usefulness of the presented solution method. Results for the nonlinear Gao beam are compared with results for the classical Euler-Bernoulli beam model.

1. Introduction

Contact problems belong to the most important industrial applications and contact problems for beams have their own significant position among them. The Euler-Bernoulli beam is the most popular model used in engineering applications. This model is linear and its validity is limited only to relative small deflections. If we suppose large deformations, we have to switch the mathematical model to a nonlinear one. One of the best nonlinear beam models was developed by Gao in [1].

In this paper we are going to deal with a beam and an elastic obstacle which are in possible contact. For convenience we start our considerations with the classical Euler-Bernoulli beam model and the Winkler foundation. Nowadays it is well known that the contact of the elastic bodies is usually modeled using the Signorini conditions and followed by variational inequalities (see, e.g., [2, 3]).

But this approach is necessary only for the case when the obstacle is rigid as it has already been published in [4]. Using the so-called normal compliance condition (for more details see, e.g., [5] or [6]) for a deformable foundation we get description in the form of variational equation. This equation is of course nonlinear.

Here we can find as very useful the so-called control variational method described in the papers [4, 7, 8] and thoroughly analyzed in [9, Chapter 6]. The functional of total potential energy is transformed in such a way that we are able to formulate an optimal control problem governed by the beam equation. Such problem was separately studied in [10], and much more general problems were considered in the excellent monographs [11, 12]. The nonlinear terms from our variational formulations are fully included in the control variable. Afterwards we solve the resulting optimal control problem to obtain a solution of our initial contact problem.

2. Contact Problems for Classical Beam

2.1. Classical Formulations

Our considerations will start by recalling the problem with a beam resting on an elastic foundation (see Figure 1). For convenience we will deal initially with the Euler-Bernoulli model on Winkler foundation. The equation reads as follows:

where is the deflection of the beam, is the length of the beam, is Young’s elastic modulus, is the area moment of inertia, is the applied vertical load, and is the foundation modulus. This situation in fact represents the case when the foundation is firmly connected to the beam. We can speak about the bilateral foundation (according to [13] or [5]). If the beam only rests on the foundation, the above equation must be modified, especially the second term in (1), and to distinguish such cases it is sometimes spoken about the unilateral foundation.

Figure 1: Beam resting on elastic foundation.

First let us define the positive and negative part of the function as follows:

It can be easily seen that

Now if at some point , we have no contact between the beam and the foundation and consequently the second term in (1) is zero. If conversely , then we get a penetration into the foundation. Hence substituting (2) and (3) into (1) we obtain the equation

Adding the boundary conditions

we obtain the problem which was studied in [14] and the analogous one in the papers [15, 16] for a nonlinear beam.

Another situation occurs when the beam is resting on a rigid foundation. Then we must instead of (1) take into account the following equation:

where represents the contact force which depends on unknown deflection. As the only possibility is in , we can consider two cases. If for some , then the beam has contact with the foundation and a contact force appears according to the third Newton law. If for some , then we have no contact at this point and therefore . From these considerations we are able to formulate special conditions for this case as follows:

in (). These relations are well known from the literature as the Signorini conditions and must be added to the final mathematical formulation of the problem with rigid foundation (see, e.g., [3, 13]).

Moreover, looking back on (4) we can see that the problem with a unilateral elastic foundation could be expressed in the form (6) too. And by comparing both equations we get

It is quite useful result for our next considerations and it is usually referred to as the so-called normal compliance condition.

Next we add some gap between the beam and the foundation and a support on the right end of the beam (see Figure 2). Generally the gap is a function describing an obstacle, but for simplicity we will consider here only a constant value .

Figure 2: Beam above elastic foundation.

Regardless of whether the foundation is rigid or elastic, the mathematical formulation for this case must be based on the same equation (6). At first we make some observations for an elastic foundation similarly as above. If for some , then we have contact and consequently a contact force as well. If for some , then we have no contact at this point and therefore . Hence we obtain the contact problem

where by analogy with (8) we have

This can be interpreted again as the normal compliance condition which generally describes a reactive normal force or pressure that depends on the penetration into the foundation.

For the case of rigid foundation we must add to (9)–(11) the slightly modified Signorini conditions (7)

in .

We can observe that both problems are quite similar and letting the gap equal zero in the second problem (9)–(11) we obtain the mathematical formulation for the first one.

Finally, adding the support to the right end of the beam, we limit our considerations without loss of generality to the case when the resultant of acting loads at is negative; that is, it is headed downward, because for such cases we can obtain a contact between the beam and the support.

2.2. Weak and Variational Formulations

After the classical formulation we now proceed to the variational one which is more suitable for a numerical solution. First we must define a suitable space for test functions. According to the boundary conditions it will be

where is the Sobolev space (see, e.g., [2, 3]). Next let us introduce for any the notations

Evidently represents a bilinear form on and the inner product in the Lebesgue space .

The well-known fact consists in the fact that contact problems can be expressed by means of variational inequalities (see, e.g., [3, 13]). For the case with a rigid foundation we introduce the set of admissible deflections (or test functions) as

Using scalar multiplication with functions we obtain from (6)

From (13) and (16) it follows that

and hence the contact problem (9)–(11) together with the Signorini conditions (13) can be now formulated as the following variational inequality:

 Find a function   such that

For the problem with an elastic deformable foundation, the formulation is obtained by taking the -inner product of (9) with test functions . After that we obtain

It is quite interesting that for the elastic foundation we received the variational equality instead of an inequality as it is usual in contact problems formulations. Equation (20) may be rearranged using the relation (12) for contact forces. Applying the notation (15) we get the weak formulation of our problem:

 Find a function such that

This is of course nonlinear problem as it shows the second term in (21).

Next let us associate the energy functional

with the Euler-Bernoulli beam. Then the problem (19) can be equivalently expressed using minimization problem:

 Find such that

We are able also to construct the functional of total potential energy for the problem (21). Let us introduce the energy of foundation reaction forces by

Then we have

and the minimization problem reads as follows:

 Find such that

This problem is convex and differentiable in the sense of Gâteaux. This implies the equivalence (26) to the previous formulation (21).

3. Nonlinear Gao Beam

In the previous section we were concerned for simplicity with the standard Euler-Bernoulli beam model. But our main interest in this paper will be to study contact for the nonlinear beam model which was proposed by Gao in [1, 17]. We will mention this model here only briefly and only the most simple version (with zero axial force) given by a fourth-order nonlinear differential equation

The beam has a constant stiffness given by an elastic modulus and a constant area moment of inertia . Its length is and thickness is measured from the -axis, and therefore the full thickness will be . The width of the beam is considered as a unit. The transverse load is denoted by and describes the deflection of the beam at a position . Equation (27) has the following built-in relationships:

where denotes the Poisson ratio. Boundary conditions remain the same as before, that is, (10)-(11). For the load we will assume the same as in the previous section for load ; that is, the resultant at must be negative.

Of course, the considerations concerning the contact with foundation that we have made before remain valid. Hence the classical formulation for the contact problem with a rigid foundation reads as follows:

in , and for an elastic deformable foundation we have

with .

Variational formulation for this problem we can get by a similar way as above. Let us define the functionals of total potential energy for our problems by adding the nonlinear term to the functional or from the previous section with substituting for and for as

This is a nonlinear convex functional. Hence we can again formulate the minimization problem for a deformable foundation as follows:

 Find such that

and this problem is due to its convexity equivalent to the following nonlinear variational equation:

 Find a function such that

where we denoted

and means the Gâteaux differential of at the point in the direction . We of course substituted for in as well.

Analogously we get for a rigid foundation the next two problems:

 Find   such that

where (again substituting for )

and (equivalently)

 Find a function such that

4. Optimal Control Problem

Our approach will be based on an optimal problem formulation and solution. The idea is due to the method applied in the papers [4], [7], or [8]; the authors call it the control variational method.

First of all, let us recall briefly the concept of the optimal control of elliptic equations by the right hand side (see, e.g., [11, 12]). Let be a subset of a space of controls , and let be a linear continuous operator from into , where is the dual space to a Hilbert space . For any given and any we will define the problem:

 Find a function   such that

Subset is called the set of admissible controls and (38) is the classical formulation of the state problem. Such type of problems belongs to problems with distributed control which means that acts in the whole interval .

Next we will proceed to the variational formulation for the state problem which will be in most cases better suited than the previous one. At the same time we will specify more precisely its setting. Let be a nonempty bounded closed convex subset of a reflexive Banach space , and let be a linear continuous mapping from to . For any we will want to

 Find a function such that

where denotes the duality pairing between and .

Let a cost (or criterion) functional be given and let it be weakly lower semicontinuous; that is,

Now let us define the resulting optimal control problem:

 Find a function such that

where solves the state problem (39).

It is important that we are able to prove the following proposition (using arguments from [11] or [12]).

Theorem 1. The problem (41) has under the abovementioned assumptions at least one solution.

Proof. Let

that is, is a minimizing sequence of (41). As is bounded and is reflexive, there exists a subsequence such that

At the same time , which are the solution of

are bounded as well and we may assume that

Passing to the limit with  + in (44) we see that

that is, . On the other hand from (40)

that is, is a solution of (41).

Finally, we will deal with optimality conditions for our optimal control problem. In the text below we will use for simplicity the abbreviation

with the linear continuous operator from to . The key idea lies in the fact that the problem (41) can be considered as a constrained minimization problem. For this purpose let us define the set

It is evident that

At the same time we can see that

Let us denote

Then it is easy to see that

is the indicator function of the set . This implies

and consequently we can define the Lagrangian with a Lagrange multiplier

Let us consider a saddle point of . The following proposition is a known characterization of saddle points.

Lemma 2 (see, e.g., [18, Chapter VI, Proposition ]). An element is a saddle point of a Lagrangian if and only if

Taking into account the definition (55) and Lemma 2 it can be seen that

From here it is evident that is a solution of (41).

Next let us suppose that is differentiable. Then from the saddle point definition it follows that

From the first condition (58) we get

whereas the last condition (60) gives

which means that . From (61) using definition of the adjoint operator to the operator we have

which is the so-called adjoint equation. We can see that its solution , called the adjoint state, is nothing else but the optimal value of the Lagrange multiplier associated with the constraint given by the set . Finally the second condition (59) can be according to (52) written as

where is the adjoint operator to . The three relations (62)–(64) form the requested optimality conditions.

Later it will be important to know the gradient of the cost functional . To this purpose let us denote

and analyze properties of the mapping . It is known (e.g., from the monographs [11, 12]) that

where is a solution to the adjoint equation, and it implies

5. Transformation to Optimal Control Problems

5.1. Problem with a Deformable Foundation

Now let us return to our contact problems (32) and (35). Our next steps will be transformations of the energy functionals and . Main attention will be focused on the first case as the second one with a rigid foundation is for many years well known in the literature (e.g., [2, 13]).

First let us perform some rearrangement of its part given by (22). For this purpose let the state problem be (39) with operator as identity; that is,

 Find   such that

Let be the (unique) solution of the state problem for the control variable ; that is, fulfills

Then substituting for and using (68) and (69) we find that

As the last term in (70) is in fact a constant, we can omit it from the viewpoint of minimization and define the new functional as

We still have not defined the set . From a physical point of view it seems to be a good choice followed by

for some positive constant (depending on the beam compliance) as we do not want to break the beam.

Thus the cost functional will be therefore the following conversion of :

Finally we get the first optimal control problem:

 Find such that

where solves the state problem (68).

We can state the following proposition about a relation between the problems (74) and (32).

Theorem 3. Let solve the optimal control problem (74). Then is a solution of the corresponding variational problem (32) providing that the constant in (72) is big enough.

Proof. Let be an optimal pair of the problem (74). Let us consider admissible pairs of the form for any , , and any satisfying the state equation (68)

As

we get after substitution that must fulfill

or

From the optimality of we have

Substituting the above mentioned expressions for and we get after some manipulations

Dividing this inequality by and letting + we obtain (using (76), (77), and (34)) that

for any satisfying (78). As we can take positive as well as negative, this relation gives us the desired equality (33) which is equivalent to the minimization problem (32).
Finally let us notice that the assumption regarding the constant from the definition of means that and were inner points of the set .

Theorem 4. The variational problem (32) has exactly one solution.

Proof. The proof follows immediately from the fact that the functional is strictly convex on . This property can be easily deduced from the inequality

where denotes the -norm; that is, .

As a consequence of Theorems 1, 3, and 4 we obtained the following proposition.

Corollary 5. The optimal control problem (74) has exactly one solution under the assumption that the constant in (72) is big enough.

Proof. As the functional is weakly lower semicontinuous and the set in (74) is a bounded subset of , we get from Theorem 1 that (74) has at least one solution. For every such solution its first component solves the associated variational problem (32) by Theorem 3. But (32) has exactly one solution by Theorem 4, whereas the mapping is linear according to the state equation (68). This completes the proof.

5.2. Problem with a Rigid Foundation

The state problem remains formally the same as before, that is, (68). But now we must respect the contact conditions (13), especially . Of course, a solution to the state equation generally does not satisfy such a condition for arbitrary . In order to guarantee it, the set has to be properly redefined.

Using notation (48) we have from (38)

and it follows that

Hence let us define

Obviously for any we have a unique solution of the problem (68) and it fulfills ; that is, it lies in .

The cost functional will be the following conversion of :

Finally we get the second optimal control problem:

 Find   such that

where solves the state equation (68).

Now it is important to establish a proposition about the equivalence of the problems (87) and (35).

Theorem 6. Let solve the optimal control problem (87). Then is a solution of the corresponding variational problem (35) providing that the constant in (85) is big enough.

Proof. Let be an optimal pair of the problem (87). Let us consider admissible pairs of the form for any , , and any satisfying the state equation (68).
Using similar arguments as in the previous proof of Theorem 3 we come to the following inequality:

As , we obtain the final result in the form

which is the requested problem (37).

Finally, we may use almost the same arguments as in the previous paragraph to obtain the following proposition.

Corollary 7. The optimal control problem (87) has exactly one solution.

Remark 8 (construction of the set ). Determination of the set (72) does not make any trouble unlike the case with somewhat abstract definition (85). But these difficulties can be solved by using a penalization, especially by means of barrier function method. Hence we add to the functional the expression

and take the limit for (see, e.g., [2, Chapter 2]).

6. Numerical Realization and Algorithms

Next some approximations of our problems are needed since exact solutions are not available. We start with the general problem (41). Instead of the state problem (39) its full finite element discretization must be taken into account, that is, not only discretization of function spaces but also approximations of the linear form and the mapping denoted by , respectively. Here denotes the discretization parameter which in the finite element interpretation is given as the maximal diameter over all elements from the current finite element mesh. For our purposes we do not need numerical computations of the bilinear form (we will discuss it a bit later).

Let and be finite-dimensional subspaces of and , respectively, and let be a nonempty closed, convex, and bounded subset of ; that is, we do not assume that necessarily .

The state problem (39) is now replaced by its finite element approximation:

 Find such that

where is the duality pairing between and .

Let be a cost functional which is lower semicontinuous on ; that is,

By the approximation of the problem (41) we will understand the problem:

 Find a function such that

where solves the discretized state problem (91). The bilinear form is continuous and -elliptic. If we are allowed to assume and , then the discretized state problem (91) has exactly one solution for any . Let moreover (92) hold, and then it is easy to prove that the optimization problem (93) has at least one solution .

Now we want to establish a relation between the continuous problem (41) and the discrete problem (93) as +. This analysis is important for computations but we do not intend to go into details here since they are somewhat comprehensive.

In this section we establish matrix formulation of the optimal control problem (74). To this purpose let discretization parameter be fixed, let two divisions of into and subintervals be done, and let

Furthermore let

that is, , are basis functions of , respectively. Then any function can be expressed as

Analogously we have

The Galerkin method is based on the following idea: we substitute expression (96) and (97) into (68) and then as test functions successively choose basis functions . This way we obtain for the given vector (representing a given control parameter value) the system of linear algebraic equations with the unknown

where we denoted

The matrix is the well-known stiffness matrix, and is the load vector.

The finite element method usually works with the so-called Courant basis functions which are characterized by the following “-property”:

Applying this to the functions from , we can see according to (96) and (97)

Next defining

it is possible to rewrite the system (98) in the new and more usual form

The previous procedure is quite general; hence for next calculations we should specify the finite element basis more precisely. The standard beam elements are chosen for the beam and it is well known that they are working with Hermitian shape functions. The system is created by linking of appropriate shape functions from adjacent subintervals. This leads to the following description of the finite element space:

where denotes th subinterval of and means the space of cubic polynomials defined on . Values of the bilinear form for constant values and are then done exactly; that is, no numerical integration is needed. For more details see, for example, [19, Chapter 5].

Construction of the set uses here the same system of subintervals as , that is, value , and works with linear shape functions defined on them. Therefore we have

Next let us formally define two isomorphisms

Then the discrete control variable belongs to the set

Similarly we are able to construct the algebraic representation of as follows:

where is a multivariable function. Then the optimization problem (93) for a fixed parameter can be equivalently expressed as follows:

 Find a vector such that

where solves the linear system (103). This formulation represents a nonlinear programming problem.

Remark 9 (construction of the set ). As the definition (107) is rather formal, some comment is needed especially for the case with rigid foundation. Most algorithms can keep the constraints like providing that an initial guess satisfies this condition. Such starting guess can be constructed in our problems rather easily. For instance, let the load be nonpositive function. Then we get as well and there are two possible cases. In the first one we have and then we can put . In the second case we have in some points . Then we can try for some value such that the total load causes the deflection . Usually just a few attempts are enough.

Remark 10 (solving the nonlinear programming problem (109)). There are two good ways of how to solve the problem (109). The first one consists in using some gradient method. For this purpose, of course, we need to compute gradient of the cost functional or at any point. It can be done by means of adjoint state as it was established in (66) and (67). Returning to the optimality condition (63) we can write the adjoint equation for in the weak form as follows:

while in the case of the equation remains almost the same; only the last term on the right hand side vanishes. Solving this equation we get the adjoint state and then the required gradient by (67)

Let us notice that, instead of using penalization techniques as it was mentioned in Remark 8, we can more preferably use interior-point methods (see, e.g., [20]) particularly for the contact with Signorini conditions.

7. Examples

This section presents four examples which are intended to demonstrate the above explained theory and methods. In addition to the Gao beam model we will present also results for the classical Euler-Bernoulli beam model for comparison. As the nonlinear beam model is tougher than the classical one, the upper point curves show the results for the Gao beam and lower point curves are the results for the Euler-Bernoulli beam. Beams have always the left end fixed and the right end simply supported, as it was considered before (see Figure 2).

The common data are as follows:

 , , , , and . Number of finite elements .

Example 1 (contact with normal compliance condition). Consider the following:

 Input data: , , , and . Significant results: see the left half of Figure 3. E-B beam: maximal deflection , contact zone: . Gao beam: maximal deflection , contact zone: .

Figure 3: Results of Examples 1 and 2.

Example 2 (contact with normal compliance condition). Consider the following:

 Smaller load than in Example 1. Input data: , , , and . Significant results: see the right half of Figure 3. E-B beam: maximal deflection , contact zone: . Gao beam: maximal deflection , contact zone: .

Example 3 (contact with normal compliance condition). Consider the following:

 Larger foundation stiffness than in Example 1. Input data: , , , and . Significant results: see the left half of Figure 4. E-B beam: maximal deflection , contact zone: . Gao beam: maximal deflection , contact zone: .

Figure 4: Results of Examples 3 and 4.

Example 4 (contact with Signorini conditions). Consider the following:

 Input data: , , and . Significant results: see the right half of Figure 4. E-B beam: maximal deflection , contact zone: , Gao beam: maximal deflection , contact zone: .

Numerical computations were realized using MATLAB.

8. Conclusions

We presented here quite promising methods to solve contact problem for the Gao beam and an obstacle, especially a deformable or a rigid foundation. As a result we received an optimal control problem with following properties:

(i)State problem is always linear; nonlinearities are “hidden” in control variable .(ii)Objective functional has a “handy” form since the total potential energy of the system is convex.

Hence both optimal control problems are very suitable from a computational point of view.

Numerical experiments show us that the Gao beam is tougher than the classical beam. Computations were, of course, somewhat longer for the nonlinear beam. Another observation is in a good agreement with theoretical results (see [5, 6]): for increasing obstacle toughness we can observe a numerical convergence of results for normal compliance contact problems to results for the corresponding Signorini contact problem.

Conflict of Interests

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

References

D. Y. Gao, “Nonlinear elastic beam theory with application in contact problems and variational approaches,” Mechanics Research Communications, vol. 23, no. 1, pp. 11–17, 1996.

R. Glowinski, J.-L. Lions, and R. Tremolieres, Numerical Analysis of Variational Inequalities, North-Holland, Amsterdam, The Netherlands, 1981.

N. Kikuchi and J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, Pa, USA, 1988.

M. Sofonea and D. Tiba, “The control variational method for contact of Euler-Bernoulli beams,” Bulletin of the Transilvania University of Braşov, Series III, Mathematics, Informatics, Physics, vol. 2, pp. 127–136, 2009.

M. Shillor, M. Sofonea, and J. J. Telega, Models and Analysis of Quasistatic Contact: Variational Methods, Springer, 2004.

M. Sofonea and A. Matei, Mathematical Models in Contact Mechanics, Lecture Note Series 398, London Mathematical Society, Cambridge University Press, 2012.

M. Sofonea and D. Tiba, “The control variational method for elastic contact problems,” Annals of the Academy of Romanian Scientists. Series on Mathematics and its Applications, vol. 2, no. 1, pp. 99–122, 2010.

M. Barboteu, M. Sofonea, and D. Tiba, “The control variational method for beams in contact with deformable obstacles,” ZAMM, vol. 92, no. 1, pp. 25–40, 2012.

P. Neittaanmaki, J. Sprekels, and D. Tiba, Optimization of Elliptic Systems, Springer Monographs in Mathematics, Springer, New York, NY, USA, 2006.

J. Machalová and H. Netuka, “Optimal control of system governed by the Gao beam equation,” in Proceedings of the 10th AIMS Conference, Madrid, Spain, July 2014.

J.-L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Springer, Berlin, Germany, 1971.

F. Tröltzsch, Optimal Control of Partial Differential Equations: Theory, Methods and Applications, American Mathematical Society, Providence, RI, USA, 2010.

I. Hlaváček, J. Haslinger, J. Nečas, and J. Lovišek, Numerical Solution of Variational Inequalities, vol. 66 of Springer Series in Applied Mathematical Sciences, Springer, New York, NY, USA, 1988.

J. Machalová and H. Netuka, “Solving the beam bending problem with an unilateral Winkler foundation,” AIP Conference Proceedings, vol. 1389, pp. 1820–1824, 2011.

D. Y. Gao, J. Machalová, and H. Netuka, “Mixed finite element solutions to contact problems of nonlinear Gao beam on elastic foundation,” Nonlinear Analysis: Real World Applications, vol. 22, pp. 537–550, 2015.

J. Machalová and H. Netuka, “Bending of a nonlinear beam reposing on an unilateral foundation,” Applied and Computational Mechanics, vol. 5, pp. 45–54, 2011.

D. Y. Gao, “Finite deformation beam models and triality theory in dynamical post-buckling analysis,” International Journal of Non-Linear Mechanics, vol. 35, no. 1, pp. 103–131, 2000.

I. Ekeland and R. Témam, Convex Analysis and Variational Problems, SIAM, Philadelphia, Pa, USA, 1999.

J. N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill, New York, NY, USA, 3rd edition, 2006.

J. Nocedal and S. J. Wright, Numerical Optimization, Springer, Berlin, Germany, 2nd edition, 2006.