Optimal prestress of cracked unilateral structures: finite element analysis of an optimal control problem for variational inequalities

cwlputermethcds in applied mechanicsand engineering ELSEVIER

Computer Methods

in Applied Mechanicsand Engineering 123 (1995) 231-246

Optimal prestress of cracked unilateral structures: finite element analysis of an optimal control problem for variational inequalities G.E. Stavroulakis

’

Institute of Applied Mechanics. Department of Engineering Sciences, Technical University of Crete, GR-73100 Chania, Greece

Received 22 April 1994; revised 5 September 1994

Abstract The structural analysis problem concerning static stabilization of cracked structures through prestress reinforcement is studied in this paper. A linearly elastic structure with given crack interfaces along which frictionless no-tension no-penetration mechanism hold (the so called unilateral contact mechanism) is assumed. The prestressing control, which in most civil engineering applications is realized by externally placed prestressed cables, is modelled by means of prestressing (control) forces. The static analysis (state) problem turns out to be a variational inequality. The goal of the stabilization is taken to be the healing of the cracks in the structure with the minimum reinforcement among an available and technologically feasible set and this is expressed by means of a suitably written combined cost function. The arising optimal control problem is a nonsmooth and nonconvex one, thus appropriate techniques recently developed in the area of nonsmooth analysis are used for its solution. Numerical examples illustrate the theory and parametric numerical experimentation with various algorithms permit us discuss certain aspects of the computer application. Appropriate condensation techniques for the reduction of the problem dimensionality and the integration of the studied algorithms in a general purpose finite element software are given.

1. Introduction Stabilization and restoration of cracked (e.g., masonry or reinforced concrete) structures is usually performed by means of appropriately placed prestressed cables. Externally placed reinforcements are normally preferred in applications due to the fact that they can be placed with minimum disturbance of the existing structure and it’s users and due to the fact that frictional losses of prestressing forces are lacking (see e.g., the case study reported in [ 1 ] ) . The goal of the restoration will be the minimization of all crack openings, which ideally means that all cracks close after the application of the prestressed cables and then the structure regains it’s structural integrity (at least partially). The cost of the restoration, which in our case is considered to be proportional to the cost of the prestressed cables should naturally be the minimum possible as well. The analysis of the aforementioned problem presents two difficulties: the cracked structure obeys to nonclassical interface laws due to the fact that only compressive forces are transmitted between the crack sides while separation (crack opening) is possible under certain loading, thus the structural analysis problem (state ’ Alexander von Humboldt Research Fellow, J&r- und Forschungsgebiet fiir Mechanik, and Lehrstuhl C fiir Mathematik, RWTH Aachen, D-52062 Aachen, Germany. 0045-7825/95/$09.50 c 1995 Elsevier Science .%A. All rights reserved SSDl0045-7825(94)00739-X

232

G.E. Stavroulakis/Comput. Methods Appl. Mech. Engrg. I23 (1995) 231-246

Fig. 1. Optimal prestressingof a cracked structure. problem) turns out to be a variational inequality problem. In turn the optimal prestressing action, leads to a nonsmooth and nonconvex optimal control problem concerning structures governed by variational inequalities. Finite element analysis techniques for the solution of this problem will be studied in the present paper based on recently developed techniques in nonsmooth analysis. The proposed methods are useful for the study of prestressing restoration of two and three dimensional structures where engineering intuition and classical trial and error techniques do not always lead to an efficient solution of the problem. A simplified discrete crack model will be used for the static analysis of the cracked structure throughout this paper: discrete cracks of known (e.g., measured in situ) position and initial crack opening will be introduced in the finite element model of the structure, which is assumed to be linearly elastic in the parts outside of the crack interfaces. No-tension and no-penetration conditions between the crack sides are modelled by means of appropriately defined monotone and possibly multivalued interface constitutive conditions (i.e. relations between the traction normal to the crack interface and the crack opening). For simplicity crack interfaces are assumed

frictionless in this study. The prestressing control action is simply modelled (in view of the technological considerations already mentioned) by means of appropriately placed control forces. The study of the arising structure with unilateral contact interfaces leads to certain variational inequality problems [2-51, which have extensively been studied in the framework of nonsmooth mechanics [6,7]. The discrete crack model which is used here and several generalizations of it is discussed in detail in [4-lo]. The variational inequality problem is given in the next section and serves as the nonlinear state model for the studied optimal control problem. For the cost function we assume here a sum of two quadratic terms: the first one corresponds to the cost of the reinforcing forces (which in turn are assumed to be directly analogous to the whole cost of the restoration) while the second one measures the crack opening along the crack interfaces of the structure (which should be equal to zero if an efficient restoration is obtained). Despite the simple cost function adopted here the optimal reinforcing problem is a nonclassical one due to the fact that the state system is a variational inequality. Thus certain results from the theory of optimal control of structures governed by variational inequalities must be used (see e.g., [ 1 l-181 ). The optimal control problem is formulated in section three in the sequel and necessary optimality conditions are derived in section five. Towards this end certain results from the sensitivity analysis theory for variational inequalities are summarized in section four, mainly due to the fact that this step have an interest of its own. Our results are based on theoretical investigations of [ 19-221 (see aLso [ 23,241). Optimal control problems for variational inequalities have been studied in [ 11,12,25-271 by the method of regulatization and in [ 28,291 by the method of conical derivatives. In this paper we prefer to formulate the optimal control problem for variational inequalities as nonsmooth optimal control problem (see also [ 14,15,30,3 I] ) because the direct (i.e. without regularization) treatment of the inequality constraints has a serious advantage in mechanical problems [ 4,8] : the information concerning the variable mechanical behaviour of the structure (e.g., contact and separation regions in a unilateral contact problem) are not destroyed by this method, as it is

the case with the method of regularization. As it will be shown in the sequel the highly

nonlinear

nature of the state variational inequality makes

233

GE. Stavroulakis/Comput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

the optimal control problem nondifferentiable and nonconvex. In this sense only local necessary optimality conditions can be written which in our case involve the generalized gradient operator in the sense of Clarke (see e.g., [ 181). Thus a bundle-type algorithm (see [ 321) seems to be the natural candidate for the numerical solution of the problem on hand, as it is already done in [ 13-161. A heuristic algorithm proposed in [ 331 will also be considered here. The description of the algorithms and the results of the numerical experiments are summarized in section six. For an efficient application of the method to large-scale structural analysis problems a boundary integral-like reformulation of the problem is presented in the appendix, based on certain results of [ 10,34-361. Appropriate static substructuring techniques (static condensation) are used there in order to reduce the size of the optimal control problem by excluding from further considerations the linearly elastic parts of the structure. It is worth noting that these matricial reformulations allow for the coupling of the optimal control algorithms, which are studied here, with general purpose finite element analysis software.

2. Simplified static modelling of cracked structures inequality state problem

with prestressing

control: a parametkl,

.&ational

Let us consider the cracked structure shown in Fig. 1 (for the two dimensional case, e.g., a cracked wall). Let us denote by 0 the region of the R3 space which occupies the structure in its undeformed state and let us assume that 0 is made up of the complementary, nonoverlapping, subregions (parts of a multi-body structure) L?,, w = 1,. . . ,p with smooth boundaries To, w = 1,. . . ,p. Parts of the boundaries r,, w = 1,. . . ,p correspond to the adjacent parts of the crack interfaces which are assumed to be given (measured). Interfaces will be denoted by Tin hereafter. The remaining parts of the boundaries r,, w = 1, . . . , p are composed of the complementary and nonoverlapping parts T,,, where the displacements of the structure are given (support boundary) and rF, where external tractions on the boundary are prescribed (boundary loading). The whole structure is referred to the fixed, right-handed, Cartesian coordinate system 0(x, ~2x3). Moreover, we assume that a finite element discretization has been done and such that the parts of the structure L$,,, w = 1, . . . , p are discretized by means of classical finite element techniques, while the nodes lying on adjacent parts of the interfaces are arranged in couples such as to permit us to model the interface behaviour by means of equivalent relative displacement vs. interface traction laws between each specific couple along the interface. This nodeto-node modelling of the interface is consistent with the small displacement and small deformation assumption and it has been proven to be adequate for this kind of problems [ 5,34,37]. Moreover, it is assumed that a fine mesh is used in the finite element discretization of the contact interfaces and/or appropriate finite elements have been used such as to capture with sufficient accuracy stress singularities that arise near the unknown, due to the unilateral contact mechanism, crack tips. Within the framework of the aforementioned discretization let u be the n-dimensional vector (displacement variables of all nodes of the discretization), p be the energy equivalent n-dimensional loading vector and K be the (n x n-dimensional) stiffness matrix of the (unsupported) structure. Here n is the number of the displacement d.o.f’s of the discretized structure. For the needs of this study we assume that the crack interfaces of the structure are represented in the discrete model by means of s couples of adjacent nodes, being placed on the adjacent sides of the crack interface ri,. Moreover, the relative displacements of the interface couples of nodes in the direction normal to the interface are written in matricial form (see e.g., [ 51)

[unl = Au,

(1)

where A is a (s x n) geometric transformation matrix and [u,] denotes the s-dimensional relative displacements. Moreover, let b denote the s-dimensional vector of initial interface unilateral contact effect at the interface is described by the inequalities Au < b. Following similar notation the classical support boundary boundary are written in the general form of equations ru = ug,

vector of interface openings [ 51. The

(2) conditions

of the structure

along the r,, part of the

(3)

234

GE. StavroulakisIComput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

where r is a (p X n) geometric transformation matrix, rqi is a p-dimensional vector of given boundary displacements and p is the number of nodes which model the support boundary r,. Note that the general form (3) includes all kinds of equality displacement boundary conditions which may arise in a structural model, as it is for instance restriction of selected displacement degrees of freedom in local coordinate systems. The prestressing control action is assumed here to be realized by induced-stress actuators and is modelled by means of equivalent loadings. To this end the contribution of the q prestress forces or groups of prestress forces (e.g., groups of cables) on the assumption that each group will be given the same prestress load, are added to the loading vector p after appropriate transfoirnation which is performed by an (n x q) transformation matrix C. Thus an enlarged loading vector p is applied on the structure, where

p=p+cz.

(4)

The induced-stress actuators model which is adopted in this study via (4) covers also the case of induced-strain actuators if the constitutive relation of the control system (the actuator or the prestressing cable, respectively, in our case) is linear. In this sense the proposed model is suitable for the study of certain kinds of intelligent structures as well which incorporate elements of piezoelectric or shape memory alloys’ controllers. The solution of the considered static structural analysis problem leads to the following potential energy (quadratic) minimization problem (see e.g., [ 4, p. 261) on the usual small displacements and small deformation assumptions min{lir(u) = !u’Ku - p’u - z’dulu E Uad},

(5)

where the space of admissible displacements uad is defined by u,, = {U E R”IAU < b,

rU

= Ug}.

(6)

Here K is the (n x n) stiffness matrix of the discretized, unsupported structure which is a symmetric, positive semidefinite matrix. An equivalent formulation for the unilateral structural analysis problem is given by the following variational inequality problem: Find u E u,d such that (F(Z,U),V-U)=(Ku-p-Cz,L’-U)>O

kE&,j,

(7)

where (., -) denotes the inner product in the R” space and F( z, u) = Ku - p - Cz. For further reference we will also write problem (7) in the form of a generalized equation (see [4,6,381), called also a convex differential inclusion 0 E F(z,u)

(8)

+Nu,(U),

where NU, (u) denotes the normal cone of the space of admissible displacements Crti at the point u which is expressed here by

NV,(u) = {r E R”lr = Pp+A’A,

PE RP, AE R: suchthatA’=O,

Vi$?Z(z,u)}.

(9)

In (9), Ai denotes the i-th element of vector A and Z( z, u) denotes the subset of interface nodes where geometric contact conditions hold, i.e. the adjacent nodes on the i-th interface couple come in contact Z(z,u)

={iE

{l,...,

s} such that (A’,u) = b’}.

Here A’ denotes the i-th row of matrix A. In order to clarify their mechanical interpretation, relations (8)-(9)

(10)

are written explicitly as fOllows

Ku-p-Cz+r’p+Atl=O,

(11)

-Ih+u()=o

(12)

-Au + b +

Nqcl\)3 0.

(13)

THUS,vector p corresponds to the discrete nodal reactions at the displacement boundary condition (3), vector A corresponds to the discrete interface tractions along interface rt, (cf. (2)) and relations (13) are more

GE. Stavroulakis/Comput.

Methods Appl. Mech. Engrg. 123 (1995) 231-246

235

frequently written in structural analysis literature in the following complementary system of inequalities (cf. [4,51)

no penetration relation:

y = -Au + b 3 0,

no tensile tractions: complementarity :

(14) (15)

I\ b 0, (J&Y) = 0.

(16)

Problems

(7)) (8) and ( 1l)-( 12)-( 13) are well known from the mathematical optimization literature optimality conditions for the quadratic minimization problem (5). Moreover, problem (1 I)-( 12)-( 13) is also known as the generalized linear complementarity problem. All the aforementioned problems which describe the mechanical behaviour of the structure under the action of the variable control vector z will be considered as parametric state models in the sequel for the analysis of the optimal control model. For further reference, we denote the solution of problem (7) (or equivalently of problems (8) and ( 11) (12)-( 13)) for a fixed control vector z by the generally multivalued relations u E S(z), (z&&h) E B(z).

(17) (18)

( 17) and ( 18) are in general set-valued relations with convex, compact image sets S(z) (or E(z) , respectively) for each fixed z, which may be the empty set if the aforementioned problems are unsolvable. On the assumption that the support boundary conditions (3) prevent every (infinitesimal) rigid body motion of the structure (or of a part of it), and on the assumption that the finite elements which have been used for the derivation of the stiffness matrix K are appropriately chosen such as to prevent zero energy deformation modes of the structure (i.e. K is strictly positive definite on the subspace of R” where (3) holds) a unique solution of the problem exists for every value of control variable z and relations ( 17), ( 18) are one-to-one relations (this result follows from the general existence and uniqueness results for the solution of convex variational inequalities that describe the unilateral contact effects involved in our model, see [ 4 ] among others, and section four of this paper where the sensitivity analysis problem is discussed). Let us consider for the moment that a certain solution of problem (7) is computed for a given value of control forces z. Then for the unilateral contact interface, described by the inequalities (2), the following subsets are defined: l the subset of binding inequality constraints Z( z, u) , i.e. where geometric contact occurs, and relation ( IO) hold, l the subset of Z( z, u), denoted by ,7( z, u), where geometric contact is accompanied with a strictly positive contact traction, i.e. the strongly active inequality constraints J(z,u) l

> O},

(19)

the subset of semiactive inequality constraints defined by 3(z,u)

l

= {i E Z(z,u)(A’

=Z(z,u)\&7(z,u),

(20)

and the subset of separating inequality constraints defined by z-,=(1,2

,...,

s}\Z.

(21)

Following [ 221 (see also [ 141) we introduce the cone K(z,u)

= {V E R”I(F(z,u),

u) =0} n 7&(u)

={uER”I(F(z,u),

Y)=O, Fv=O,

(Aj, v)
Vj’jEz(z,u)},

and the linearity space of K( z, u) L(z,u)

= {t E R”((F(z,u),

Y) =0, I? =0, (Aj,t) =0

Vj ~Z(z,u)}

=K(z,u) fl(-K(z,u)). (22)

236

G.E. Stavroulakis/Comput.

Methods Appl. Mech. Engrg. 123 (1995) 231-246

In the presence of infinitesimal rigid body displacements and rotations for the structure or parts of it the stiffness matrix K is symmetric and positive semidefinite and the space of infinitesimal rigid body displacements and rotations for the unsupported structure is written (cf. [ 4,391) !R = {r E RNJr’Kr = 0). Using the definition sets

(23)

of F( z, u) given in (7) and (22)) we get the following expressions

K(z,u)={r~!J?\fr=O,

(A+

fZ(z,u)={r~%]rr=O,

(Aj,r)=O

3. Optimal prestressing

GO

for the aforementioned

Vj’jEZ(z,u),

(p+Cz,r)=O),

(24)

VjJjEZ(z,u),

(p+Cz,r)=O}.

(25)

reinforcement

A cost function f( z, Y(u(z))) is assumed which additively measures the cost of the control prestressing forces z and the cost of the crack opening vector y (see ( 14) ) . For the most simple case of a quadratic, additively decomposable cost function the optimal control problem reads

min f( z, y(u(z))) = i z’Mz + $y(u)‘Ny(u)

(26)

under the state and control space restrictions Y(U) = -Au

+ 6,

(27)

u(z) E S(z) 9

(28)

z E &d.

(29)

Here weight matrices M( q x q) and N( s x s) are assumed to be symmetric and positive definite. The admissible control space &d in (29) is assumed to be a closed, convex subset of R*. Frequently the case of box constraints will be considered, i.e. Zad={ZE

R’lLi<

Zi
(30)

i=I,...,y’}

where for prestressing control variables lower bounds Li are usually set equal to zero while upper bounds Ui will be dictated from technological restrictions of the prestressing system or from local strenght considerations (e.g., cracking of a concrete or masonry structure near the anchorange points of the tendons). An engineering interpretation of weighting matrices M and N in (26) is as follows: the two constants measure the relative gain of control forces and crack destruction while the values assigned to them define the relative emphasis posed on minimizing control effort versus minimizing allowable structural destruction. In fact cost function f in (26) is a compromise between the need to optimize the control action, represented by the first term of f, and the goal of healing the cracked structure, the destruction of which is measured by the second term. More complicated cost functions may also be considered (cf. e.g., [ 27,401) In the sequel the necessary optimality conditions are derived for the optimal control problem (26)-(29) following the general theory of optimal control of systems governed by variational inequalities. It should be emphasized that the problem is in general a nonconvex and nondifferentiable one due to the fact that the feasible space of the (z, u) variables is nonconvex (recall that restriction (28) is equivalent to the whole set of equations, inequalities and complementarity conditions of problems (5)-( 8)). Moreover, relation (28) is generally multivalued and nondifferentiable. Therefore the classical adjoint structure technique of smooth optimal control theory can not be directly applied in our case.

4. Sensitivity analysis of the state variational

inequality

Let us consider how the solution of the state variational inequality problem (7)) denoted by S(z) or in general 8( Z) (see ( 17)) ( 18) ) changes when the control variables z undergo small variations. To this end certain results from the sensitivity and stability analysis theory of variational inequalities and convex optimization will be used here.

GE. Stavroulakis/Comput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

I,

:

231

I

Fig. 2. Active, semi-active and non-active unilateral constraints.

Let us assume that the solution ua of problem (7) for a given value of control vector ~0 is given and accordingly the subsets of binding Z( a, uo), strongly-active J( ~0, ~a), and semi-active constraints J( ~0, ~0) are known (cf. (lo), (19), (20), and Fig. 2). Following [38] (see also [ 14,15,21,22,41]) the following regularity condition is introduced: The Strong Regularity Condition (SRC) holds at a point (~0, ~0, m, ha) for the generalized equation (8) if the generalized equation

(31)

has a unique solution (vi(~), Ki (A) is used K,(A) = {A E R”IA’= 0

/+(I?),

Al(C)) for each L E R‘+ws . In (31) the following polyhedral cone

Vi E J’(za,ua),

A’ 3 0

Vi E J(za,ua)}.

(32)

On the assumption that the aforementioned SRC holds at (~0, ~0, w, Aa), the relations E(U) (and accordingly S(U), see (17) and (18)) are locally single-valued and Lipschitz [ 14,15,38,41]. REMARK 4.1. A necessary and sufficient condition for the validity of the SRC at (~0, ua, ,uo, &) is that the matrix Bii(za,uc) is nonsingular and the matrix -B~~(zo,uo)(BII(zo,uo))-~B~~(zo,uo) has positive principal minors, where [V&(y) ~ll(zo,~o)

=

p”!. (33)

(34) (35) Thus SRC holds true in our case if the engineering precautions discussed after relations ( 17) and ( 18) are taken because V,F( ~0. ~0) = K and the matrix --&I (Z&i)-i&2 has the meaning of a flexibility matrix for unit interface loads on the interface nodes along the part 3 of the interface (cf. [ 351 and the Appendix). REMARK 4.2. Solution stability of the state variational inequality which is guaranteed by the SRC is a desirable feature for the optimal control problem since it guarantees that a slight change of the control variables data will not change drastically the solution [41, p. 6591. With respect to the optimal control problem studied here, solution stability of the state model permits the use of Lipschitz nonsmooth optimization algorithms.

G.E. Stavroulakis/Comput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

238

Furthermore, the relation E( 2) is directionally differentiable at za and its directional derivative along the direction h E RS, S’ (zo, h) = (S’(zo, h), A’( zo, h)) = (~2, ~2, AZ) solves the following sensitivity generalized equation r

1

0

I

E

r’

A>

v2

A>

-r

E.c2

0

-47

h2,s

0

-A3

A,,.7

where sr is the cardinality of A2 and Az,r,,,g = 0. REMARK

-7

l-

V,F(zo,uo)

-V,F(zo,uo)h

4.3. Problem

min{$v:Kv2

of the index set 3, where Az,i, i E (3,

(36) is equivalent

- h’Cv2lrv2

to the following

= 0, A902 = 0, Ap2

(36)

+

3, L’in\T} denotes the appropriate subregions

quadratic minimization

problem

(cf. [ 421)

f 0},

(37)

which, compared with the state structural analysis problem (5), reveals the efficiency of this approach: in fact the same structural model with the same boundary conditions and modified interface conditions has to be solved for an auxiliary loading vector equal to Ch for the calculation of the sensitivity information. However, on the assumption that the strict complementarity holds at a point (~0, ua), i.e. Z( ~0, un) = z( ~0, ua) (equivalently J( ZO,ua) = 0)) the operator E(z) is differentiable at za and the derivative is given by (38) where

In general the construction of an element of the generalized Jacobian matrix aa (or 8s) in the sense of Clarke [ 181 is considered in [ 14-161. Under certain conditions which are given in the aforementioned publications but which are difficult to be checked in general the final result is obtained: an element of the generalized gradient aE(za) is constructed by considering an index set 21 with

(41) where matrix Dg, is defined as in equation

5. Optimality

(39) and is assumed to be nonsingular.

conditions

The most general first order optimality conditions can be written for the minimization of the composite function f( zn,y(q)) defined by (26), (27), (28) over the closed convex set Zad, restriction (29)) i.e. 0 E af(zoAzo))

+JVI(Z,).

cost (42)

G.E. Stavroulakis/Comput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

239

Under the assumption discussed in the previous section (see also [ 13,141) the generalized cost function f is a Lipschitz continuous function, thus (42) does have a meaning with a denoting the generalized gradient in the sense of Clarke, and can be exploited further [ 18, p. 721 to give 0 E V,~(ZO~Y(ZO)) + c~{~~Vyf(z~,~(z~))} 1

+Nz,(zo)

VPi E ~Y(~(zo))

= -ADS>

(43)

or explicitly 0 E Mz + co{p:N(b

- Auo)} +nlz,(zo)

i

Vpi E -AaS(

uo E S(zo).

(44)

Here coi{Xi} denotes the convex hull of vectors Xi. Note that the multivaluedness of the operator aS produces the vectors Xi in (43) and (44) through its elements pi. A combination of the necessary optimality conditions (42)-(44) with the sensitivity analysis results (see (36) and (37) ) permit us express the former in the following adjoint variable like form. Let (~0, ~0) be a local optimality solution of the optimal control problem. Moreover, let SRC holds at (zo, uo) and S(m) be differentiable at zo. Let p. be the unique solution of the adjoint linear equation 0 E K’po -A’N(b

- Auo) +&-L(zo,uo),

(45)

where L(zo, uo) = {Z E R”JrZ = 0, (Aj, Z) = 0, j E c7(zo,uo)}.

(46)

0 E Mzo -cp,

(47)

Then +Nz,(zo).

If, at the point (~0, ua) , S( ZO) is not differentiable at ~0 (i.e. the set of semi-active constraints 3 = Z\,? $0) then we consider all possible subsets of 3 i.e. the index-subsets (48)

{Z(zo) Ii E K(zo)} = P(z\J), and let the vectors pi, i E K( ZO) solve the adjoint generalized equations OE K’pi-A’N(b-Auo)

+~/L,(~o,uo),

(49)

where L(zo,

UO)

=

{Z

E

R”II’Z

=

0, (Aj, Z) = 0, j E J(zIJ,uo)(A~,

Z) = 0, j E Z(ZO, ~0)).

(50)

Then 0E

MZO

-c’

co{pi,i

i

E

K(z0))

+N~(zo>.

(51)

The proofs of relations (45)-( 51) are given in [ 13-161 and will not be repeated here. In the appendix a boundary variable reformulation of the aforementioned problems is presented which facilitates the coupling of the here presented method with various nonsmooth optimization packages, with general purpose finite element software featuring substructuring options and with unilateral contact solvers. 6. Numerical solution of the problem and examples Since the studied problem is in general a nondifferentiable and a nonconvex one, appropriate numerical techniques must be used for its numerical solution. Our goal will be the numerical solution of the necessary optimality conditions (44). Three different approaches have been tested on the problem studied here: ( 1) Classical subgradient optimization algorithms (also known as the Kiev method, see e.g., [ 43]), a direct outgrowth of the gradient optimization algorithm has been the first choice in our investigation. This algorithm reads: (a) Initialization: k = 0, Starting point zk; (b) k = k + 1, Iteration counter;

240

GE. Stavroulakis/Cotnput. MethodsAppl. Mech. Engrg. I23 (1995) 231-246 Cc> zk+’

= zk - mgk ( zk>, where ffk is an appropriately chosen step length (usually obtained by a line search method), and gk( zk) E df( zk) is an element of the subgradient of f; (d) Stopping rule: If 1Izk+l - zkl 1 < E then exit, else continue with step (b). Note that the convergence of the subgradient algorithm is not guaranteed due to the nondifferentiability of function f. Furthermore, a reliable stopping criterion is lacking, due to the fact that the subgradient at a minimum point need not be zero, instead it could be a nonzero set that contains the element zero. The speed of convergence is also dramatically slow. Despite its deficiencies the subgradient algorithm is sometimes used due to its easy implementation. Our limited numerical experience shows that the convergence of a subgradient type algorithm depends on the starting point in an unpredictable way and the speed of convergence is extremely slow. (2) Use of an appropriate nondifferentiable optimization algorithm, from which methods based on the bundle approximation of the generalized gradient af( ~0) in the (42) has also been considered. Under the minimal assumption that the value of the function f and one member of its’ generalized gradient can be computed, one of the currently available implementations of the bundle method [ 3244-471 can be coupled with the finite element unilateral contact solver with the ability to provide the required sensitivity information. First applications of this technique have been reported in [ 13-161. The bundle concept is in some sense an extension of the subgradient algorithm such as to cope with the nondifferentiability of the cost function: a set valued, polygonal approximation is used as a discrete substitute of the set valued subgradient. This approximation is constructed iteratively along the iteration steps of the algorithm. Based on this approximation a piecewise linear approximation of the function is used in each step. Two families of bundle algorithms have been developed until now, to the authors best knowledge: the pioneering bundle algorithm of C. Lemarechal [32,44,48] and the Bundle Trust algorithm of H. Schramm and J. Zowe [46,47]. The basic version of the Bundle Trust algorithm for convex, nondifferentiable functions and for nonconvex, nondifferentiable functions under box type constraints (cf. (30)) (codes BT and BTNCBC, respectively [ 46,471) have been used in the numerical experiments which are reported in the sequel. (3) A heuristic optimization algorithm proposed in [ 331 under the title generalized gradient decomposition is another possibility. This algorithm is restricted to box-type constraints for the control variables (see (30) >. The steps of the algorithm read (a) Box constraints: .ZL”= Li, Vi; (b) Starting point: Zs = Zp = i ( Li + Ui) . Set initial box: Ly = Li, @ = Ui; (c) k = k + 1, Iteration counter; (d) Update box constraints: If g f 0 then Z,ik+’= $, Uik+’ = Uf, if 2 2 0 then Uf+’ = z:, Ly’ = Lf;

(e)

New point: zik+’ = i(Li”

+ U!+‘);

( f) Stopping rule: If 1I$+’ - z/l ) < E stop, else continue with step (c) . Since the latter algorithm is based on heuristic concepts no convergence results can be obtained. Nevertheless, if this scheme converges at all this happens in much less number of iterations than in the previous exact algorithms. Thus in very complicated (e.g., highly nonconvex, with a large number of local minima) optimization problems the last heuristic scheme could be considered as an alternative with engineering value. Let us consider the plane stress problem shown in Fig. 3. For the plate fi (height = 2.70, length = 12.0, all quantities are assumed dimensionless in this section) discretized by means of quadrilateral isoparametric finite elements, classical boundary (support) conditions are assumed to hold along boundary r,. Along the boundary r’” unilateral contact conditions are assumed to hold. The initial gap between boundary r’, and the rigid support is assumed to be linearly varying with a maximum value equal to 0.24 below node 12. For the plate material Hooke’s law is assumed with elasticity modulus E = 1.00 * E3 and Poisson’s ratio Y = 0.3. Let the plate be loaded by the constant loading vector P and the variable control loading 2 = [ 21, . . . , Z6I, as it is depicted in Fig. 3. The aim of the numerical investigations will be to find the optimal control of the mechanical behaviour along the unilateral contact boundary Tin. In particular the minimization of the gaps between boundary r’” and the rigid support is sought by using minimum control effort, in the sense of cost function f(z) of (26). For simplicity weight matrices M and N in cost function f are assumed to be diagonal matrices. Moreover, all

241

G.E. Stavroulakis/Comput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

P = 3 x 10.

z,-Z,J

Fig. 3. Elastic structure with unilateral contact boundary and control forces.

-

(b)

1 .- . heuristic

1

Cc)

2

. -- - iterative

3

control

- -* - no_control

4

5

heuristic

-

iterative

-

no_control

interface

(

6

variables

Fig. 4. Effect of control action on the example of Fig. 3. Comparison between heuristic and iterative algorithms. (a) Contact stresses, (b) interface relative displacements, and (c) optimal control variables.

diagonal elements of M (or N, respectively) are assumed to be equal to m (or it, respectively). A comparison between the results obtained by the heuristic algorith and the iterative optimization algorithms is presented first. For weight values m = 0.005, R = 1000.0 the results of the optimal control problem computed with the heuristic algorithm of [ 331, the subgradient optimization algorithm and the bundle trust optimization algorithm are shown in Fig. 4. For reference the contact stress distribution and the initial crack opening (gap) for the case with zero control forces are also given in Fig. 4. A comparison between optimal control problems with free and constrained control variables (resourses) is presented in the sequel. For m = 0.0, n = 10000.0 the results are given in Fig. 5. Here case a corresponds to the unconstrained optimal control problem, i.e. Z,, = R6, case b corresponds to the case of box constrained control variables with upper bound equal to 20.0 (i.e. 0 Q zi < 20.0, i = 1,. . . ,6) and for the cases c and d the upper bound is taken to be equal to 10.0 and 5.0, respectively. The minimum cost calculated is equal to

GE. StavroulakisIComput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

242 -a

-b-c

-d

Fig. 5. Effect of box constraints in the control space. Case a: unbounded, case b: upper bound = 20.0, case c: upper bound = 10.0, case d:

upper bound = 5.0. (m = 0.0,

n =

lONO.0) (a) Contact stresses, (b) interface relative displacements, and (c) optimal control variables.

0.535.10m3 for case a (code BT), 0.6872. lo-* for case b (code BTNCBC), 0.3420 for case c and 1.7397 for case d. Solutions obtained with slightly different values of m do not differ essentially, although different values for the optimal cost are calculated. The last series of numerical experiments demonstrates the role of the weight matrices M and iV in the optimal control problem. For m = 0.0005results for case a: n = 1000.0, case b: PZ= 10000.0 and case c: n = 100000.0 are presented in Fig. 6. The choice of the values to be used for the weight matrices is case dependent and it is dictated from the requirements that both the solution of the resulting optimal control problem should be possible with bounded numerical errors and that the compromise cost function f should be representative of the modelled real-life application,

7. Discussion The feasibility of solving optimal prestress problems for unilateral contact structures in the framework of the finite element method has been demonstrated in this paper. A quadratic compromise cost function has been considered for the optimization problem. Recent advances in the study of the optimal control problems for variational inequalities and nondifferentiable optimization software have been used. Efficient implementation with appropriate substructuring techniques has been discussed as well (see also the Appendix). A small number of numerical experiments have been reported. It may be briefly concluded that: the problem of optimal prestress of unilateral structures (e.g., cracked structures) can be studied with the finite element method in connection with nondifferentiable (nonsmooth) optimization techniques, further investigations concerning the sensitivity of the solution with respect to the weight coefficients which are used in the cost function, as well as the dependence of the solution on the mesh discretization may be performed in order to exploit various practical aspects of the proposed method, and the appropriateness of the quadratic compromise cost function is case dependent. The use of other, more general cost functions and the extension to control problems under multiple loading cases remain open questions for future research.

G.E. StavrouhkidContput.

Methods Appl. Mech. Engrg. 123 (I 995) 231-246

243

-a-b-c

;

45

-

40

.’

f 35 0 so

,I ._._........ ,*__..-~~

..

,,,_

,,I’

y r 25’ _,: ,_: r c 20 ,I’ Oe,s .._____....-’ ,I’ 1 ,’___... _.~’-._ ._. , s lo-. ,’ , ....I. -. - : .I.....*<. . 5.’ ..;,...:I_- ~-.. . . _. . I.__<-‘-4 ol 1 4 6 2 3 5 control variables Cc) Fig. 6. Effect of relative weight in the cost function. m = 0.0005 aad for case a: n = 1000.0, for case b: n = 10000.0 and for case c: n = 100000.0 (a) Contact stresses, (b) interface relative displacements, and (c) optimal control variables.

Acknowledgement

The author would like to express his gratitude to Professor P.D. Panagiotopoulos, Aristotle University, Thessaloniki, Greece, and RWTH Aachen, Germany, for his continuous encouragement and support, and to Professor 3.V. Outrata, Czech Academy of Sciences, Praha, Czech Republic, and University of Jena, Germany, for helpful discussions. The permission of Professor J. Zowe, University of Jena, Germany, to use the bundle trust optimization software is gratefully acknowledged. Thanks are also due to the Alexander von Humboldt Foundation, Bad Godesberg, Bonn, Germany, for financial support.

Appendix A. Boundary

integral reformulation

for the state and sensitivity problems

Boundary integral type reformulations of the unilateral contact and the sensitivity analysis problems in the spirit of [ 10,34-361, are summarized here. These techniques drastically reduce the size of the aforementioned subproblems, which constitutes a significant gain since these subproblems have to be solved several times in the course of the numerical optimization algorithm. These reformulations are generalizations of the well-known static condensation (substructuring) techniques of linear finite element analysis. The static analysis of the state unilateral contact problem leads to the following Linear Complementurity Problem (LCP) (see ( 1 1)-( 13)). In this section classical equality boundary conditions ( 12) are dropped for the sake of notational simplicity. Their effect is nevertheless assumed to be incorporated in the derivation of the (effective) stiffness matrix K which in turn is assumed to be positive definite. Ku+A’A=p+Cz

(A.11

Au+y=b

(A.21

y 2 0, A 2 0, y’l = 0

(A.31

GE. Stavroulakis/Comput. Methods Appl. Mech. Engrg. I23 (1995) 231-246

244

A.I. State LCP with respect to intelface or boundary stresses Solution of (A.l) for u and substitution in (A.2) results in y = b - yN.p - FN.~ z + FN.N& y 2 0, A > 0, y’h = 0,

(A.4) (A-5)

where FN.N = AK-‘At is a boundary flexibility matrix (cf. [ 10,35]), FN.~ = AK-‘C is an influence matrix which gives the interface (or boundary) displacements for unity control forces, and yN.p = AK-‘p is the vector of the initial relative displacements normal to the unilateral interface or boundary (effects of external loading, thermal stresses, etc.). A.2. Solution differentiability: sensitivity of interface stresses Besides the general formulation of problem (36)) LCP (A.4), (AS) can also be treated with the method of [41, p. 6751. The directional derivative of A along direction q E Rq, A’( zo, q) solves the following mixed LCP

(--FN.~)TT+ (-FN.~ )JV+ z>[(-F~.~)jq+

(FN.N)~.zz~ + (FN.N)~.~z~ =& (FN.N)~.~z~ + (FN.N)~.~z~ =O, (FN.N)~.~zJ + (FN.N)~.JzJ] ~0,

zy b 0, Z& = 0.

(A.61 (A.7)

(A.81 (A-9)

Moreover, A(zc) is F-differentiable at zo if either the semi-active index set ,? is empty, or the following relation holds true (-FN.~)~

- (FN.N)~.~~FN.N)~!~~(-FN.~)~

~0.

(A.lO)

A.3. State LCP with respect to interface or boundary relative displacements Without restriction of the generality we assume here that the displacement vector u (and all dependent variables in the discretized formulation of the problem) is rearranged such that u = [ut 1~21t and that the geometric transformation matrix A has the form A = [ZIO] (in fact variables ut are the nodal displacements normal to a prescribed unilateral boundary). After the transformation of the variables U{ = b - u1 we get the following state LCP A = PN.p - KN.Nb + KN.~ z + A /> 0 9 u’l/V> 0 Atu’1 = 0 ’

(A.ll)

KN.Nu:,

(A.12)

where pN.p = p, - K12K;‘p2 Fe the initial stresses along the unilateral boundary due to external loading of the structure, KN.N = K11 - K12KG1K21 is the stiffness matrix of the substructure defined by ut (cf. [ 10,35]), and KN.* = CI - K12K;,‘C2 is the influence matrix which contains the interface tractions due to unitary control forces. For the general case of relative displacement variables the derivation follows [ 341. A.4. Solution differentiability: crack opening sensitivity As previously the structure of LCP (A. 11)) (A. 12) is used for the formulation of the following mixed LCP for the directional derivative of y along direction q E H, y’( ~0,q)

(KN.~)I,v + (KN.N)z,.z,YT~ + (KN.N)~,.TY~ = 0, (KN.~)T~+ (KN.N),J.Z,YZ, + (KN.N)J=.~Y~ ~0, y’j [(KN.~)JT+ (KN.N)J .GY~~ + (KN.N)~*~Y~] YJ 30,

yg=o.

(A.13) (A.14)

=O,

(A.15) (A.16)

245

GE. StavroulakisIComput. Methods Appl. Mech. Engrg. 123 (1995) 231-246

y( ~0) is F-differentiable at ~0 if either the semi-active set ,? is empty, or the following relation holds true (KN.~)~

-

(KN.N)~.Z,(KN.N)Z

.C'(--KN.~)Z,

~0.

(A.17)

References [ll V. Scholz, Anwendung extemen Spannglieder am Beispiel der Isarbtilcke Unterflihring, Bauingenieur 68 (1993) 151-158. [21 G. Duvaut and 1.L. Lions, Les inequations en mechanique et en physique (Dunod, Paris, 1972; English translation: Springer Verlag, Wien, 1974). 131 PD. Panagiotopoulos, Nonmonotone friction laws on interfaces. Applications to the mechanics of masonry structures, in: A. Giuffre and A. Grimaldi, eds., Studi Italiani Sulla Meccanica delle Murature (Roma, 1985) l-48. [41 PD. Panagiotopoulos, Inequality problems in mechanics and applications-Convex and nonconvex energy functionals (Birkhlluser Verlag, Boston, 1985; Russian translation: MIR, Moscow, 1989). [51 G.E. Stavroulakis, On the static behaviour of cracked masonry walls, Eur. J. Mech., A/Solids, 9(4) (1990) 341-358. [61 J.J. Moreau, PD. Panagiotopoulos and G. Strang, eds., Topics in nonsmooth mechanics (Birkh&user Verlag, Boston, 1988). [71 J.J. Moreau and PD. Panagiotopoulos, eds., Nonsmooth mechanics and applications, CISM Lecture Notes 302 (Springer Verlag, Wien, 1988). [81 PD. Panagiotopoulos, Hemivariational inequalities and their applications in mechanics and engineering (Springer Verlag, New York, 1993). t91 PD. Panagiotopoulos and C.C. Baniotopoulos, A hemivariational inequality and substationary approach to the interface problem, theory and prospects of applications, Engrg. Anal. 1 ( 1984) 20-3 1. [lOI G.E. Stavroulakis, Analysis of structures with interfaces. Formulation and study of variational-hemivariational inequality problems, Ph.D. Dissertation, Aristotle University, Thessaloniki, Greece, 1991. IllI PD. Panagiotopoulos, Subdifferentials and optimal control in unilateral elasticity, Mech. Res. Commun. 3 (1976) 91-96. [I21 PD. Panagiotopoulos, Optimal control and parameter identification of structures with convex and non-convex strain energy density. Applications to elastoplasticity and to contact problems, SM Arch. 8 (1983) 160-183. I131 J.V. Outrata, Necessary optimality conditions for Stackelbetg problems, J. Optim. Theor. Appl. 76(2) (1993) 205-320. [141 J.V. Outrata, On optimal control of systems governed by variational inequalities, DFG Report No. 351, University of Bayreuth, Germany, 1992; SIAM J. Optim. 2 (1994) in press. I151 J.V. Outrata and J. Zowe, A nondifferentiable approach to the solution of control problems with variational inequalities, DFG Report No. 368, University of Bayreuth, Germany, 1992. [I61 J.V. Outrata and J. Zowe, A numerical approach to optimization problems with variational inequality constraints, DFG Report No. 463, University of Bayreuth, Germany, 1993. [I71 J.-P Aubin and H. Frankowska, Set-valued analysis (Birkhguser Verlag, Boston, 1990). t181 F.H. Clarke, Optimization and nonsmooth analysis. (John Wiley & Sons, New York, 1983). [191 A.V. Fiacco, Sensitivity analysis for nonlinear programming using penalty methods, Math. Program. 10 (1976) 287-311. [201 M.P. Bendsoe, N. Olhoff and J. Sokolowski, Sensitivity analysis of problems of elasticity with unilateral constraints, J. Struct. Mech. 13 (1985) 201-222. [211 J. Kyparisis, Uniqueness and differentiability of solutions of parametric nonlinear complementatity problems, Math. Program. 36 (1986) 105-113. 1221 J. Kyparisis, Solution differentiability of variational inequalities, Math. Program. 48 (1990) 285-301. [231 1. Rentzeperis, Empfmdlichkeitsanalyse von Tragwerken mit einseitigen Bindungen, PhD Dissertation, RWTH Aachen, Germany, 1988. 1241 C.C. Baniotopoulos, A contribution to the sensitivity analysis of the sea-bed-structure interaction problem for underwater pipelines, Comput. Struct. 40(6) (1991) 1421-1427. [251 J.P. Yvon, Etude des quelques probltmes de contr6le pour des systemes distribues, These de Doctorat d’Etat, Universite Paris VI, France, 1973. 1261 V. Barbu, Necessary conditions for nonconvex distributed control problems governed by elliptic variational inequalities, J. Math. Anal. Appl. 80 (1981) 566-597. 1271 1. Haslinger and P NeittaanmIiki, Finite element approximation of optimal shape design: theory and applications (John Wiley & Sons, Chichester, 1988). t281 E MignSt, Contr6le dans les int!quations variationnelles elliptiques, J. Func. Anal. 22 (1976) 130-185. ~291 F. Mignh and J.P Puel, Optimal control in some variational inequalities, SIAM J. Control Optim. 22( 3) ( 1984) 466-476. 1301 J. Haslinger and T. RoubiEek, Optimal control of variational inequalities, Appl. Math. Optim. 14 (1986) 187-201. [311 M.M M%kell, Nonsmooth optimization. Theory and algorithms with applications to optimal control, Report 47, University of Jyvkkyla, Department of Mathematics, JyvLkyla, Finland, 1990. 1321 C. Lemarechal, J.J. Strodiot and A. Bihain, On a bundle algorithm for nonsmooth optimization, in: O.L. Mangasarian, R.R. Meyer and S.M. Robinson, eds., Nonlinear Programming, Vol. 4 (Academic Press, New York, 1981) 245-282. [ 33 I T.L. Friesz, R.L. Tobin, H.-J. Cho and N.J. Mehta, Sensitivity analysis based heuristic algorithms for mathematical programmes with variational inequality constraints, Math. Program. 48 (1990) 265-284. 1341 K.C. Khubring, Quadratic programs in frictionless contact problems, hit. J. Engrg. Sci. 24(7) (1986) 1207-1217.

246

G.E. Stavroulakis/Comput. Methods Appl. Mech. Engrg. I23 (1995) 231-246

L3.51 E. Mitsopoulou and IN. Doudoumis, A contribution to the analysis of unilateral contact problems with friction, SM Arch. 12 ( 1987)

165-186. ]36] H. Antes and PD. Panagiotopoulos, An integral equation approach to the static and dynamic contact problems. Equality and inequality Methods (Birkhluser Verlag, Boston, 1992). 1371 M.Ap. Tzaferopoulos, On an efficient new numerical method for the frictional contact problem of structures with convex energy density, Comput. Struct. 48( 1) (1993) 87-106. [38] S.M. Robinson, Strongly regular generalized equations, Math. Op. Res. 5 (1980) 43-62. [39] GE. Stavroulakis, PD. Panagiotopoulos and A.M. Al-Fahed, On the rigid body displacements and rotations in unilateral contact problems and applications, Comput. Struct. 40(3) ( 1991) 599-614. [40] K. Miettinen, M.M. MakeIs, An iterative method for nonsmooth multiobjective optimization with an application to optimaI control, Optim. Methods Software 2 ( 1993) 31-44. [41] R.W. Cottle, J.S. Pang and R.E. Stone, The linear complementarity problem (Academic Press, New York, 1992). [42] M. Jittorutrum, Solution point differentiability without strict complementarity in mathematical programming, Math. Program. 21 (1984) 127-138. [ 431 N.Z. Shor, Minimization methods for nondifferentiable functions (Springer Verlag, Berlin, 1985). [44] C. LemarechaI, M. Bancora Imbert, Le Module MlFCl (INRIA, Le Chesnay, 1985). [45] J.-J. Strodiot and V.H. Nguyen, On the numerical treatment of the equation 0 E af(x), in: J.J. Moreau and PD. Panagiotopoulos, eds., Topics in nonsmooth mechanics (Birkhauser Verlag, Basel, 1988) 267-294, [ 461 H. Schramm, Eine Kombination von Bundle- und Trust-Region-Verfahren zur L&ung nichtdifferenzierbarer Optimierungsprobleme, Bayreuther Mathematische Schriften, Heft 30, Bayreuth, 1989. [47] H. Schramm und J. Zowe, A version of the bundle idea for minimizing a nonsmooth function: conceptual idea, convergence analysis, numerical results, SIAM J. Optim. 2 (1992) 121-152. [ 481 K.C. Kiwiel, Methods of descent for nondifferentiable optimization (Springer Verlag, Berlin, 1985).