Interval elastoplastic analysis of structures

Computers and Structures 151 (2015) 1–10

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Interval elastoplastic analysis of structures C. Yang, S. Tangaramvong ⇑, W. Gao, F. Tin-Loi Centre for Infrastructure Engineering and Safety, School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia

a r t i c l e

i n f o

Article history: Received 10 June 2014 Accepted 3 December 2014

Keywords: Structural safety Interval uncertainty Elastoplastic analysis

a b s t r a c t The present work considers the effects of uncertainties in a classical elastoplastic analysis. The uncertainties are related to the applied forces and the plastic material capacities, both of which are taken to lie within deterministic but bounded intervals. The primary aim is to obtain both upper and lower bound responses (e.g. of the displacement history). Our interval elastoplastic problem is formulated as an optimization problem which we refer to as an interval mathematical program with equilibrium constraints or interval-MPEC. We develop an efficient and robust algorithm to process such interval-MPEC accurately with any readily available nonlinear programming solver. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Traditional elastoplastic analysis, with its ability to provide a rich spectrum of stress and displacement responses, undoubtedly helps engineers to assess the safety and integrity of structures. As evidenced by the literature, the classical analysis approach, which is based on an ideal premise that all input data are precisely known, is well developed. Much less studied, in spite of its importance, is an assessment of the effects of uncertainties inherent in real-world engineering applications, such as manufacturing defects, construction misalignments and so on [1]. The main concepts and methodologies developed to handle uncertainties in engineering mechanics problems can be found in [2–4]. A well-known technique, and one that is used in the present study, to take into account the influences of uncertain data is an interval approach (or so-called ‘‘convex’’ model) [5–7] that assumes the uncertainties are nonprobabilistic but are described by a set of independent parameters lying explicitly within known continuous interval ranges. Typically, one would then compute, for any given load level, the two bound solutions (namely most favorable and least favorable) to some response variable of interest, thus enabling a more reliable assessment of the structural safety of the structure under consideration. An interval analysis is founded on the hypothesis that the interval of the values of any uncertain parameter is known. In the framework provided by Information Theory [8]—a main tool for constructing probability models—the support of the probability ⇑ Corresponding author. Tel.: +61 2 9385 5146; fax: +61 2 9385 6139. E-mail addresses: (S. Tangaramvong).

[email protected],

http://dx.doi.org/10.1016/j.compstruc.2014.12.004 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

[email protected]

distribution of the random variable that models this uncertain parameter is also known. If there is no additional statistical knowledge concerning the uncertain parameter, then the use of the maximum entropy principle [9] with the constraints defined by such an available information yields a uniform distribution for this parameter. The maximum entropy principle is adopted to efficiently construct priori stochastic models of random variables [10] and the uncertainty quantification [11]. Consequently, the interval approach is a particular case of a probability approach. In general, additional information are available when constructing a probabilistic model of an uncertain parameter (a nominal value, some algebraic properties, etc.) and, therefore, without any experimental data, the probabilistic approach is considered to be richer than the interval approach. The interval approach, as adopted herein, is appealing to engineers since it merely requires knowledge of the upper and lower bounds of the uncertain data, and yet can provide fruitful preliminary sensitivity-type solutions to some response variable, prior to any comprehensive probabilistic and/or imprecise probabilistic analysis. In the light of the aforementioned comments, we propose, in this paper, a novel mathematical programming based approach to perform the interval elastoplastic analysis of structures. Without undue loss of generality, we assume that both the externally applied forces and the yield capacities of the elastic perfectlyplastic materials employed are interval data. Assuming small deformations and monotonically applied loads, we extend a standard path-independent (so-called ‘‘holonomic’’) elastoplastic analysis (see e.g. [12]) to include interval data. We note that a holonomic assumption (instead of nonholonomic or path-dependent assumption) enables a sufficiently accurate prediction of the

2

C. Yang et al. / Computers and Structures 151 (2015) 1–10

elastoplastic behavior of the class of structures under consideration [13]. The work presented herein is an extended version of the conference paper [14]. The formulation adopted in the present work is, however, simpler in that it is expressed in terms of total variables, rather than finite step quantities as in [14]. Both methods, moreover, provide identical interval elastoplastic bound responses. A related reference on the effect of interval uncertainty in elastoplasticity, albeit only on its effect on cross-section moment capacities at first yield and at full plasticity, is the work of Erdolena [15]. That work, however, contains some physical inconsistencies in assuming that the interval parameters for the same geometric and material data are independent for the two moment capacities of a cross-section. The main features of the proposed interval elastoplastic analysis are as follows: (a) The interval holonomic elastoplastic analysis is carried out under a stepwise load control. For each incremental load step, the algorithm captures the two extreme most favorable and least favorable limits to some response variable (e.g. the two displacement bounds at some specified location), thus enabling complete upper and lower response curves to be mapped out. (b) Our solution scheme eliminates the need to carry out any expensive enumerative-type search for the extreme response limits, as typically suggested in the interval literature [16,17]. Such a direct determination is enabled by forming and solving a pair of nonstandard optimization problems which we refer to as interval mathematical programs with equilibrium constraints or interval-MPECs. The equilibrium constraints represent the complementarity conditions typically used to describe the inelastic material properties. (c) In addition to the well-known numerical difficulties associated with solving standard MPECs [18] (e.g. disjunctivity, nonconvexity and nonsmoothness), interval-MPECs contain additionally interval data that need to be accommodated. In spite of these, we propose an efficient, robust and accurate algorithm to process our particular interval-MPEC. We first transform it into a standard (noninterval) MPEC by some additional bound constraints representing the interval data. The standard MPEC is then converted into a conventional NLP problem, by using a penalty regularization technique, and solved iteratively to enforce the original complementarity conditions. The organization of the paper is as follows. In the next Section 2, we review the classical (noninterval) holonomic elastoplastic analysis, formulated as an instance of a mathematical programming problem known as a mixed complementarity problem (MCP). Section 3 transparently extends, as an interval-MCP, this holonomic state problem to incorporate the influence of interval applied loads and yield capacities. In Section 4, we present a pair of interval-MPECs, extracted from the interval-MCP with the explicit aim of obtaining the upper bound in one case, and the lower bound in the other, for a given load level. Section 5 presents a penalty NLP-based regularization technique to solve the pair of standard MPECs. Three numerical examples are given in Section 6 to illustrate application of the proposed interval elastoplastic analysis approach. Finally, some concluding remarks are drawn in Section 7. Vectors and matrices are indicated in bold. A real vector x of size m is indicated by x 2 Rm and a real m  n matrix A by A 2 Rmn . For brevity, a vector of functions fðxÞ : Rm ! Rn is written simply as f 2 Rn .

2. Noninterval holonomic elastoplastic formulation 2.1. Generic FE model For simplicity, we use the common and important class of frame structures to detail our approach. As is well-known, a frame can be modeled as an aggregate of ‘‘line’’ finite elements. The generic beam FE model used is shown in Fig. 1. Rigid body motion is excluded by retaining only the independent generalized stresses Q i 2 R3 (namely, one axial force Q i1 and two end moments Q i2 and Q i3 ). Associated with the stresses are the generalized strains qi 2 R3 (i.e. an axial deformation qi1 and two end rotations qi2 and qi3 , respectively). It is useful to note that for this natural FE [12] the material behavior is described directly by the elemental-level response. A conventional plastic hinge concept is adopted, where the potential locations of plastic hinges are a priori assumed to be at the two ends of each element, while the response between these ends remains elastic. Such a model is appropriate for structures composed of ductile standard steel sections since plasticity is strongly localized at a limited number of critical zones. The external loads af are monotonically applied solely at nodes, where a and f denote a (positive scalar) load multiplier and a vector of given basic nodal forces, respectively. A distributed load can be approximated as a series of point loads located appropriately at nodal points. Under such a proportionally applied load regime, a holonomic (path-independent) behavior provides a sufficiently accurate prediction to the actual response [13] and is thus assumed in this work. The holonomic problem is described next. For clarity of exposition, a simple frame example (Fig. 2) is used to illustrate in detail key expressions and formulations throughout this paper. The two supports are fully restrained in all directions, and u defines a vertical displacement (m) associated with the applied force. For all members, steel sections with elastic perfectly-plastic material properties are employed, with Young’s elastic modulus of E ¼ 2  105 MPa, second moment of area of I ¼ 5  106 m4, cross-sectional area of A ¼ 2  103 m2, nominal flexural force yield capacity of Q 2u ¼ 20 kNm and nominal axial force yield capacity of Q 1u ¼ 500 kN. Nodal forces with a nominal magnitude of 10a kN each are horizontally and vertically applied as shown. 2.2. Holonomic state problem For a discretized frame of n elements, m generalized stresses (and strains), d degrees of freedom and y yield functions, the holonomic elastoplastic response of the structure can be extracted from the three basic static, kinematic and constitutive relations. These can be written in terms of standard notation and description as follows [12]:

CT Q ¼ af;

ð1Þ

q ¼ Cu;

ð2Þ

(a)

(b)

Fig. 1. Generic in-plane frame element i (a) stresses, (b) strains.

3

C. Yang et al. / Computers and Structures 151 (2015) 1–10

10α, u

2m

10α

2m

2m Fig. 2. Simple frame example.

q ¼ e þ p; Q ¼ Se;

ð3Þ ð4Þ

p ¼ Nz;

ð5Þ

w ¼ NT Q þ r P 0; z P 0; wT z ¼ 0:

ð6Þ

Fig. 3. PWL hexagonal yield locus; solid line denotes nominal limit and dashed lines interval bound limits.

The two relations given in (1) and (2), respectively, express linear equilibrium between the nodal forces af 2 Rd and the generalized stresses Q 2 Rm , and linear compatibility between the nodal

where h is an inclined angle of an undeformed frame member i, and l a member length. As an example of the plasticity conditions, consider the hexagonal PWL yield diagram in Fig. 3 for hinge a under the influence of

displacements u 2 Rd and the generalized strains q 2 Rm through a constant compatibility matrix C 2 Rmd . The total strains q are obtained in (3) from the additivity of elastic e 2 Rm and plastic p 2 Rm strains. The elastic material relations between the stresses Q and the elastic strains e are given in (4), where S 2 Rmm is the conventional (positive definite) elastic stiffness matrix. The assumed piecewise linear (PWL) holonomic elastoplastic constitutive laws are given in (5) and (6). An associative flow rule defines plastic strains p 2 Rm as functions of plastic multipliers z 2 Rm through a (constant) normality matrix N 2 Rmy that collects all unit normal directions to the yield hyperplanes. The PWL yield functions w 2 Ry in (6) are expressed in terms of the stresses Q and the yield limits r 2 Ry . Finally, (6) represents a complementarity condition (wT z ¼ 0) between the two signconstrained vectors w P 0 and z P 0. Such a condition describes the componentwise relationship wj zj ¼ 0; wj P 0 and zj P 0 for all j ¼ 1 to y, that implies either yielding (wj ¼ 0 and zj P 0) or elasticity (wj P 0 and zj ¼ 0). For a generic element i, the two key matrices Ci 2 R36 and Si 2 R33 are respectively

2

combined axial Q i1 and flexural Q i2 forces (as is typical of a standard steel I-section [19]), where rb ¼ 0:15; Q 1u and Q 2u denote the yield capacities associated with Q i1 and Q i2 , respectively. The complementarity conditions between the yield functions wa 2 R6 and the associated plastic multipliers za 2 R6 are

wa ¼ NaT Q a þ ra P 0; za P 0; waT za ¼ 0; where

waT ¼ ½ w1

EA=l

6 Si ¼ 4 0 0

0

0

0

6 Ni 2 R312 ¼ 4 1 0

w4

w5

w6 ;

~ ¼ Q 2u =Q 1u and s ¼ 1 þ r b tan c. For a pure bending model, where n the yield functions are obviously waT 2 R2 ¼ ½ w1 w4 . At the element i; wi 2 R12 simply assembles wa and wb for the two hinges a and b (with ra ¼ rb ), and hence can be written as

ð7Þ

wi ¼ NiT Q i þ ri P 0; zi P 0; wiT zi ¼ 0;

ð10Þ

where

3

7 4EI=l 2EI=l 5: 2EI=l 4EI=l

2

w3

z2 z3 z4 z5 z6 ; h i aT 2 Q 2 R ¼ Q i1 Q i2 ;   ~ tan c n ~ tan c 0 n ~ tan c n ~ tan c 0 n ; Na 2 R26 ¼ 1 1 1 1 1 1 raT 2 R6 ¼ ½ Q 2u sQ 2u sQ 2u Q 2u sQ 2u sQ 2u ;

and

2

w2

zaT ¼ ½ z1

3

cos h sin h 0  cos h  sin h 0 6 7 Ci ¼ 4  sin h=l cos h=l 1 sin h=l  cos h=l 0 5;  sin h=l cos h=l 0 sin h=l  cos h=l 1

ð9Þ

 wiT ¼ waT

ð8Þ

Q iT

~ tan c n ~ tan c n

0

~ tan c n ~ tan c 0 n

   12 wbT ; ziT 2 R ¼ zaT zbT ; h i 2 R3 ¼ Q i1 Q i2 Q i3 ; riT 2 R12 ¼ ½ raT

~ tan c n ~ tan c n

0

~ tan c n ~ tan c n

1

1

1

1

1

0

0

0

0

0

0

0

0

0

0

0

1

1

1

1

1

1

3 7 5:

rbT ;

4

C. Yang et al. / Computers and Structures 151 (2015) 1–10

Therefore, for the whole structural system composed of n generic elements, the governing vectors w; z; r and Q , and matrices N and S can be assembled as concatenated vectors, e.g. wT ¼ ½w1T ; . . . ; wnT , and block diagonal matrices, e.g. N ¼ diagðN1 ; . . . ; Nn Þ, respectively. The compatibility matrix C is assembled, as usual, from the nodal degrees of freedom. After some straightforward algebraic manipulations, we formulate the holonomic state problem in the unknown variables ðu; zÞ as follows:

CT SCu  CT SNz  af ¼ 0;

ð11Þ

w ¼ NT SCu þ NT SNz þ r P 0; z P 0; wT z ¼ 0:

The state problem as written in (11) is an instance of the class of mathematical programs, known as a mixed complementarity problem or MCP [20]. The complete holonomic elastoplastic response of the structure can be obtained by successfully processing a series of MCPs (11), each of which is associated with a specified load level a. A commonly used complementarity solver is PATH [21], that is available within the mathematical programming environment GAMS (an acronym for general algebraic modeling system) [22]. For the simple frame in Fig. 2, the vectors r 2 R16 ; f 2 R3 , and matrices C 2 R36 ; S 2 R66 ; N 2 R616 are

2

0

1

6 C¼4 0

0

1

200000

6 6 6 6 6 S¼6 6 6 6 4

0:707

0:250 0:250

3

7 0:500 0:500 0:707 0:250 0:250 5;

0 2

0

0

0

0

1

0

0

0

0

0

0

3

0

0

0

0

0

141420

0

0

0

0

0

1414:2

7 7 7 7 0 7 7; 7 0 7 7 707:107 5

0

0

0

0

707:107

1414:2

0

2000 1000

0

1000 2000

0

0

0

0

0

0

0

0

The extension of the standard holonomic elastoplastic formulation given by MCP (11) to include interval uncertainties is obvious. We assume that the externally applied forces and yield limits are L

L

L

U

0

0

0

0

0

0

0

0

U

L

U

U

and

U

r 6 ½r ; r  6 r , respectively. Thus, the interval forces ½f ; f  are uncertain but bounded within a convex domain (with componentL

L

U

U

wise relations implying f k 6 ½f k ; f k  6 f k for all k ¼ 1; . . . ; d, where U

L

f k and f k are the two upper and lower bound values associated with the interval data, respectively), and similarly for the interval yield limits ½rL ; rU . The uncertainties in the yield capacities are described by the combined stress diagram in Fig. 3, where the interval bounds (dotted lines) imply isotropically shrinking or expanding yield loci without altering the original shape. Of course, in an analysis, members with assumed identical cross-sectional materials can be grouped to reflect this. The interval elastoplastic analysis in variables ðu; zÞ, obtained by simply replacing f and r, respectively, with the intervals U

L

U

CT SCu  CT SNz  a½f ; f  ¼ 0; T

T

w ¼ N SCu þ N SNz þ ½rL ; rU  P 0; z P 0; wT z ¼ 0:

ð12Þ

The interval holonomic state problem in (12) is a nonstandard MCP, which for each specified load a contains the complete feasible set of response variables u and z under all possible

0

0

0

0

6 1 1 0 0 0 0 0 0 0 0 0 0 0 6 6 6 0 0 1 1 0 0 0 0 0 0 0 0 0 N¼6 6 0 0 0 0 0 0:0471 0:0471 0 0:0471 0:0471 0 0:0471 0:0471 6 6 40 0 0 0 1 1 1 1 1 1 0 0 0 0

L

represented by interval quantities such that f 6 ½f ; f  6 f

L

f ¼ ½ 10 10 0 ;

0

3. Interval holonomic elastoplastic formulation

½f ; f  and ½rL ; rU  then takes the form of the following MCP involving intervals:

0

T

2

The elastoplastic analysis of this structure was carried out using the PATH solver as a series of MCP (11) solves with an initial incremental load size of Da ¼ 5 that was progressively decreased by Da ¼ Da=10 when the limit load level was approached. The algorithm terminated when u ¼ 26:658 mm with a ¼ 36:355. The complete holonomic elastoplastic a  u response is displayed in Fig. 4.

0

1

1

1

0 0 0 0 0 1

0

0

3

7 7 7 7 0 0 7; 0:0471 0:0471 7 7 7 5 0 0 0

0

1

1

rT ¼ ½ 20 20 20 20 20 23:529 23:529 20 23:529 23:529 20 23:529 23:529 20 23:529 23:529 :

L

50 40

α

30 20 10

0

0.005

0.01

0.015

u (m) Fig. 4. Simple frame: noninterval holonomic elastoplastic a  u responses.

U

combinations of the interval data ½f ; f  and ½rL ; rU . Following the interval literature, one can try to process such problems using exhaustive exploration (e.g. Monte Carlo) simulations. Thus, the aim is to exhaustively explore the solution space in order to obtain approximate bounds to the response variable of interest (e.g. displacement uj at some specified point). From the physical reasoning that each interval parameter would be either at its least or its greatest value when extreme lower or upper bounds responses are sought, we could obtain the critical maximum and minimum response values by simply solving all MCP (12) instances associated with the 2k combinations of the k interval parameters solely at their extreme (upper and lower) bound limits. Unfortunately, the computational effort to successfully process practical-size problems is often intractable. An alternative, though still computationally demanding approach, is to use a 0–1

C. Yang et al. / Computers and Structures 151 (2015) 1–10

programming formulation as was done for interval limit analysis [6]. The method we propose in the next section overcomes these computational difficulties; and, for large scale problems, a Monte Carlo approach, as we have used in this paper, appears to be reasonable for verifying to some extent the solutions obtained. We should also note that owing to the nonconvex nature of the problem, it is possible that nonuniqueness of optimal parameter sets exists in that the same extreme bound could be caused by different combinations of interval parameters, even by not necessarily extreme value parameters. We describe in the next section a quick and reliable method to obtain these response bounds. L U We assume for our simple frame in Fig. 2, interval forces ½f ; f  L U of [0.9, 1.1] and interval yield limits ½r ; r  of ½0:95; 1:05, both with nominal values of 1. Hence, the associated interval vectors L

U

f 2 R3 ; f 2 R3 ; rL 2 R16 and rU 2 R16 are

f

LT

¼ ½ 9 11 0 ;

f

UT

¼ ½ 11 9 0 ;

5

In the following, we describe some of the numerical techniques, for which we have had considerable success with, aimed at processing the MPECs that arise from similar engineering mechanics problems (e.g. [23–25]).

5. Solution algorithms We focus on solving both MPECs in (14) for a given load level a. The entire response bounds can be obtained by solving these MPECs for various other load multipliers. There exists a family of well-known solution algorithms which hold the promise for solving standard MPECs. Each particular algorithm is based on some regularization technique, in which the MPEC is suitably reformulated into a standard NLP problem involving a (positive scalar) regularization parameter q. A series of reformulated NLP subproblems are iteratively processed, with successively increasing (or decreasing) parameter q, until the preset tolerance on the original complementarity condition (e.g. wT z 6 106 ) has been satisfied.

rLT ¼ ½ 19 19 19 19 19 22:353 22:353 19 22:353 22:353 19 22:353 22:353 19 22:353 22:353 ; rUT ¼ ½ 21 21 21 21 21 24:706 24:706 21 24:706 24:706 21 24:706 24:706 21 24:706 24:706 :

4. Obtaining bounds for the interval problem Our aim is not to exhaustively explore MCP (12), but to extract from it, for any specified load multiplier a, the two extreme maximum and minimum bound limits related to the desired response variable (e.g. uj ). In order to achieve this, we set up a pair of optimization problems, one maximizing the response variable uj and the other minimizing the variable uj subject to the constraints describing the interval holonomic elastoplastic problem. These two problems in variables ðu; zÞ can be written as

maximize uj ðor minimize uj Þ L

Several parameterized functions, categorized by the way the complementarity term is treated, exist, three of which are described below. 1. Relaxation: The original complementarity conditions ðwT z ¼ 0Þ are relaxed by wT z 6 q. The relaxed problem is iteratively processed with successively decreasing the value q (e.g. q ¼ q=10) until the tolerance on the original complementarity condition has been achieved [26]. 2. Smoothing: The complementarity conditions are smoothed using the appropriate functions uq ðwl ; zl Þ ¼ 0 for all l ¼ 1; . . . ; y. A well-known smoothing function is the following Fischer–Burmeister function [27]:

U

subject to CT SCu  CT SNz  a½f ; f  ¼ 0; w ¼ NT SCu þ NT SNz þ ½rL ; rU  P 0; z P 0; wT z ¼ 0:

uq ðwl ; zl Þ ¼

ð13Þ These optimization problems appear to be new and we refer to then as interval-MPECs. In addition to the underlying difficulties associated standard MPECs, the interval quantities make them more challenging to solve in their present form. We propose in the following a simple and efficient reformulation of the interval-MPECs (13) as standard MPECs by simply L

U

replacing all interval data ½f ; f  and ½rL ; rU  by the two respective additional force ~f and yield limit ~r bounded variables, namely f 6 ~f 6 f and rL 6 ~r 6 rU . These relations imply the componentL U wise relationships f 6 ~f 6 f for all k ¼ 1; . . . ; d and r L 6 ~r 6 rU L

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi w2l þ z2l þ 2q  ðwl þ zl Þ:

ð15Þ

The function uq has the property that uq ðwl ; zl Þ ¼ 0 if and only if wl P 0; zl P 0 and wl zl ¼ q. The smoothing algorithm then consists of solving a series of smoothed NLP sub-problems with iteratively decreasing parameter q (e.g. q ¼ q=10), until the complementarity tolerance is achieved [23]. 3. Penalty: The penalty algorithm transfers the penalized complementarity condition to the objective function [28]. The original complementarity constraint is progressively satisfied by iteratively increasing the value of q. The penalized NLP sub-problem is

U

k

k

k

l

l

l

for all l ¼ 1; . . . ; y, respectively. The resulting standard MPECs in variables ðu; z; ~f; ~rÞ can be written as follows:

maximize uj ðor minimize uj Þ subject to CT SCu  CT SNz  a~f ¼ 0; w ¼ NT SCu þ NT SNz þ ~r P 0; z P 0; wT z ¼ 0; L U f 6 ~f 6 f ; rL 6 ~r 6 rU :

ð14Þ

maximize uj  qwT z ðor minimize uj þ qwT zÞ subject to CT SCu  CT SNz  a~f ¼ 0; w ¼ NT SCu þ NT SNz þ ~r P 0; z P 0; L U f 6 ~f 6 f ; rL 6 ~r 6 rU :

ð16Þ

For this study, these three parameterized schemes have been implemented and our preliminary testing has shown that all can successfully process the MPEC (14) concerned. In the present work, we have chosen the penalty NLP-based technique.

6

C. Yang et al. / Computers and Structures 151 (2015) 1–10

The algorithmic steps adopted to map out the complete interval elastoplastic bound responses using the penalty regularization are outlined in the following: Step (0): Initialization  Set stopping criteria, e.g. maximum number of incremental load steps, minimum incremental load size (e.g. 104 ), etc.  Initialize: initial incremental load step Ma (positive scalar); a0 ¼ 0; u ¼ 0; and z ¼ 0.  Go to Step (a).

Table 1 Simple frame: percentage differences %Diff between interval and noninterval responses, where MCS denotes 200,000 Monte Carlo simulations.

a

Method

umin

%Diff

umax

%Diff

15

MPEC MCS

1.708 1.719

16.99 16.51

2.409 2.399

16.99 16.54

30

MPEC MCS

3.417 3.437

16.99 16.51

5.310 5.243

28.88 27.34

Step (a): Optimization  At loading step i, update ai ¼ ai1 þ Ma.  Process the penalty algorithm for obtaining the maximum and minimum bound solutions for some specified response variable.  If a pair of MPEC solves are successful, collect the optimal solutions and repeat Step (a).  Else, go to Step (b). Step (b): Termination  If any of the stopping criteria has been met, terminate the algorithm.  Else, decrease the loading step size (e.g. Ma ¼ Ma=10) and return to Step (a). It should be noted that the algorithm at Step (b) automatically decreases the incremental load size Ma when the elastoplastic structural responses approach the collapse limit. The proposed algorithm terminates when the incremental step size is too small (namely less than a minimum tolerance) or any of the preset termination criteria at initial Step (0) has been met. The interval elastoplastic bound analysis for the frame in Fig. 2 was performed using the penalty procedure with an initial finite step of Ma ¼ 5, which was automatically adjusted when the limit load was approached. The desired responses are the two critical displacements maximum umax and minimum umin . The complete interval elastoplastic a  umax (dashed line) and a  umin (thick line) behaviors generated from 7 and 31 MPEC solutions respectively are shown in Fig. 5, corresponding to those points for u 6 0:018 m. Each of the MPEC runs employed an initial regularization parameter of q ¼ 1, which was iteratively increased by q ¼ 10q until wT z 6 106 . This Fig. 5 also presents the noninterval elastoplastic a  u behavior (thin line).

50 u = 0.015 m

40 α = 30

α

30

(b)

(a)

Fig. 6. Simple frame: hinge dispositions at (a) umin ða¼15Þ ; umax ða¼15Þ ; umin umax ða¼30Þ ; umax ¼ 0:015 m, umin ¼ 0:015 m;  denotes plastic hinge.

ða¼30Þ ;

Clearly, as the number of plastic hinges developed with increasing a, the interval responses in Fig. 5 showed the significant influences of interval forces and plastic capacities. For instance, the interval responses deviated from the noninterval counterpart behavior up to 42% at u ¼ 0:015 m. The interval a  umax and a  umin responses gradually approached the two maximum load bound limits, namely a lower bound limit of alo ¼ 31:398 and an upper bound limit of aup ¼ 42:415, respectively. When umax ¼ umin ¼ 0:015 m, aumax ¼ 30:500 and aumin ¼ 41:690. In Fig. 5 we also show the interval a  umax and a  umin responses (open dots) found from a series of complete enumerative searches using extreme interval value parameters, and for the same a levels as for the MPEC approach. For each a; 24 MCP (12) solves are required (associated with 24 combinations of 2 interval forces and 2 interval member plastic limits at their extreme bound values). Exact agreement between the bound responses computed by the interval-MPEC solves and the enumerative MCP searches was obtained, hence validating the accuracy of the present MPEC approach. The possibility of nonunique optimal sets of interval parameters, that we mentioned previously in Section 3, was evident in this example. At optimality, 6 of 7 interval-MPEC umax results and 30 of 31 umin results reported all critical interval parameters at their extreme bounds. For those two MPEC results with parameters not entirely at their extreme limits, the enumerative approach (which assumes all extreme value parameters) gave exactly the same bound values, thus indicating nonuniqueness of parameter sets. We also noticed from the enumerative optimal solutions that there were a few cases when different extreme value parameter sets provided the same optimal bound, once again indicating nonuniqueness. Table 1 details the percentage differences (%Diff) between umax and u, as well as umin and u at the two load levels of a ¼ 15 and 30, and the corresponding hinge dispositions are depicted in Fig. 6. The accuracy of the obtained interval response limit solutions was also partially validated by 200,000 Monte Carlo

20 α−u α − umin

10

u

α − umax

0

0.005

0.01

0.015

u (m) Fig. 5. Simple frame: interval holonomic elastoplastic a  umax and a  umin responses,  denote the responses obtained by complete enumerative searches.

[email protected] = 127.26 ft Fig. 7. Example 1: three-span arch bridge.

15 ft 8.785 ft

α

α = 15

(b)

7

C. Yang et al. / Computers and Structures 151 (2015) 1–10

simulated runs which gave umax and umin results within the associated interval-MPEC solutions.

8α 7.5α 7α 6.5α 6α 5.5α 5α 4.5α

6. Illustrative examples

6.1. Descriptions 6.1.1. Example 1: three-span arch bridge The first example is the three-span arch bridge in Fig. 7 with a span length of 42.42 ft, where all base supports are fully fixed and u defines a vertical displacement (ft) at the center of the mid bay. This bridge was subjected to a uniformly distributed load with a nominal magnitude of a kip/ft. The material properties adopted were those of a standard I-steel section with the nominal values of E ¼ 4:176  106 kip/ft2; A = 2 ft2; I = 2 ft4; Q 2u ¼ 5184 kip-ft; and Q 1u ¼ 10; 368 kip. All potential plastic hinges obeyed the combined stress hexagonal PWL yield diagram in Fig. 3. The interval

5α u

[email protected] m = 15 m

4α 3α 2α

6m

6m

Fig. 8. Example 2: five-storey frame.

2.25 m

1m

6m

1m

α

4α 3.5α 3α 2.5α 2α 1.5α α 0.5α

16@4 = 64 m

Three numerical examples are provided to illustrate application of the proposed interval elastoplastic scheme. For all examples, elastic perfectly-plastic I-steel sections were used for both beams and columns. Beams were assumed to yield in pure bending and columns under the combined bending-axial force PWL yield diagram in Fig. 3. Positive and negative plastic properties were chosen to be identical. Each of the three examples considered a large number of interval force and member yield limit parameters, which in effect made it computationally intractable to perform complete enumerative MCP searches associated with the interval parameters at their extreme bound values (as done for the previous small-size example); for each a level, the number of combinations are 2103 ; 2129 and 2494 for Examples 1 to 3, respectively. Instead, we attempted to (partially) validate the accuracy of the computed interval elastoplastic responses by using 200,000 Monte Carlo deterministic simulations for some selected load multipliers a. Additionally, we carried out interval limit analyses [6] to provide some reference results on the most maximum collapse load limit aup and the most minimum collapse limit alo . We implemented the interval holonomic elasoplastic analysis algorithm as a MATLAB code that was interfaced with a GAMS mathematical modeling framework. A transparent communication between these two modeling environments was enabled by a GDXMRW interface program [29]. Each NLP sub-problem (16) was processed using the general purpose NLP solver GAMS/CONOPT [30]. The CPU times reported were for a 3.4 GHz Intel Core i7 computer with 16 GB RAM, running a Win7 operating system.

u

8@8 = 64 m Fig. 9. Example 3: twin-tower structure.

ranges associated with the nominal values of applied forces and plastic limits were ½0:9; 1:1 and ½0:95; 1:05, respectively. The discrete model consists of 63 nodes, 62 elements, 177 degrees of freedom and 992 yield functions. The uniformly distributed load was approximated as 41 nodal forces. The associated interval forces were assumed to be independent, while the yield limits were taken to be the same over each span. The interval elastoplastic a  umax and a  umin responses were mapped out with an initial incremental load step of Ma ¼ 2. 6.1.2. Example 2: five-storey frame The second example concerns the five-storey frame in Fig. 8 [31], where all base supports are fully restrained and u defines a top-right sway displacement (m). The structure was monotonically subjected to interval horizontal forces (kN) with the nominal magnitudes shown in Fig. 8 and interval vertical forces with the nominal values of 6a, with the corresponding interval load variation of ½0:9; 1:1 being applied throughout. Material properties employed were those for standard I-steel sections with E ¼ 2  105 MPa: for all columns, 400WC328 with the nominal values of Q 2u ¼ 1988 kNm and Q 1u ¼ 11; 704 kN; and for all beams and inclined members, 460UB82.1 with the nominal values of Q 2u ¼ 552 kNm and Q 1u ¼ 3150 kN. Pure-bending plasticity was assumed for all beams, and the combined stress yield locus in Fig. 3 was adopted for all columns and braces. The interval yield limit variation of ½0:95; 1:05 was adopted. The structural model consists of 56 nodes, 74 elements, 150 degrees of freedom and 648 yield functions, in which the beam over each bay was discretized into 2 elements. The interval elastoplastic responses were traced using an initial load step of 4a ¼ 5. 6.1.3. Example 3: twin-tower structure The last example involves the interval elastoplastic analysis of the practical-sized twin-tower steel structure shown in Fig. 9, where all ground supports are fully fixed and u denotes a top-storey sway displacement (m). The structure was subjected to interval vertical forces with a typical nominal magnitude of 5a (kN) and interval lateral loads (kN) shown in this Fig. 9. An interval load variation of ½0:9; 1:1 was adopted. The standard I-steel sections with E ¼ 2  105 MPa employed were: for all columns, 500WC440 with nominal values of Q 2u ¼ 2912 kNm, Q 1u ¼ 15; 680 kN; and for all beams and braces, 460UB82.1 with nominal values of Q 2u ¼ 552 kNm,

8

C. Yang et al. / Computers and Structures 151 (2015) 1–10

140

Table 2 Examples 1–3: percentage differences (%Diff) between interval and noninterval responses, where MCS denotes 200,000 Monte Carlo simulations.

u = 0.25 ft

120 100

a

Method

umin

%Diff

umax

%Diff

1

40

MPEC MCS

0.0157 ft 0.0179 ft

19.14 7.946

0.0232 ft 0.0210 ft

19.14 7.94

80

MPEC MCS

0.0366 ft 0.0466 ft

32.48 14.020

0.0756 ft 0.0626 ft

39.59 15.61

40

MPEC MCS

0.00676 m 0.00745 m

18.55 8.36

0.00984 m 0.00918 m

18.55 10.64

80

MPEC MCS

0.0137 m 0.0154 m

23.90 9.24

0.0242 m 0.0206 m

34.38 14.26

18

MPEC MCS

0.283 m 0.295 m

10.55 6.80

0.350 m 0.335 m

10.55 5.86

36

MPEC MCS

0.692 m 0.789 m

22.01 11.09

1.261 m 1.000 m

42.08 12.74

α = 80

α

80

Example

60 2

40

α = 40

α−u α−u

min

20

α−u

max

3

0

0.05

0.1

0.15

0.2

0.25

u (ft) Fig. 10. Example 1: interval holonomic elastoplastic a  umax and a  umin responses.

(a)

(b)

(c)

(d)

(e) Fig. 11. Example 1: hinge dispositions at (a) umin ða¼40Þ ; umax ða¼40Þ ; (b) umin umax ða¼80Þ ; (d) umin ¼ 0:25 ft; (e) umax ¼ 0:25 ft;  denotes plastic hinge.

ða¼80Þ ;

(c)

Q 1u ¼ 3150 kN. A pure-bending plasticity was assumed for all beams, while for all columns and braces the combined stress yield locus in Fig. 3 was adopted. The interval variation associated with the uncertain yield limits of all members was ½0:95; 1:05. The structural discretization consists of 282 elements, 204 nodes, 588 degrees of freedom and 2040 yield functions. The interval elastoplastic responses were computed using an initial incremental load step of 4a ¼ 2. 6.2. Results For all of the three examples presented, the interval elastoplastic a  umax and a  umin bound responses were successfully mapped out using the proposed analysis procedure. The CPU times reported were: in Example 1 some 21 s for the complete a  umin

responses, some 28 s for the a  umax responses; in Example 2 some 22 s for the a  umin responses, some 25 s for the a  umax responses; and in Example 3 some 56 s for the a  umin responses, some 90 s for the a  umax responses. The full spectra of the a  umax (dashed line) and a  umin (thick line) responses are displayed in Figs. 10, 12 and 14 for Examples 1 to 3, respectively. In these figures, the noninterval (nominal) a  u elastoplastic responses of the structures are also plotted. Hinge dispositions corresponding to the critical response solutions at some selected load levels a are depicted in Figs. 11, 13 and 15 for the three respective Examples 1 to 3. The percentage differences (%Diff) between the interval displacement responses (umin and umax ) and the noninterval results (u) at some selected load levels a are given for each of the three examples in Table 2. The accuracy of the computed interval displacement bound solutions was partially validated through comparisons with the results obtained by the computationally expensive 200,000 Monte Carlo simulations, denoted as ‘‘MCS’’ in the same Table 2. Direct interval limit analysis computations [6] provided the two critical (viz. lower alo and upper aup ) collapse load bound limits of the structures involved: in Example 1; alo ¼ 88:833 and aup ¼ 120:003; in Example 2; alo ¼ 95:592 and aup ¼ 129:133; and in Example 3; alo ¼ 39:336 and aup ¼ 53:138. Clearly, the obtained interval elastoplastic a  umax and a  umin responses all approached the bound load limits given by these reference interval limit analysis solutions, namely in Example 1 when umax ¼ umin ¼ 0:25 ft, aumax ¼ 88:460 and aumin ¼ 118:196; in Example 2 when umax ¼ umin ¼ 0:08 m, aumax ¼ 90:711 and aumin ¼ 121:650; and in Example 3 when umax ¼ umin ¼ 2:5 m, aumax ¼ 38:183 and aumin ¼ 51:098. These results further validate the accuracy of the interval elastoplastic results. 6.3. Comments Some pertinent comments in relation to the computation of the interval elastoplastic analysis approaches adopted for the three typical examples presented are worthwhile mentioning: (a) The interval elastoplastic analysis approach we propose is an efficient and reliable numerical method to successfully map out under load control the complete bound spectra of the interval holonomic elastoplastic (a  umax and a  umin ) structural responses. (b) The penalty NLP-based algorithm furnishes (at least local) optimal solutions to all of the underlying MPECs (14), and is recommended as a solution algorithm for such problems.

9

C. Yang et al. / Computers and Structures 151 (2015) 1–10

140

60

u = 0.08 m

u = 2.5 m

120

50

100

40 α = 80

α

α

80

α = 36

30

60 40

20

α = 40

α−u α − umin

20

α = 18

α−u α−u

10

min

α−u

α−u

max

max

0

0.02

0.04

0.06

0.08

0.1

0

0.5

u (m)

1

1.5

2

2.5

3

u (m)

Fig. 12. Example 2: interval holonomic elastoplastic a  umax and a  umin responses.

Fig. 14. Example 3: interval holonomic elastoplastic a  umax and a  umin responses.

(a) (a)

(b)

(d)

(c)

(b)

(c)

(d)

(e)

(e)

Fig. 13. Example 2: hinge dispositions at (a) umin ða¼40Þ ; umax ða¼40Þ ; (b) umin umax ða¼80Þ ; (d) umin ¼ 0:08 m; (e) umax ¼ 0:08 m;  denotes plastic hinge.

ða¼80Þ ;

(c)

(c) The automatic reduction of the incremental load step size Da enables the maximum load limit for each of the interval elastoplastic bound responses to be approached. In particular, the interval responses in Figs. 10, 12 and 14 attain the maximum load capacities of the structures involved, which coincide with the critical collapse load bound limits reported by the interval limit analysis approach [6]. This supports the accuracy of our interval elastoplastic algorithm. (d) The interval elastoplastic structural responses in Figs. 10, 12 and 14 illustrate the importance of incorporating the influences of uncertain applied loads and yield capacities for a reliable structural safety assessment. The effects of interval

Fig. 15. Example 3: hinge dispositions at (a) umin ða¼18Þ ; umax ða¼18Þ ; (b) umin umax ða¼36Þ ; (d) umin ¼ 2:5 m; (e) umax ¼ 2:5 m;  denotes plastic hinge.

ða¼36Þ ;

(c)

data become increasing important (viz. up to some 40% deviation from the noninterval behaviors even with a small interval range of applied forces ½0:9; 1:1 and plastic capacities ½0:95; 1:05) as the magnitude of the load multiplier a increases and a sufficient number of plastic hinges have developed.

10

C. Yang et al. / Computers and Structures 151 (2015) 1–10

(e) The accuracy of the computed interval holonomic elastoplastic responses was partially validated using 200,000 Monte Carlo simulated runs for some selected values of the load multiplier a. Moreover, these enumerative-type simulations fail by some significant margin to reach the critical response bound limits obtained by the direct MPEC (14) computations presently proposed. (f) Finally, the results are of practical value in that they provide engineers with the range of possible responses in the presence of interval uncertainties, and thus in highlighting possible design risks involved in using other than nominal ones. For instance, in the case of Example 3, for a load of a ¼ 36 (as indicated in Fig. 14), the MPEC solutions give a range of u lying within umin ¼ 0:692 m and umax ¼ 1:261 m, where the nominal value (i.e. obtained by simply replacing the interval data with their average values) is u ¼ 0:887 m. Assuming that the deflection limit is ulim 6 1 m, both nominal u and minimum umin responses indicate satisfaction of the current design under this serviceability criterion, while it is not for the maximum result (viz. umax > u). The designer will need to assess which displacement response to adopt. If the value of umax were chosen, the structure concerned requires further strengthening, and on the other extreme end if the value of umin were used, a more economical design with some smaller member sizes could be adopted. 7. Concluding remarks Motivated by the need to account for uncertainties in an elastoplastic analysis, we propose a novel interval approach to obtain the full spectra of the extreme responses under interval force and yield capacity data. For the monotonically applied load regime, the interval elastoplastic behaviors are traced in a holonomic fashion by directly computing a series of critical bound limits to some response variable associated with each of the predefined load levels a. The proposed method thus by-passes the need to carry out some typically intractable enumerative analysis procedure. The formulation is cast initially as an interval-MCP that contains all possible solutions for all combinations of interval data. Since we aim to capture only a pair of extreme bounds for each load level, we transform this interval-MCP as an interval-MPEC and finally as a standard MPEC which we process using a penalty regularization technique. The efficiency and robustness of the proposed interval analysis method have been illustrated through a number of practicallymotivated engineering structure examples (three of which are provided). The results have been partially validated by Monte Carlo simulations and interval limit analyses. Finally, these examples also highlight the importance of assessing the influence of uncertain applied forces and yield limits for practical application, especially at the higher load levels when a sufficient number of plastic hinges have developed. Acknowledgements This research was supported by the Australian Research Council – Australia. The authors would like to thank the reviewers for their constructive comments on an earlier version of the manuscript.

References [1] Möller B, Beer M. Engineering computation under uncertainty–capabilities of non-traditional models. Comput Struct 2008;86:1024–41. [2] Elishakoff I, Soize C. Nondeterministic mechanics. New York: Springer Wien; 2012. [3] Ghanem R, Spanos PD. Stochastic finite elements: a spectral approach. New York: Springer-Verlag; 1991. [4] Schuëller GI, Pradlwarter HJ. Computational stochastic structural analysis (COSSAN) – a software tool. Struct Saf 2006;28:68–82. [5] Ben-Haim Y, Elishakoff I. Convex models of uncertainty in applied mechanics. Amsterdam: Elsevier Science; 1990. [6] Tangaramvong S, Tin-Loi F, Wu D, Gao W. Mathematical programming approaches for obtaining sharp collapse load bounds in interval limit analysis. Comput Struct 2013;125:114–26. [7] Zhang H. Interval importance sampling method for finite element-based structural reliability assessment under parameter uncertainties. Struct Saf 2012;38:1–10. [8] Shannon CE. A mathematical theory of communication. Bell Syst Technol J 1948;27:379–423. 623–59. [9] Jaynes ET. Information theory and statistical mechanics. Phys Rev 1957;108:171–90. [10] Kapur JN, Kesavan HK. Entropy optimization principles with applications. San Diego: Academic Press; 1992. [11] Soize C. Stochastic models of uncertainties in computational mechanics. Reston: American Society of Civil Engineers (ASCE); 2012. [12] Maier G. A matrix structural theory of piecewise linear elastoplasticity with interacting yield planes. Meccanica 1970;5:54–66. [13] Tangaramvong S, Tin-Loi F. A complementarity approach for elastoplastic analysis of strain softening frames under combined bending and axial force. Eng Struct 2007;29:742–53. [14] Yang C, Tin-Loi F, Tangaramvong S, Gao W. A complementarity approach for elastoplastic analysis of frames with uncertainties. In: Steven GP, Li Q, Zhang Z. editor, Applied mechanics and materials. vol. 553; 2014: p. 594–99, doi:http:// dx.doi.org/10.4028/www.scientific.net/AMM.553.594. [15] Erdolen A. Uncertainty definition in structural systems with elasto-plastic materials under the bending moment effect by using interval analysis. Mech Based Des Struct Mach 2013;41:111–22. [16] Chinneck JW, Ramadan K. Linear programming with interval coefficients. J Oper Res Soc 2000;51:209–20. [17] Zhang H, Mullen RL, Muhanna RL. Interval Monte Carlo methods for structural reliability. Struct Saf 2010;32:183–90. [18] Luo ZQ, Pang JS, Ralph D. Mathematical programs with equilibrium constraints. Cambridge: Cambridge University Press; 1996. [19] Cohn MZ, Rafay T. Collapse load analysis of frames considering axial forces. ASCE J Eng Mech Div 1974;100:773–94. [20] Ferris MC, Munson TS. Complementarity problems in GAMS and the PATH solver. J Econ Dyn Control 2000;24:165–88. [21] Dirkse SP, Ferris MC. PATH solver: a nonmonotone stabilization scheme for mixed complementarity problems. Optim Methods Softw 1995;5:123–56. [22] Brooke A, Kendrick D, Meeraus A, Raman R. GAMS: a user’s guide. Washington, DC: GAMS Development Corporation; 1998. [23] Tangaramvong S, Tin-Loi F. Collapse load evaluation of structures with frictional contact supports under combined stresses. Comput Struct 2011;89:1050–8. [24] Tin-Loi F, Que NS. Parameter identification of quasibrittle materials as a mathematical program with equilibrium constraints. Comput Methods Appl Mech Eng 2001;190:5819–36. [25] Tangaramvong S, Tin-Loi F. Mathematical programming approaches for the safety assessment of semirigid elastoplastic frames. Int J Solids Struct 2011;48:1011–23. [26] Ferris MC, Tin-Loi F. Limit analysis of frictional block assemblies as a mathematical program with complementarity constraints. Int J Mech Sci 2001;43:209–24. [27] Kanzow C. Some noninterior continuation methods for linear complementarity problems. SIAM J Mat Anal Appl 1996;17:851–68. [28] Tangaramvong S, Tin-Loi F. Simultaneous ultimate load and deformation analysis of strain softening frames under combined stresses. Eng Struct 2008;30:664–74. [29] Ferris MC, Jain R, Dirkse S. GDXMRW: interfacing GAMS and MATLAB. In: Technical Report. Madison, Wisconsin: Computer Sciences Department, University of Wisconsin; 2010. [30] Drud AS. CONOPT-a large-scale GRG code. ORSA J Comput 1994;6:207–16. [31] Spiliopoulos KV, Patsios TN. An efficient mathematical programming method for the elastoplastic analysis of frames. Eng Struct 2010;32:1199–214.