# Constrained and Unconstrained Optimization Formulations for Structural Elements in Unilateral Contact with an Elastic Foundation

Authors: Ricardo A. M. Silveira; Wellington L. A. Pereira; Paulo B. Gonçalves
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

In this work, two numerical methodologies are proposed for the solution of unilateral contact problems between a structural member (beam or arch) and an elastic foundation. In the first approach, the finite element method is used to discretize the structure and elastic foundation and the contact problem is formulated as a constrained optimization problem. Only the original variables of the problem are used, subjected to inequality constraints, and the relevant equations are written as a linear complementary problem (LCP). The second approach is based on the Ritz method, where the coordinates defining the limits of the contact regions are considered as additional variables of the problem. The contact problem here is treated as an unconstrained optimum design problem. These proposed methodologies are then tested and compared using results from specific problems involving structures under unilateral contact constraints.


---

## Body

## 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 Figure1. 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: Sc, Su, and Sf. The displacements
are prescribed in Su and the
forces in Sf. Sc is the region where boundary conditions
are “ambiguous.” For the structure, the equilibrium equations, the cinematic
relations, and the constitutive equations are given by(2.1a)σij,j=0,(2.1b)εij=(ui,j+uj,i)2,(2.1c)σij,j=Cijklεkl,
where σij are the Cauchy stress components,
εij are the infinitesimal
strain components, ui are the displacements, and Cijkl are the material parameters.Structure under unilateral contact constraints imposed by an elastic foundation.
(a)
Engineering problem(b)
Engineering model(c)
Deformation patternIf the elastic foundation is described by Winkler’s model, then the following
constitutive relation can be written:(2.2)rb=Cbub,
where ub and rb are the displacement and reaction, respectively, of the
elastic foundation and Cb is the foundation modulus.For the structural system studied here, the following boundary conditions must be
satisfied:(2.3a)ui=u¯ionSu,(2.3b)Fi=σijnjonSf,(2.3c)φ=ub−ui≥0onSc,
where ui 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 Sc. Inequality (2.3c) is
the compatibility condition that represents the impenetrability between the
bodies. Contact occurs when φ=0. If φ > 0, there is no contact.As the elastic foundation reacts only to compression, the
following inequality must also be satisfied onSc:
(2.4)rb≥0,
where compressive reactions are assumed to be positive. Finally,
in the potential contact region, the following complementarity relationship
between φ and rb should be verified:(2.5)∫ScrbφdSc=0.These restrictions define in a
complete way the contact as being unilateral. Figure2 shows the domain of
validity of these three relations and the contact law.Domain of validity of the contact constraints.
(a)
φ⩾0(b)
rb⩾0(c)
rbφ=0(d)
Contact lawThe 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]
(2.6)MinΠ(u,ub)(2.7)Subject  to:−φ≤0,onSc,
where(2.8)Π=12∫VCijklεklεijdV+12∫ScCbub2dSc−∫SfFiuidSf,is 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(3.1)∫ScτφdSc≥0,
where τ belongs to the positive cone J, in which the admissible reactions rb are the
elements of(3.2)J={rb∈Y′,∫ScrbwdSc≥0,∀w∈Y,w≥0}
and Y′ and Y
are the vectorial spaces that contain the solutions of rb 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(3.3)Π1=12∫VCijklεklεijdV+12∫ScCbub2dSc−∫ScrbϕdSc−∫SfFiuidSf.The first variation ofΠ1, after eliminating ub from the previous equation
by way of relation (2.3c) and according to (3.1), is given by the following
variational inequality:(3.4)δΠ1=∫VCijklεklδεijdV+∫ScCb(ϕ+u)δudSc+∫Sc[Cb(ϕ+u)−rb]δϕdSc−∫ScϕδrbdSc−∫SfFiδuidSf≤0.Elimination
ofφ from (3.4), by the use of (2.2) and (2.3c),
leads to a variational inequality in terms of u and rb only, which corresponds to the first variation
of the following functional:(3.5)Π1=12∫vCijklεklεijdv−12∫ScDbrb2dSc+∫ScrbudSc−∫SfFiuidSf,
where Db=Cb−1.
The variables u and rb must be obtained in such a way that the first variation of the functional Π2 satisfies the
inequality δΠ1(u,rb)≤0.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 rb, are related to the nodal
displacements u and nodal reactions rb
by
(3.6)u=Nu,rb=Hbrb,
where N and Hb
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(3.7)Π¯1=12uTKu−12rbTTrb+rbTCu−uTFe,
where Fe 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 by(3.8)C=∑mc∫ScHbTNdSc,T=∑mc∫ScHbTDbHbdSc,where mc 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]:(3.9)Ku+CTrb−Fe=0,(3.10)Cu−Trb≤0,rb≥0,(Cu−Trb)Trb=0.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:
(3.11)u=u++u−,z1=Trb−Au,z2=Ku+ATrb−R,z3=−z2,
where u+≥0 and u−≥0 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:(3.12)w=q+Mz,(3.13)w≥0,z≥0,wTz=0with(3.14)M=[K-KCT-KK-CT-CCT],q={-FeFe0},z={u+u-rb},w={z2z3z1}.Equation (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 matrixK in (3.9) is positive definite, it is
possible to establish a relationship between u and rb
as
follows:(3.15)u=K−1(Fe−CTrb).Substituting
(3.15) in (3.7), one arrives at a variational expression that is function of
the nodal values of the base reaction rb
only, that is,(3.16)Π¯2=−12rbTPrb+rbTJ−v,
in which P is a
symmetric positive definite matrix, J is
a vector, and v is a constant, and
are defined as follows:
(3.17)P=CK−1CT+T,J=CK−1Fe,v=(FeTK−1Fe)2.Equation (3.16), with the foundation reaction constraint,
characterizes the following quadratic programming problem (QPP):
(3.18)MaxΠ¯2(rb)(3.19)Subjectto:rb≥0onSc.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,(3.20)M=P,q=−J,z=rb
with w defined as
a Lagrange multiplier introduced to represent the impenetrability condition
between the bodies. Once known rb, u can be obtained from (3.15).


## 3.1. The Dual Formulation

If the stiffness matrixK in (3.9) is positive definite, it is
possible to establish a relationship between u and rb
as
follows:(3.15)u=K−1(Fe−CTrb).Substituting
(3.15) in (3.7), one arrives at a variational expression that is function of
the nodal values of the base reaction rb
only, that is,(3.16)Π¯2=−12rbTPrb+rbTJ−v,
in which P is a
symmetric positive definite matrix, J is
a vector, and v is a constant, and
are defined as follows:
(3.17)P=CK−1CT+T,J=CK−1Fe,v=(FeTK−1Fe)2.Equation (3.16), with the foundation reaction constraint,
characterizes the following quadratic programming problem (QPP):
(3.18)MaxΠ¯2(rb)(3.19)Subjectto:rb≥0onSc.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,(3.20)M=P,q=−J,z=rb
with w defined as
a Lagrange multiplier introduced to represent the impenetrability condition
between the bodies. Once known rb, 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 Sc  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 (tk) 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(4.1)Π¯3(u,Sc)=U(u,Sc)−FeTu,
where U
is the strain energy which is a function of displacement vector u and of the vector Sc which contains the coordinates defining the limits
of the contact regions (tk). Note that the length of each contact
region is a function of the system parameters and is not known a priori. Fe is the external load vector.If
the Ritz method is applied, the following displacement field, written in matrix
form, can be used to approximateΠ¯3:(4.2)u=ΨA,where the
matrix Ψ contains the functions that satisfy the
boundary conditions on Su 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:(4.3)δΠ¯3=∂Π¯3∂AδA+∂Π¯3∂ScδSc+12[δAT∂2Π¯3∂A2δA+2δAT∂2Π¯3∂A∂ScδSc+δScT∂2Π¯3∂Sc2δSc]+0(δA3,δSc3)
or using a more compact matrix notation,
(4.4)δΠ¯3=gTδW+12(δWTMδW)+0(δW3),
in which the vector W contains the unknown variables of the problem (A and Sc), 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:(4.5)M=[KCCTS]=[∂2Π¯3∂A2∂2Π¯3∂A∂Sc∂2Π¯3∂Sc∂A∂2Π¯3δSc2],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) are(4.6)∂Π¯3∂W=gT={∂Π¯3∂A∂Π¯3∂Sc}={00}which
consist of a highly nonlinear set of algebraic equations, involving polynomial or
transcendental functions of tk.Algorithm1 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)
apredictor phase, where
initialization of A and Sc (W0)
are necessary;(ii)
acorrector 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.
(1) Initialize:W0 (A and Sc)(2) Iterations (equilibrium problem): k=1,2,…,Imax(i) Compute:M(k-1)(ii) Compute:g(k-1)(iii) Check convergence:∥gk∥/∥Fe∥≤ξYes: Go to step (3)No: Compute the corrections: δWk=−[M(k−1)]−1g(k−1)Updated variables: Wk=W(k−1)+δWk, go to step (2)(3) Print the variables and stop


## 5. Numerical Examples

Figure3
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 ξ=10−3 is adopted and consistent
units are used in all examples.Structural elements under unilateral contact constraints and corresponding deformation patterns.
(a)(b)(c)(d)(e)(f)Consider as a first problem a beam of lengthL, 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 = KL4/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:
(5.1)w=∑inWisen(iπxL)
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 Wi 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 (Figure3(a)) for several values of the foundation stiffness k=KL4/EI.Now, consider the beam shown in Figure3(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, ti and tf 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 = KL4/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 (Figure3(c)) for several values of the foundation stiffness k=KL4/EI.The last
example, illustrated in Figure3(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:(5.2)u=∑imUisin((2i−1)πθ2γ),w=∑jnWjcos((2j−1)πθ2γ),where i and j
are the number of half waves, and Ui and Wj are the modal
amplitudes.Figure 6
Deformed shape of the arch (Figure3(e)) for the foundation stiffness parameter
k = KR4/EI = 106.Figure 7
Elastic foundation reactionrb for P = 0.1 and k = KR4/EI =106.In
Figure6, the variation of the lateral displacement w, considering half the
arc, is plotted as a function of the angle θ  for the foundation parameter k = KR4/EI = 106. 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 Algorithm1, 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.



---

*Source: 786520-2008-03-05.xml*