1. Introduction

Foundations of structures can be divided into shallow and deep foundations. In the modeling of shallow foundations, structural elements such as beams, arches, plates, and shells are supported by a continuous substrate. The primary difficulty in the analysis of the structure-substrate system lies in the determination of the contact pressure in the interface. Usually the substrate is modeled as an elastic foundation. Most constitutive models consider that the foundation reacts both under tension and compression. However, certain types of soil and most liquids only react under compression. In these circumstances, the structure may lose contact with the foundation in certain regions, leading to unexpected contact pressure concentration with a consequent variation of the structure’s internal forces. Problems where the structure can enter in or lose contact with other bodies, or even slide on its support, are usually found in the literature under the denomination of “unilateral contact problems” [1, 2].

The first step in obtaining the numerical solutions of contact problems involving contin-uous systems generally consists in reformulating the problems in approximation spaces. To this end, numerical techniques such as the Ritz, finite element, or boundary element methods are employed. After discretization, a proper methodology to adequately treat the unilateral contact constraints must be selected. This usually requires the problem to have a finite dimension. Among the options found in the literature, two are noteworthy.(i)Transformation of the contact problem into a minimization problem without restriction by applying usual formulations of the structural mechanics—differenti-able functional and bilateral constraints—to the case of unilateral contact constraints. There is no guarantee of convergence of these procedures, which by their nature are unavoidably incremental-iterative. These procedures, however, introduce no new concepts. So, existing codes for nonlinear analyses that can be adapted to this particular case, resulting in efficient computational time, granted that there is no change in the contact region between two load steps [3–7].(ii)Use of mathematical programming techniques. This approach allows the solution of the contact problem with or without explicit elimination of unilateral constraints. Methods such as Lagrange's multipliers or penalties allow the elimination of the unilateral constraints. Usually these methods are based on the use of special finite elements derived to simulate the impenetrability condition between two surfaces [8]. On the other hand, the unilateral constraints can be maintained in the formulation, retaining the original mathematical characteristics of the problem, and alternative linear complementary problems (LCP) can be obtained and solved by specific algorithms such as those developed by Lemke and Dantzig [1, 2, 8–24].

Since the seventies, these two general approaches have been employed for the numerical treatment of different types of contact problems. Following the first approach, Stadter and Weiss [3] used special gap elements and an iterative procedure which adjust the modulus of these elements to model the contact process between structures. An analytical model of the postbuckling behavior of a beam under unilateral contact constraints imposed by a rigid smooth surface was derived by Adan et al. [4], and Holmes et al. [5] studied the elastic buckling of an inextensible beam with hinged ends in the presence of rigid, frictionless sidewalls which constraint the lateral displacements. A nonlinear modal solution methodology capable of solving equilibrium problems of structural systems (beams and arches) with contact constraints was presented by Silveira and Gonçalves [6]. More recently, Li and Berger [7] presented a semi-analytical approach for 3D elastostatic normal contact problems with friction.

Numerical simulation of contact problems based on mathematical programming can also be easily found in literature. Fundamentals of the unilateral contact boundary value problem, including friction, together with finite element applications to the solution of the variational inequalities arising in static and dynamic structural contact problems can be found in Panagiotopoulos’ articles [9–12]. Stavroulakis et al. [12] give additionally the mechanical interpretation of Lemke’s algorithm [13, 14] for the contact problem. More specific applications of numerical methods for the unilateral contact problem between plates and an elastic foundation are presented in [15–17]. Contact constraints present in the stability of rods and large postbuckling behavior are presented, respectively, by Stein and Wriggers [18] and Simo et al. [19]. Algorithms designed to trace the nonlinear response using path-following techniques, such as arc-length procedures [20], are proposed by Wriggers et al. [21], Koo and Kwak [22], and Silveira and Gonçalves [23]. A good review and discussion of the complexity of physical effects that can occur at the contacting surfaces are presented by Barber and Ciavarella [25]. Finally, algorithms and discretization techniques for structural contact problems can be found on Wriggers’ book [24].

These two general approaches are employed in this work for the treatment of unilateral contact constraints. Since 1990, the authors have analyzed several contact problems involving a deformable structure and an elastic foundation [2, 6, 17, 23]. These works attempted to establish reliable numerical methodologies for the analysis of structures with unilateral boundary conditions. This article is an extension of these previous works where the comparison of alternative formulations for the analysis of structures resting unilaterally on an elastic foundation may be considered as a new contribution to the field.

The basic structural unilateral contact equations and inequations are stated in the next section. Two numerical methodologies are then proposed to solve the contact problem. In the first methodology, the contact problem is treated as a constrained optimization problem and the finite element method is used for the structure and for the Winkler-type elastic foundation. Two alternative linear complementary problems (LCPs) are derived and solved by Lemke’s algorithm [13, 14]. The second proposed approach for treating contact problems is a semi-analytical model that uses a Ritz-type methodology with moveable boundaries where the coordinates defining the limits of the contact regions are considered as additional variables. The unilateral contact problem is treated as an unconstrained optimum design problem. At the end, three particular structural elements under unilateral contact constraints (tensionless foundation) are used to validate the proposed formulations. The first method is more general and can be extended to a large number of problems, including three-dimensional problems. However, the second approach may lead to accurate low-dimensional models, as shown herein, with analytic coefficients which can be used in, for example, sensitivity analyses, optimization, feedback control problems, and parametric analyses.

2. Structural Unilateral Contact Problem: Basic Equations

Consider the structural contact problem shown in Figure 1. Note that the elastic foundation reacts only in compression. The structure is defined as a solid elastic continuum which occupies a domain V, limited by three different regular boundaries: S_c, S_u, and S_f. The displacements are prescribed in S_u and the forces in S_f. S_c is the region where boundary conditions are “ambiguous.” For the structure, the equilibrium equations, the cinematic relations, and the constitutive equations are given by where σ_ij are the Cauchy stress components, are the infinitesimal strain components, u_i are the displacements, and are the material parameters.

Figure 1: Structure under unilateral contact constraints imposed by an elastic foundation.

If the elastic foundation is described by Winkler’s model, then the following constitutive relation can be written: where u_b and r_b are the displacement and reaction, respectively, of the elastic foundation and C_b is the foundation modulus.

For the structural system studied here, the following boundary conditions must be satisfied: where u_i is the deflection of the structure orthogonal to the foundation and is the gap between the structure and the foundation in the potential contact region S_c. Inequality (2.3c) is the compatibility condition that represents the impenetrability between the bodies. Contact occurs when . If > 0, there is no contact.

As the elastic foundation reacts only to compression, the following inequality must also be satisfied on S_c: where compressive reactions are assumed to be positive. Finally, in the potential contact region, the following complementarity relationship between and r_b should be verified:

These restrictions define in a complete way the contact as being unilateral. Figure 2 shows the domain of validity of these three relations and the contact law.

Figure 2: Domain of validity of the contact constraints.

The solution of the unilateral contact problem can be obtained by solving (2.1a), using (2.1b) and (2.1c), together with boundary conditions (2.3a) and (2.3b), contact inequalities (2.3c) and (2.4), and the complementarity condition (2.5). The nonlinearity due to the unilateral constraints makes it difficult to solve the contact problem directly. For this reason, an equivalent minimization problem is formulated which is particularly suitable for numerical analysis. It can be shown that the optimization problem [2, 16, 17] whereis equivalent to the contact problem described above. Based on these equations, two alternative solution strategies are proposed in the next sections for the numerical analysis of structures resting on a tensionless foundation.

3. Constrained Contact Problem Formulation

According to Ascione and Grimaldi [15], restrictions (2.3c), (2.4), and (2.5) can be substituted by the variational inequality where belongs to the positive cone J, in which the admissible reactions r_b are the elements of and and Y are the vectorial spaces that contain the solutions of r_b and , respectively. See that the satisfaction of the variational inequality (3.1) for all belonging to J is equivalent to the contact constraints (2.3c), (2.4), and (2.5). Thus, the contact constraint may be eliminated from the analysis by writing

The first variation of , after eliminating u_b from the previous equation by way of relation (2.3c) and according to (3.1), is given by the following variational inequality:

Elimination of from (3.4), by the use of (2.2) and (2.3c), leads to a variational inequality in terms of u and r_b only, which corresponds to the first variation of the following functional: where . The variables u and r_b must be obtained in such a way that the first variation of the functional satisfies the inequality .

Using the finite element method, one can assume that for a generic structure and a foundation finite element the displacement and reaction fields within the element, u and r_b, are related to the nodal displacements u and nodal reactions r _b by where N and H _b are the matrices that contain the interpolation functions that describe, respectively, the behavior of the structure and elastic base.

From these definitions and adding the contributions of each finite element, one arrives at the discretized functional of the problem in the global form where F _e is the nodal load vector, K is the stiffness matrix of the structure, C is the joining matrix between the structure and the elastic foundation, and T is the flexibility matrix of the elastic foundation. These last two matrices are defined bywhere m_c is the number of elements of the contact region.

After the first variation of (3.7), one arrives at the following linear complementarity problem (LCP) in terms of the structure displacements and foundation reaction [15, 17]:

Equation (3.9), when considering the constraints (3.10), may be solved using mathematical programming methods, in particular, pivoting techniques developed for complementary problems [13, 14]. However, it is first necessary to reduce the previous relations to a standard LCP form. This can be obtained through the use of the following definitions: where and are the positive and negative parts of the vector u [26].

Using these new variables, it is possible to rewrite (3.9) and (3.10) in the following form:withEquation (3.12) and constraints (3.13) correspond to a standard LCP which is solved here using Lemke’s algorithm [14].

3.1. The Dual Formulation

If the stiffness matrix K in (3.9) is positive definite, it is possible to establish a relationship between u and as follows:

Substituting (3.15) in (3.7), one arrives at a variational expression that is function of the nodal values of the base reaction only, that is, in which P is a symmetric positive definite matrix, J is a vector, and v is a constant, and are defined as follows:

Equation (3.16), with the foundation reaction constraint, characterizes the following quadratic programming problem (QPP):

Thus, considering the Kuhn-Tucker conditions of this QPP, one can derive a standard LCP similar to the one described by (3.12) and (3.13), where now, with w defined as a Lagrange multiplier introduced to represent the impenetrability condition between the bodies. Once known , u can be obtained from (3.15).

4. Unconstrained Contact Problem Formulation

A different strategy of solution is now proposed. This strategy assumes that the contact constraints (2.3c), (2.4), and (2.5) on S_c as well as the elastic foundation displacements can be indirectly introduced into the analysis by explicitly considering the coordinates defining the limits of the contact regions (t_k) as additional variables of the problem (see Figure 1(c)). Thus, for the structural element in contact with a tensionless Winkler-type elastic foundation and subject to conservative loads, the total potential energy functional can be rewritten as where U is the strain energy which is a function of displacement vector u and of the vector which contains the coordinates defining the limits of the contact regions (t_k). Note that the length of each contact region is a function of the system parameters and is not known a priori. is the external load vector.

If the Ritz method is applied, the following displacement field, written in matrix form, can be used to approximate :where the matrix contains the functions that satisfy the boundary conditions on S_u and the vector A contains the unknown coefficients.

Thus, substituting (4.2) into (4.1), considering a small variation in the total potential energy, and expanding (4.1) in a Taylor series, one obtain: or using a more compact matrix notation, in which the vector W contains the unknown variables of the problem (A and ), g is the gradient vector (out-of-balance load vector), and M is the Hessian matrix (representative stiffness matrix) of the structural system, which can be organized as follows:where K, S, and C are the stiffness, contact, and joining matrices, respectively. For equilibrium, the change in (4.4) should be stationary irrespective of δ A and hence, the equilibrium equations (gradient vector) arewhich consist of a highly nonlinear set of algebraic equations, involving polynomial or transcendental functions of t_k.

Algorithm 1 shows the iterative solution strategy used to solve (4.6). This numerical solution is based on Newton’s method which uses the Hessian matrix (second-order derivatives) to obtain the direction of the solution. In this scheme, there are two distinct phases required for the successful completion of a single load step:(i)a predictor phase, where initialization of A and (W ⁰) are necessary;(ii)a corrector phase, where these approximations are corrected to satisfy the equilibrium equations, the convergence criterion adopted here is based on the norm of the gradient vector g.

Algorithm 1: Nonlinear solution strategy based on the second-order Newton method.

5. Numerical Examples

Figure 3 shows the structural elements under unilateral contact constraints used to test the constrained and unconstrained solution strategies proposed in this research. A Winkler tensionless elastic foundation of modulus K is assumed in all the examples. In the proposed unconstrained nonlinear solution strategy, a convergence factor of is adopted and consistent units are used in all examples.

Figure 3: Structural elements under unilateral contact constraints and corresponding deformation patterns.

Consider as a first problem a beam of length L, bending stiffness EI, and moments M applied at the ends, as illustrated in Figure 3(a). For such a load system, the expected behavior of the beam is shown in Figure 3(b), where one region of contact is expected. Note that t, which represents the length of the contact region between the bar and the elastic foundation and is not known a priori, constitutes the additional variable of the problem. The results of this contact problem are presented in Figure 4, where the deflection along the beam axis is given for different values of the foundation modulus k = KL⁴/EI. The unconstrained results (Ritz method) are compared with those obtained using a finite element model and mathematical programming techniques (primal and dual constrained solutions), and are obtained using the following linear combination of harmonic functions: to approximate the displacement field. Here, i is the number of half-waves, n is the total number of modes necessary to achieve convergence, and W_i are the modal amplitudes. Excellent agreement is observed between these three different solution strategies. Note also that the contact region (and the corresponding displacements) decreases steadily as k increases, while the displacements of the noncontact region increase. One of the main characteristics of tensionless foundation is the dependence of the contact area on the foundation stiffness.

Figure 4: Deformed shape of the beam (Figure 3(a)) for several values of the foundation stiffness .

Now, consider the beam shown in Figure 3(c). It is subjected to moments at the ends and a load P at the center and rests on a Winkler tensionless elastic foundation of modulus K. The corresponding deformation pattern is shown in Figure 3(d), where two contact regions and one central noncontact region are observed. Here, t_i and t_f constitute the additional variables of the problem. Figure 5 shows the influence of the foundation stiffness on the behavior of the beam, where the lateral deflection along the beam axis is given for different values of the k = KL⁴/EI. Once more, there is good agreement between the results obtained by the dual-constrained numerical formulation and the Ritz method, in which the same displacement field (5.1) is adopted. Under unilateral constraints, the beam displacement w increases in the central region while the contact regions between the bodies decrease.

Figure 5: Deformed shape of the beam (Figure 3(c)) for several values of the foundation stiffness .

The last example, illustrated in Figure 3(e), is a slender circular arch of radius R, length 2γ R, bending stiffness EI, and membrane stiffness EA in contact with a Winkler tensionless elastic foundation of modulus K. For this structural system, the expected deformation pattern is shown in Figure 3(f), where a central contact region defined by the angles ± t is expected. Results obtained using the unconstrained and constrained solution strategies are presented in Figures 6 and 7 for an arch with [27] R/h = 500, γ 10^°, EI = 1.4, EA = 420, and P = 0.1. For the unconstrained Ritz solution approach, the following approximations are used for the tangential and transversal displacements, respectively, of a point along the arch centroidal axis:where i and j are the number of half waves, and U_i and W_j are the modal amplitudes.

Figure 6: Deformed shape of the arch (Figure 3(e)) for the foundation stiffness parameter k = KR⁴/EI = 10⁶.

Figure 7: Elastic foundation reaction for P = 0.1 and k = KR⁴/EI =10⁶.

In Figure 6, the variation of the lateral displacement w, considering half the arc, is plotted as a function of the angle for the foundation parameter k = KR⁴/EI = 10⁶. The conventional foundation model (bilateral contact) was also considered. One can observe that under unilateral contact constraints the arch displacement w increases on the noncontact region. The contact angle between the bodies 2 7.46^°is the same for all formulations. Figure 7 shows the variation of elastic foundation reaction with the angle θ.

6. Conclusions

Two numerical approaches to solve the unilateral contact problem between a structure and a tensionless elastic foundation are proposed in this work. The first proposed formulation is based on the finite element method and mathematical programming techniques. Two alternative linear complementary problems (LCPs) are derived—primal and dual—and solved by Lemke’s algorithm. These formulations consider explicitly the inequality constraints that characterize the unilateral contact problem. The primal formulation can be used without restrictions in several structural problems while the dual formulation is restricted to structural systems in which the structure’s stiffness matrix is positive definite. However, the primal formulation leads to higher-order matrices and, consequently, to lower-computational efficiency and more processing time and sometimes present numerical instability.

The second formulation, named here unconstrained formulation, is based on the Ritz method and an iterative solution strategy. This methodology is particularly suited for the analysis of simple problems where the number and location of the contact regions are known a priori but not its length. This methodology can substitute in these cases large and time-consuming finite element packages. It may also be used as a benchmark for more general and complex formulations. However, this methodology leads to highly nonlinear equilibrium equations (due to the existence of unknown boundaries). So efficient iterative strategies, as proposed in Algorithm 1, should be used to solve such problems.

The examples involving beams and arches under contact constraints show that there is an excellent agreement among the results obtained from the different formulations, corroborating the effectiveness of all procedures in the analysis of unilateral contact problems. Finally, the results clarify the influence of the unilateral contact on the behavior of the structure.

Acknowledgments

The authors are grateful for the financial support from the Brazilian National Council for Scientific and Technological Development (CNPq/MCT). They also acknowledge the support from Brazilian Steel Company USIMINAS, Foundation Gorceix, FAPEMIG and FAPERJ-CNE. Special thanks are also due to Professor Michael Engelhardt and his department at the University of Texas for their hospitality during the manuscript’s preparation.