Interval field model and interval finite element analysis

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. xxx (xxxx) xxx www.elsevier.com/locate/cma

Interval field model and interval finite element analysis B.Y. Ni, C. Jiang ∗ State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, College of Mechanical and Vehicle Engineering, Hunan University, Changsha City, 410082, PR China Received 28 November 2018; received in revised form 23 October 2019; accepted 24 October 2019 Available online xxxx

Abstract Spatially uncertain parameters are traditionally represented by random field models. However, the large amount of information required for construction of precise probability distribution is often difficult to obtain for many practical engineering problems. In this paper, an interval field model is proposed to represent spatial uncertainties with insufficient information, in which the variation of the parameter at any location is quantified by an interval with upper and lower bounds. The spatial dependency is measured by a covariance function or a correlation coefficient function, defined for the interval variables at different locations. An interval Karhunen–Loève expansion is formulated for the proposed interval field model, in which the continuous spatial uncertainty is expressed through a series of deterministic functions with uncorrelated interval coefficients. Furthermore, by incorporating the interval field model into the finite element method, interval finite element analysis with spatially uncertain parameters is carried out. Perturbation-based interval finite element methods are developed to evaluate the upper and lower bounds of structural responses such as displacement and stress. A Monte Carlo simulation method is also presented to provide a reference solution for the structural analysis with interval fields. Finally, two numerical examples are investigated to demonstrate the effectiveness of the interval field model and the interval finite element methods. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Interval field model; Spatial uncertainty; Interval Karhunen–Loève expansion; Interval finite element method

1. Introduction The modeling of uncertain parameters with inherent spatial variability is commonly encountered in engineering. These uncertainties include material properties of the heterogeneous media such as concrete or porous rock, external loads applied on structures such as wind or snow loads, geographical parameters over the scale of meters such as soil permeability. This type of uncertainties are generally present in spatially variable properties, and are traditionally quantified by a random field model [1–4]. In random field models, the property at each spatial location is treated as a random variable, and the correlation between these random variables is taken into account. Based on random field models, the stochastic finite element method (SFEM) [5–9] was developed for structural analysis in the presence of spatial uncertainties. Moreover, random fields have been widely applied in various areas including civil engineering, geological engineering, materials engineering, etc. However, a large amount of information is generally required to determine ∗ Corresponding author.

E-mail address: [email protected] (C. Jiang). https://doi.org/10.1016/j.cma.2019.112713 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

2

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

the actual probability density function (PDF) for construction of a random field model, which is impractical or very costly to obtain in many engineering problems. For this reason, a series of non-probability or imprecise probability approaches such as evidence theory [10–12], fuzzy set theory [13,14], interval analysis [15–18], convex models [19–23], etc., are proposed and developed in recent decades to handle uncertainties with insufficient data. Among these approaches, the interval method plays a fundamental role for the following two reasons. First, despite nonprobability and imprecise probability methods quantify uncertainty by different mathematical descriptions, they are all closely related to an interval model. For example, in evidence theory, the uncertainty of a parameter can be recognized as a combination of interval focal elements with basic probability assignments (BPAs) [24,25]. And in fuzzy set theory, a fuzzy number can also be viewed as intervals with different membership levels [26]. Second, most of the scientific problems in uncertainty analysis that encountered in the interval method are also faced in other non-probability or imprecise probability methods. For example, in evidence-based reliability analysis, interval analysis of response function over each joint focal element is generally necessary for computation of the belief and plausibility functions of the structural reliability [27]. And in convex model theory, uncertainty propagation of the convex model variables eventually also lead to propagation analysis of an interval model [28]. Therefore, uncertainty analysis of structures using interval methods can be regarded as a basic requirement for further formulation of the other non-probability and imprecise probability methods [29]. In recent decades, the application of the interval method for structural uncertainty analysis has been studied extensively. Among them, the interval finite element method (IFEM) has been a very important development. It can be regarded as an uncertainty propagation method, which aims to calculate the response bounds of a structure with the input parametric bounds. In many cases, IFEM can also provide efficient computational tools for reliability analysis and optimization under uncertainty. In the 1990s, the problem of determining response bounds of structures with material and load uncertainties was studied using the finite element method (FEM) and the interval model [30,31]. It has been pointed out that the key problem in IFEM is to solve the non-deterministic equilibrium equations, where the interval stiffness matrix and the interval load vector are involved. As discussed in [15,16,32–35], the solution of this type of interval equations formulates an NP-hard (non-deterministic polynomial-time hard) problem. After that, different efforts have been made on IFEM to derive an accurate envelop of the solution set for the interval equations. Some of them are associated with fuzzy FEM [36–39], since the problem of solving the fuzzy equations can be transformed to the solution of a series of linear equations with interval variables. Muhanna and Mullen [40] proposed the element-by-element (EBE) formulation technique to avoid the overestimation problem in interval analysis. By modifying the displacement constraints used in the EBE formulation, Rama Rao et al. [41] proposed a new formulation for IFEM, which allows the simultaneous computation of primary and derived variables at the same level of accuracy. Moens and Hanss [42], Moens and Vandepitte [43] discussed the development of non-probabilistic finite element analysis and its applications in applied mechanics, through which the probability-, fuzzy sets- and interval-based finite element methods were investigated and compared. To account for the effect of fuzzy input uncertainties on structural responses, Rao and Sawyer [39], Rao and Chen [44] developed different versions of IFEMs for problems with small uncertainties. Regarding to problems where relations between responses and uncertain parameters are monotonic, Pownuk [45] presented an efficient algorithm for solution of fuzzy partial differential equations based on sensitivity analysis and FEM. Neumaier and Pownuk [46] proposed an iterative enclosure method to provide rigorous bounds on the solution of the linear systems whose coefficients have large uncertainties, which was demonstrated to be a fast algorithm for uncertain linear systems and has been selected as an effective approach for verification of some other newly proposed methods [47,48]. Qiu and Elishakoff [49], Qiu et al. [50,51] and Chen and Yang [52] applied the interval perturbation methods to estimate the upper and lower bounds of natural frequencies and static displacements of structures with bounded uncertain parameters. The above research mainly concentrates on interval finite element analysis of structures or systems with scalar uncertain parameters, which are quantified by interval variables. However, as mentioned before, many uncertain parameters in practical engineering problems possess an inherent spatial variability. Thus for problems with limited sample data, the development of non-probabilistic quantification models for spatial uncertainty and corresponding structural uncertainty analysis is very important. Moens et al. [53] proposed an interval field model which utilizes only upper and lower bounds for representation of the variability of a spatial parameter. The emergence of the interval field model, on the one hand provides an approach for uncertainty quantification of spatial parameters with probability statistics unavailable, and on the other hand, it also improves the theory of interval analysis and greatly extends its applicability. In the interval field model [53,54], an uncertain-but-bounded spatial parameter is Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

3

represented as a superposition of a series of base vectors using interval factors, through which the spatial dependency of the parameter can be considered. Afterwards, Muscolino et al. [55] extended this concept and proposed another interval field model based on the improved interval analysis via an extra unitary interval. Subsequently, this interval field model was applied to static analysis of beams with uncertain Young’s modulus under deterministic static loads [56–58]. Wu and Gao [59] carried out a static plane stress analysis of continuous structures involving interval fields, by which the upper and lower bounds of structural responses can be determined. The interval field model provides an effective quantification for spatial uncertainty when lacking sufficient experimental samples. It has been recognized that a proper formulation of the interval field model plays an extremely important role in relevant structural uncertainty analysis and reliability evaluation. Unfortunately, although some basic concepts of the interval field model have been established, an effective quantification method for describing and constructing the interval field conveniently and efficiently is still in absence. From the current research, it can be found that an expansion with interval variables is preferred in the representation of the interval field, where a series of base vectors or functions require to be determined. However, no consensus has yet been reached on the formulation of these base vectors or functions. In the interval field model proposed in [53], both the explicit interval field and the implicit interval field are discussed. In the explicit interval field, the basis vectors are reference patterns defined based on expert knowledge of the modeled system. And in the implicit interval field, the expansion form of an interval field seems more like a response surface, such as a Kriging response surface model. In the model given in [55], a Karhunen–Lo`eve (KL)-like decomposition is formulated as an expansion of the interval field. However, the upper and lower bounds of the interval field have to be passively determined by the KL-like expansion, which may cause conflicts to the predetermined bounds of the spatial uncertain parameter. Therefore, to develop an effective interval field model and establish a systematic approach for spatial uncertainty analysis have become a pivotal task, through which the practicability of the interval field model can be promoted to a certain extent. In this paper, a new interval field model is proposed to quantify the spatially uncertain parameters. An interval K–L expansion method is also developed for this interval field model, in which the continuous spatial uncertainty can be effectively represented by multiple interval variables, making it much more convenient to use in the subsequent structural uncertainty analysis. To perform the response analysis of structures with spatial uncertainties, the proposed interval field model is then introduced in the FEM, and the perturbation-based interval finite element methods (P-IFEMs) are developed for evaluation of the structural response bounds. The remainder of the paper is organized as follows: In Section 2, some basic concepts including the interval variable, correlation between intervals, and the interval process model which can also be regarded as the one-dimensional interval field model, are introduced. In Section 3, the introduction of the newly proposed interval field model is given. In Section 4, the interval K–L expansion and the simulation method for the interval field is presented. In Section 5, the interval finite element analysis is formulated in which the interval field model is employed to quantify the spatial uncertainties, and the first order and the second order P-IFEMs are developed for response bounds analysis of structures. A Monte Carlo simulation (MCS) method is also presented to give a reference solution for the P-IFEMs. Section 6 is devoted to numerical examples, and conclusions are remarked in Section 7. 2. Interval variables and interval processes 2.1. Interval variables As shown in Fig. 1, an uncertain parameter X is represented as an interval variable X I if its variation range can be determined, and it can be denoted as [15]: X ∈ X I = [X L , X U ]

(1)

where the superscripts I, L, and U represent interval, lower bound, and upper bound, respectively. The midpoint X m and the radius X r of the interval variable X I are defined as: XU + X L XU − X L , Xr = (2) Xm = 2 2 The variance D(X I ) of the interval variable X I is defined as: D(X I ) = (X r )2

(3)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

4

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 1. An interval variable.

Fig. 2. Correlation of two interval variables.

In the multidimensional case, the correlations among the variables are required since they can have considerable effect on the uncertainty analysis results. Currently, convex models [23,28,60–65] are used as the primary quantification measures for correlation of interval variables. In the ellipsoid convex model [23,66], the correlation of the two interval variables X 1 ∈ X 1I = [X 1L , X 1U ] and X 2 ∈ X 2I = [X 2L , X 2U ] is quantified through the ellipse enclosing all their bivariate samples, as shown in Fig. 2. The covariance Cov(X 1I , X 2I ) of the interval variables X 1I and X 2I is defined as: Cov(X 1I , X 2I ) = sin(θ ) cos(θ )(r12 − r22 )

(4)

where θ is the relevant angle, which is composed by the semi-major axis and the horizontal parameter axis, r1 and r2 are lengths of the semi-axes of the ellipse. The correlation coefficient ρ(X 1I , X 2I ), or ρ X I X I in brief, is defined 1 2 as: ρ(X 1I , X 2I ) =

Cov(X 1I , X 2I ) X 1r X 2r

(5)

I I The magnitude of the correlation coefficient ρ(X ⏐ ⏐ 1 , X 2 ) reflects ⏐the degree ⏐ of linear correlation of the interval I I I I ⏐ ⏐ variables X 1 and X 2 , and it satisfies ρ(X 1 , X 2 ) ≤ 1. A larger ⏐ρ(X 1I , X 2I )⏐ means a stronger linear correlation, ρ(X 1I , X 2I ) > 0 indicates positive correlation of the interval variables such as the cases in Fig. 3(a), and ρ(X 1I , X 2I ) < 0 indicates negative correlation of the interval variables such as the cases in Fig. 3(b). When the ellipse degenerates ⏐ ⏐ I I ⏐ ⏐ into a line, ρ(X 1 , X 2 ) reaches a maximum 1 which means a complete linear correlation. With the covariance defined by Eq. (4), the ellipse-shaped uncertainty domain Ω of the interval variables X 1 ∈ X 1I and X 2 ∈ X 2I can be analytically expressed as: ⎧ ⎫ ⏐ ]−1 [ ]T [ ] ⎬ ⎨[ ] ⏐⏐[ X 1 ⏐ X 1 − X 1m Cov(X 1I , X 1I ) Cov(X 1I , X 2I ) X 1 − X 1m ≤ 1 (6) Ω= m m X2 − X2 ⎭ ⎩ X 2 ⏐⏐ X 2 − X 2 Cov(X I , X I ) Cov(X I , X I ) 2

1

2

2

As shown in [23], Cov(X 1I , X 2I ) = Cov(X 2I , X 1I )

(7)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

5

Fig. 3. Different correlation cases of interval variables in the ellipsoid model.

and Cov(X iI , X iI ) = D(X iI ),

i = 1, 2.

(8)

where D(X iI ) = (X ir )2 is the variance of the interval variable X iI . For a multivariate problem with n interval variables X i ∈ X iI = [X iL , X iU ], i = 1, 2, . . . , n, the joint uncertainty domain ΩX of X = [X 1 , X 2 , . . . , X n ]T is described as an n-dimensional ellipsoid: { ⏐( } ) )T ( ⏐ (9) ΩX = X ⏐ X − Xm C−1 X − Xm ≤ 1 where Xm = [X 1m , X 2m , . . . , X nm ]T denotes the vector of midpoints, and the n-by-n matrix C is the covariance matrix: ⎡ ⎤ Cov(X 1I , X 1I ) Cov(X 1I , X 2I ) · · · Cov(X 1I , X nI ) ⎢ ⎥ ⎢Cov(X I , X I ) Cov(X I , X I ) · · · Cov(X I , X I )⎥ ⎢ n ⎥ 2 1 2 2 2 C=⎢ ⎥ .. .. .. .. ⎢ ⎥ . . . . ⎣ ⎦

(10)

Cov(X nI , X 1I ) Cov(X nI , X 2I ) · · · Cov(X nI , X nI )

2.2. Interval processes Recently, the authors extended the above interval variable case to the time-variant problems, thereby developing the interval process model [66,67] to quantify the time-variant or dynamic uncertain parameters such as road Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

6

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 4. Correlation for the interval process model.

excitations applied to vehicles. In the interval process model, an interval rather than a precise probability distribution is employed to describe the parametric uncertainty at each time point. The variation range of the whole timevariant uncertain parameter is described by a pair of upper and lower bounds. Thus the interval process can be also interpreted as a ‘time-variant interval variable’. By use of the upper and lower bounds instead of the accurate probability distributions, the interval process model provides a useful supplement to the traditional stochastic process for time-variant uncertainty quantification, and is expected to be applied to many complex structures with limited data information. A brief introduction to the interval process model is given in the following. Definition 1. A time-variant uncertain parameter {X (t), t ∈ T } is an interval process if for arbitrary time tk ∈ T, k = 1, 2, . . ., the possible values of X (tk ) can be represented by an interval X I (tk ) = [X L (tk ), X U (tk )], where T is a parameter set for the time t. } { An interval process can be denoted as X (t) ∈ X I (t), t ∈ T , or X I (t) in brief. Definition 2. For an interval process X I (t) with the upper bound function X U (t) and the lower bound function X L (t), the midpoint function X m (t) of X I (t) is defined as: X U (t) + X L (t) 2 The radius function X r (t) and the variance function D X I (t) of X I (t) are defined as: X m (t) =

X r (t) =

X U (t) − X L (t) 2

D X I (t) = (X r (t))2

(11)

(12) (13)

Definition 3. For an interval process X I (t), the auto-covariance function of interval variables X I (ti ) and X I (t j ) at any two time points ti and t j is defined as: Cov X I X I (ti , t j ) = Cov(X I (ti ), X I (t j )) = sin(θ) cos(θ )(r12 − r22 )

(14)

in which the relevant angle θ , the lengths r1 and r2 of the major axis and the minor axis can be obtained based on the samples x (s) (t), s = 1, 2, . . . , Ns of X I (t), as shown in Fig. 4. Herein Ns denotes the number of the samples. Definition 4. For an interval process X I (t), the auto-correlation coefficient function of interval variables X I (ti ) and X I (t j ) at any two time points ti and t j is defined as: ρ X I X I (ti , t j ) = √

Cov X I X I (ti , t j ) √ D X I (ti ) D X I (t j )

(15)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

7

Fig. 5. A stationary interval process.

As in the case of stochastic processes, interval processes can also be classified into two categories, stationary interval processes and non-stationary interval processes. Those processes whose uncertainty characteristics do not change with time are called as stationary processes. Stationary processes play a primary role in practical structural dynamics, since they are widely encountered in engineering. The magnitude of seismic waves during strong earthquakes, the turbulence process of ships, the fluctuation of voltage in lighting grids, and various noises and interferences etc., generally are all required to be treated as stationary processes. In the following, the definition of stationary interval process is provided. Definition 5. If the midpoint function X m (t) and the radius function X r (t) of an interval process X I⏐(t) are both ⏐ ⏐ constants and furthermore its auto-covariance function is only related to the time interval τ = ti − t j ⏐ rather than the time points ti and t j , namely, X m (t) = m X I

(16)

X r (t) = r X I

(17)

Cov X I X I (ti , t j ) = Cov X I X I (0, τ ) = Cov X I X I (τ )

(18)

where m X I and r X I are both constants, Cov X I X I (τ ) is a function of τ , then the interval process X I (t) is called as stationary interval process. Fig. 5 shows a stationary interval process, whose middle point function and bound functions are both constant with time. The interval processes which do not satisfy Definition 5 are termed as non-stationary interval process. As an effective mathematical tool for time-variant uncertainty modeling, the interval process model has been successfully applied in areas related to time-variant or dynamic uncertainty analysis of structures, such as uncertain vibration analysis [67,68], time-variant reliability analysis [66,69,70], etc. 3. The proposed interval field model In probability theory, random fields [1] are generally used to quantify the uncertainty of a spatially uncertain parameter, in which the quantity at arbitrary location is defined as a random variable with a probability distribution assigned. Different from the random field model, the interval field model employs variation bounds, namely a pair of upper and lower bounds, to describe the spatial uncertainty, which are generally much easier to obtain in practical engineering. At arbitrary location, the spatial uncertain parameter is described by an interval variable rather than a random variable. The interval field model can also be recognized as a generalization of the above interval model, from time { ⏐ process } ⏐t ∈ R to an n dimensional domain to spatial{domain; or more exactly, from one dimensional parameter space t ⏐ } parameter space u⏐u ∈ Rn . Specifically, when n = 1, the parameter vector u degenerates to a scalar, thus the interval field degenerates into an interval process. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

8

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 6. An interval field with constant bounds.

3.1. Interval field model Definition 6. A spatially uncertain parameter H (u) is an interval field if for an arbitrary spatial location uk ∈ D, k = 1, 2, . . ., the possible values of H (uk ) can be represented by an interval H I (uk ) = [H L (uk ), H U (uk )], where D is a bounded domain of the spatial location u. { } An interval field can be denoted as H (u) ∈ H I (u), u ∈ D , or H I (u) in brief. The dimension of the {interval fieldI is determined by that } of the domain D, namely n D . Fig. 6 gives a two-dimensional interval field H (u) ∈ H (u), u = (u, v) ∈ D with constant upper and lower bounds. For arbitrary location uk in the domain D, the variation of the variable H (uk ) is strictly limited within the interval H I (uk ) = [H L (uk ), H U (uk )]. Herein the two-dimensional interval field can be designated as a spatially uncertain property of a planar structure, such as the Young’s modulus of a thin plate, the temperature field of a sheet metal, etc. From the definition of the interval field we know that all the possible realizations of an interval field, depicted by the variation surfaces, will locate between the upper bound H U (u, v) and the lower bound H L (u, v). { } Definition 7. For an interval field H (u) ∈ H I (u), u ∈ D with the upper bound function H U (u) and the lower bound function H L (u), its midpoint function H m (u) and radius function H r (u) are defined as: H m (u) =

H U (u) + H L (u) , 2

H r (u) =

H U (u) − H L (u) 2

(19)

3.2. Correlation of the interval field model In many practical circumstances, dependence exists between the variables at different locations of a spatially uncertain parameter. In the following, the spatial dependence of the proposed interval field model is characterized, through which the correlation between interval variables at different locations can be specified. { } Definition 8. For an interval field H (u) ∈ H I (u), u ∈ D , the covariance function C H I (u1 , u2 ), u1 , u2 ∈ D of the interval field is defined as: C H I (u1 , u2 ) = Cov(H I (u1 ), H I (u2 )) = sin(θ ) cos(θ )(r12 − r22 )

(20)

where Cov(H I (u1 ), H I (u2 )) denotes the covariance of interval variables H I (u1 ) and H I (u2 ). The relevant angle θ , the lengths r1 and r2 of the major axis and the minor axis can be obtained based on the samples of the interval field H I (u), as shown in Fig. 7. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

9

Fig. 7. Construction of the covariance function of an interval field.

From the above definition we know that the covariance function C H I (u1 , u2 ) can reflect the correlation of the interval variables H I (u1 ) and H I (u2 ) at arbitrary locations u1 and u2 within the bounded domain D. The procedure for construction of the covariance function C H I (u1 , u2 ) in Fig. 7 can be summarized in the following steps: Step 1: Collect samples h (s) (u), s = 1, 2, . . . , Ns of the interval field H I (u), where Ns denotes the number of the samples. Step 2: For each h (s) (u), ( (s)sample (s) ) extract its values at u1 and u2 and map into the Ω -space to formulate a bivariate sample point h (u1 ), h (u2 ) . (s) Step 3: Repeat the {( ) operation of Step} 2 for all h (u), s = 1, 2, . . . , Ns , and create the set of points V = (s) (s) h (u1 ), h (u2 ) , s = 1, 2, . . . , Ns .

Step 4: Construct a minimum-area ellipse enclosing all the points in this set; Step 5: Obtain the geometric characteristics of the ellipse, namely, the relevant angle θ and the lengths of the semi-axes r1 and r2 . Step 6: Calculate the covariance function C H I (u1 , u2 ) through Eq. (20). Note that constructing a minimum-area ellipse through a point set in Step 4 can be conducted through a simple optimization procedure. For more details on how to construct the minimum-area ellipse based on the samples, the reader can refer to previous works [23,66]. { } Definition 9. For an interval field H (u) ∈ H I (u), u ∈ D with covariance function C H I (u1 , u2 ), u1 , u2 ∈ D, the correlation coefficient function R H I (u1 , u2 ) is defined as: ( ) C I (u1 , u2 ) (21) R H I (u1 , u2 ) = R H I (u1 ), H I (u2 ) = r H H (u1 )H r (u2 ) The correlation coefficient function R H I (u1 , u2 ) is a normalized indicator that specifies the degree of linear correlation of⏐ an interval field at any two spatial locations. It is not difficult to show that for arbitrary u1 , u2 ∈ D, ⏐ ⏐ R H I (u1 , u2 )⏐ ≤ 1. The larger is the absolute value of the correlation coefficient function, the stronger is the linear correlation between the interval variables H I (u1 ) and H I (u2 ). 3.3. Homogeneous and isotropic interval fields According to the characteristics of the correlation functions and the variation bounds, two typical interval fields, namely, “the homogeneous interval field” and “the isotropic interval field” can be defined. { } Definition 10. An interval field H (u) ∈ H I (u), u ∈ D with the upper bound function H U (u), lower bound function H L (u) and covariance function C H I (u1 , u2 ) is homogeneous in the direction of ∆u, if the interval field Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

10

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 8. Two-dimensional homogeneous interval field.

satisfies: H U (u) = H U (u + c∆u) H L (u) = H L (u + c∆u)

(22)

and C H I (u, u + µ∆u) = C H I (u′ , u′ + µ∆u)

(23) ′

where c ∈ R is an arbitrary real value such that u + c∆u ∈ D, u, u ∈ D are different locations within the domain D, and µ ∈ R is an arbitrary real value which satisfies u + µ∆u, u′ + µ∆u ∈ D. From Definition 10 we know that a homogeneous interval field is stationary in a specific direction. In this direction, both the variation bounds and the spatial correlation of the interval field have the property of translation invariance. Fig. 8 shows a two-dimensional interval field H (u, v) ∈ H I (u, v) that is homogeneous in the direction of u. From Fig. 8 it can been seen that although the upper and lower bounds can be fluctuant over the whole domain D, the bounds keep constant in the direction of u. Furthermore, the covariance function is only related to the distance along the direction of u, rather than depending on the positions of the two different points (u 1 , v1 ) and (u 2 , v2 ). The characteristics of the sample surfaces shown in Fig. 8 also exhibit stationarity in the direction of u, despite that they vary severely along the direction of v. In some cases, the homogeneity of an interval field is not limited to the translation invariance along a certain direction. A generalized homogeneity can thus be defined for an interval field if its variation bounds and spatial correlations are translation invariant in a subspace of the domain D. { } Definition 11. An n D -dimensional interval field H (u) ∈ H I (u), u ∈ D with the upper bound function H U (u), lower bound function H L (u) and covariance function C H I (u1 , u2 ) is homogeneous in the subspace SD = Span{∆v1 , ∆v2 , . . . ,∆vq }, q < n D of the domain D, if this interval field satisfies that ( ) q ∑ U U H (u) = H u+ ak ∆vk k=1 ( ) (24) q ∑ L L H (u) = H u+ ak ∆vk k=1

and ( CH I

u, u +

q ∑ k=1

) ak ∆vk

( = CH I

u ,u + ′

′

q ∑

) ak ∆vk

(25)

k=1

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

11

′ where points ∑qu, u ∈ D are different ∑q within the domain D, and ak ∈ R are arbitrary real values which satisfy u + k=1 ak ∆vk ∈ D and u′ + k=1 ak ∆vk ∈ D.

There exists a special type of homogeneous interval fields which requires more restrictions, namely, the isotropic interval field. An isotropic interval field has not only the translation invariance property, but also the rotation invariance. The variation bounds of an isotropic interval field are constants, and the correlation functions are only functions related to the spatial distance between u1 and u2 . { } Definition 12. An interval field H (u) ∈ H I (u), u ∈ D is isotropic if its upper and lower bound functions H U (u) and H L (u) are constants, and the covariance function C H I (u1 , u2 ) is a function only related to the spatial distance τu = ∥u1 − u2 ∥, namely, H U (u) = Uc H L (u) = L c C H I (u1 , u2 ) = C H I (τu )

(26)

where Uc and L c are two constants. Consider a two-dimensional interval field H (u, v) ∈ H I (u, v) as an example. Suppose that its upper and L lower bound functions are H U (u, v) ) respectively, and the covariance function is ( ≡√ 5 and H (u, v) ≡ 2, 2 2 C H I ((u 1 , v1 ), (u 2 , v2 )) = 2.25 exp −2 (u 1 − u 2 ) + (v1 − v2 ) , then this interval field is isotropic. For an isotropic interval field, its correlation coefficient function R H I (u1 , u2 ) is also a function related only to the spatial distance τu = ∥u1 − u2 ∥, namely, R H I (u1 , u2 ) = R H I (τu ). 4. Interval K–L expansion for the interval field model In many engineering applications, we generally hope to describe a spatially or temporally uncertain parameter with a series of variables. In probability theory, the well-known Karhunen–Lo`eve (K–L) expansion method [5,71–74] serves as a powerful tool for representation of random fields with explicitly known covariance functions. It is a series expansion method for random fields and stochastic processes, based on a spectral decomposition of the covariance function. Through the K–L expansion, a spatially continuous random field can be represented by a series of orthogonal deterministic basis functions with uncorrelated random coefficients. In general, only a few random variables are required to capture most of the uncertainty characteristics of a random field, thus immensely reducing complexity and computational cost of the subsequent uncertainty analysis. Based on the K–L expansion, in this paper a representation method of the interval field model is proposed, namely, the interval K–L expansion. Through this technique, we aim at representing the interval field accurately with multiple uncorrelated interval variables, thereby making the interval field model more convenient to use in practical uncertainty analysis. 4.1. The interval K–L expansion { } The interval K–L expansion of an interval field H (u) ∈ H I (u), u ∈ D can be expressed as: H(u) = H m (u) +

∞ ∑

√ H r (u) λ j ϕ j (u)ζ j ,

ζ j ∈ [−1, 1]

(27)

j=1

where H m (u) and H r (u) are the midpoint function and the radius function ∑of the2 interval field, respectively. ζ j ∈ ζ I = [−1, 1], j = 1, 2, . . . are uncorrelated interval variables that satisfy ∞ j=1 ζ j ≤ 1. Moreover, λ j ∈ [0, ∞) and ϕ j (u) : D → R, j = 1, 2, . . . are respectively the eigenvalues and eigenfunctions of the correlation coefficient function R H I (u1 , u2 ), which, based on Mercer’s Theorem [75], satisfy the following spectral decomposition: ∞ ∑ R H I (u1 , u2 ) = λ j ϕ j (u1 )ϕ j (u2 ) (28) j=1

From Eq. (28) it can be derived that: ∞ ∑ λ j = |D|

(29)

j=1

where |D| means the area or volume of the domain D. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

12

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

For a given interval field, the midpoint function H m (u) and the radius function H r (u) are known in advance. Therefore according to Eq. (27), the main task in formulating the interval K–L expansion for an interval field is to determine the eigenvalues λ j and the eigenfunctions ϕ j (u), j = 1, 2, . . ., which can be obtained by solving the following homogeneous Fredholm integral equation of the second kind: ∫ R H I (u1 , u2 )ϕ j (u1 )du1 = λ j ϕ j (u2 ) (30) D

∫ Herein the orthogonal eigenfunctions are normalized, i.e., D ϕi (u)ϕ j (u)du = δi j , where δi j is the Kronecker-delta function. Analytical solutions of the above Fredholm integral eigenvalue problem can be obtained only for specific types of covariance functions defined on simple domains. More details can be consulted in [76]. The Nystr¨om method [77] is provided in Section 4.3 for more general cases. { } Lemma 1. For an interval field H (u) ∈ H I (u), u ∈ D with the midpoint function H m (u), radius function Hr (u), and correlation coefficient function R H I (u1 , u2 ), the interval K–L expansion can accurately represent the upper and lower bounds, the correlation coefficient function and the covariance function of the interval field. Lemma 1 reveals that both the variation bounds and correlation functions of the interval field can be accurately quantified through the interval K–L expansion. The proof is given in Appendix A. Note that the interval K–L expansion is also applicable to the interval process X (t) ∈ X I (t): m

X (t) = X (t) +

∞ ∑

√ X r (t) λ j ϕ j (t)ζ j ,

ζ j ∈ [−1, 1]

(31)

j=1

in which X m (t) and X r (t) are the midpoint function and the radius function of the interval process, respectively. 4.2. Truncated interval K–L expansion and error analysis The interval K–L expansion allows us to represent the continuous spatial uncertainty of an interval field as a superposition of infinite deterministic functions with uncorrelated interval coefficients. However, it is not only impossible but also unnecessary to use an infinite number of terms for representation of this type of uncertainty. Generally, most of the characteristics of a spatially uncertain parameter can be reflected considerably by those principle terms. In the interval K–L expansion, the contribution of a term to uncertainty representation is indicated by the eigenvalue. Therefore for practical implementation, it can be approximated by sorting the eigenvalues λi in a descending order and truncating the expansion after M terms. Thus a truncated interval K–L expansion can be expressed as: ˜ H(u) = H m (u) +

M ∑

√ H r (u) λ j ϕ j (u)ζ j

(32)

j=1

{ The truncated form with M terms generates an error in the representation of the interval field H (u) ∈ H I (u), } u ∈ D . In the following error analysis, two kinds of error indexes for the truncated interval K–L expansion are given. Lemma 2. The error εm (u) of the midpoint function and the average relative error εr 2 (u) of squared radius function satisfy: ˜ m (u) ≡ 0 εm (u) = H m (u) − H

(33)

( )2 ∑M ∫ (H r (u))2 − H˜ r (u) M 1 1 ∑ j=1 λ j ∑ εr 2 (u) = du = 1 − = 1 − λj ∞ |D| D |D| (H r (u))2 j=1 λ j

(34)

j=1

˜ m (u) and H ˜ r (u) are respectively the midpoint function and the radius function of the truncated interval where H { } ˜ K–L expansion H(u) of the interval field H (u) ∈ H I (u), u ∈ D . Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

13

The proof is given in Appendix B. Lemma 2 points out that the truncated form of the interval K–L expansion does not produce an error in the midpoint function representation of the interval field. Error arises in the radius function. The average relative error εr 2 (u) over the domain D depends on the number M of terms in the truncated form. The more terms are retained in the expansion in Eq. (32), the higher precision is reached for the radius function. ˜ over the domain D satisfies that: Lemma 3. The average square error ε 2h (u) of h(u) ∫ ( ∞ M )2 1 1 ∑ 1 ∑ 2 2 ˜ ε h (u) = h(u) − h(u) du = λjζj ≤ 1 − λj (35) |D| D |D| |D| j=M+1 j=1 { } ˜ where h(u) and h(u) are the standardized forms of the interval field H (u) ∈ H I (u), u ∈ D and its truncation ˜ H(u), respectively. Namely: H(u) − H m (u) (36) h(u) = H r (u) ˜ H(u) − H m (u) ˜ h(u) = (37) H r (u) The proof is given in Appendix C. Lemmas 2 and 3 provide two kinds of spatial error indexes for the truncated interval K–L expansion. Based on the above error indexes (Eqs. (34) and (35)), the approximation degree κ is defined ˜ to degree of} approximation of the truncated interval K–L expansion H(u) for the original interval field { indicate the H (u) ∈ H I (u), u ∈ D , namely: M

κ=

1 ∑ λj |D|

(38)

j=1

For a specific spatially uncertain parameter, if more terms in the interval K–L expansion are included, the approximation degree κ approaches to 1. Ideally, in the limit M → ∞, the approximation degree reaches 1, namely 1 ∑∞ κ∞ = |D| j=1 λ j = 1. 4.3. Nyström method for the interval K–L expansion A continuous interval field can be represented by a series of orthogonal functions with interval coefficients. However, the Fredholm integral equation in Eq. (30) can be solved analytically only for a few kernel functions and geometries of D. In most cases, numerical methods such as the Nystr¨om method [77,78], the collocation method [79] or Galerkin method [80] are required. In this work, the Nystr¨om method is adopted as the numerical integration scheme for approximation of the Fredholm integral eigenvalue problem. The truncated interval K–L expansion in Eq. (32) can be approximated as: √ M ∑ ˆ H(u) = H m (u) + H r (u) λˆ j ϕˆ j (u)ζ j (39) j=1

where λˆ j and ϕˆ j (u) are approximations to the true eigenvalues λ j and eigenfunctions ϕ j (u). In the Nystr¨om method [77], the eigenvalue problem of Eq. (30) can be approximated as: N ∑

wi R H I (ui , u)ϕˆ j (ui ) = λˆ j ϕˆ j (u)

(40)

i=1

where ui ∈ D with i ∈ {1, 2, . . . , N }, represent a set of N ∈ N integration points, and wi is the integration weight associated with ui . Solving the problem of Eq. (40) at the integration points leads to the following group of equations: N ∑

wi R H I (ui , uk )ϕˆ j (ui ) = λˆ j ϕˆ j (uk ),

k = 1, 2, . . . , N

(41)

i=1 Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

14

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

which can be formulated in matrix notation as: RWy j = λˆ j y j

(42)

where R is a symmetric positive semi-definite N × N matrix with elements rki = R H I (ui , uk ), W is a diagonal matrix of size N with entries Wii = wi , and y j is an N-dimensional vector whose kth entry is y j,k = ϕˆ j (uk ). √ By substituting y ∗j,k = wk y j,k into Eq. (42), and denote y j = W−1/2 y∗j where W−1/2 denotes the inverse of the √ 1/2 diagonal matrix W1/2 whose entries are Wii = wi , thus the matrix eigenvalue problem can be reformulated to an equivalent matrix eigenvalue problem: R∗ y∗j = λˆ j y∗j

(43)

in which R∗ = W1/2 RW1/2 is a symmetric positive semi-definite matrix, and y∗j is the vector whose kth entry is y ∗j,k . The eigenvalues λˆ j and the eigenvectors y∗j , j = 1, 2, . . . , M thus can be derived, by which the values of the eigenfunctions at the integration points, namely ϕˆ j (ui ), i ∈ {1, 2, . . . , N }, can be then determined. Based on Eq. (40) we have: ϕˆ j (u) =

N 1 ∑ wi R H I (ui , u)ϕˆ j (ui ) λˆ j

(44)

i=1

Taking into account that ϕˆ j (ui ) = then expressed as: ϕˆ j (u) =

√1 wi

y ∗j,i , the Nystr¨om interpolation formula of the eigenfunction ϕˆ j (u) can be

N 1 ∑√ wi y ∗j,i R H I (ui , u) ˆλ j

(45)

i=1

By substituting Eq. (45) into Eq. (39), the truncated interval K–L expansion can be finally approximated as: ⎛ ⎞ M N ∑ ∑ ζ √ ⎝ H r (u) √ j ˆ H(u) = H m (u) + wi y ∗j,i R H I (ui , u)⎠ (46) j=1 λˆ j i=1 ∫ The eigenfunctions have to be normalized such that D ϕˆi (u)ϕˆ j (u)du = δi j . By numerical integration and using ∫ the same integration points and integration weights as the ones used in Eq. (41), the inner product D ϕˆi (u)ϕˆ j (u)du can be then approximated as: N ∑

( )T wk ϕˆi (u)ϕˆ j (u) ≈ yi∗ y∗j

(47)

k=1

    Therefore, it is required that the eigenvectors y∗j , j = 1, 2, . . . , M are orthogonal and normalized with y∗j  = 1 2 for all j. The integration points ui , i = 1, 2, . . . , N in the above formulation can be adopted as Gaussian quadrature points or uniformly distributed points, etc. If these integration points distribute uniformly over the domain D, the matrix W in Eq. (42) can be then written as W = wI, where I is the identity matrix and w = |D| N . In addition, note that herein the Nystr¨om method is employed for approximation of the eigenvalues and eigenfunctions in the Fredholm integral equation in Eq. (30). Generally, other numerical methods such as the collocation method [79], the Galerkin method [80], or the finite cell method [81], etc. can also be applied. 4.4. Simulation of the interval field model Based on the interval K–L expansion and its truncated form in Eq. (32), the simulation of an interval field H I (u) can be achieved through sampling only in the M-dimensional space of the interval variables ζ j ∈ ζ jI = [−1, 1], j = ∑ 2 M 1, 2, . . . , M with the condition M j=1 ζ j ≤ 1. Denote Ωunit as the M-dimensional unit sphere: ⎧ ⏐ ⎫ ⏐∑ ⎨ ⎬ ⏐ M 2 M Ωunit = ζ ∈ R M ⏐⏐ ζj ≤ 1 (48) ⎩ ⎭ ⏐ j=1 Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

15

Thus the simulation of the proposed interval field is converted to generating samples in the M-dimensional unit M sphere Ωunit and then mapping into the original space to obtain the realizations of the interval field H I (u). The detailed procedure can be summarized by the following steps. r Step 1: Obtain the midpoint function H m (u), the radius } function H (u), the correlation coefficient function { I R H I (u1 , u2 ) of the interval field H (u) ∈ H (u), u ∈ D .

Step 2: Set the required value κ of the approximation degree κ. Step 3: Solve the Fredholm equation in Eq. (30), obtain the largest M eigenvalues λ j , j = 1, 2, . . . , M that satisfy 1 ∑M κ = |D| j=1 λ j > κ and corresponding eigenfunctions ϕ j (u), j = 1, 2, . . . , M. M Step 4: Generate samples ζ = [ζ1 , ζ2 , · · · , ζ M ]T from the M-dimensional unit sphere Ωunit . The details on how to generate the samples of ζ are given in Appendix D.

Step 5: Use the samples of ζ to generate realizations of the interval field H I (u) with the truncated interval K–L expansion in Eq. (32). { In }the following, the simulation of a two-dimensional homogeneous interval field H (u, v) ∈ H I (u, v), (u, v) ∈ D is given as an example. The upper and lower bound functions are set as constants, namely, H U (u, v) ≡ 103 and H L (u, v) ≡ 97. The correlation coefficient function is: ( ) |u 1 − u 2 | |v1 − v2 | R H I ((u 1 , v1 ), (u 2 , v2 )) = exp − − (49) lu lv in which lu and lv are correlation lengths in the directions of coordinates u and v, respectively. The bounded domain D of the two-dimensional interval field is D = {(u, v) |u ∈ [0, 6] , v ∈ [0, 8]}. Fig. 9 gives three sample realizations of the two-dimensional interval field when the correlation lengths take different values. The correlation lengths are first set with relatively large values, i.e., lu = lv = 32, which implies a strong degree of correlation between the variables at adjacent spatial locations. In other words, the values of the variables that nearby in location are strongly related to each other. Consequently, the sample realization shown in Fig. 9(a) present relatively smooth surfaces. While in Fig. 9(b), the sample realization is coupled to the interval field with correlation lengths lu = lv = 8, thus exhibiting rough surfaces such that the correlation lengths indicate relatively weak relationship between variables at different locations. By comparing Fig. 9(a) with (b), it can be verified that for a certain type of correlation coefficient function, a larger correlation length implies a stronger spatial correlation degree for a spatially uncertain parameter, whose samples are more smooth. On the contrary, a smaller correlation length means a weaker spatial correlation degree, thus the samples of the interval field vary more severely. Fig. 9(c) shows the sample realization of the interval field with correlation lengths lu = 32, lv = 4. In the coordinate of u, the correlation length lu = 32 indicates a strong correlation; while in the coordinate of v, the correlation length lv = 4 means a relatively weak correlation. Therefore the sample surface shown in Fig. 9(c) exhibits relatively large variation along the direction of coordinate of v, while keeps smooth along the direction of coordinate of u. This obvious difference of correlation lengths in different directions makes the sample surface exhibits quite different smoothness along these directions. Although the correlation lengths can have significant impacts on the shapes of the sample realizations of the interval field, any of them will definitely limited within the predetermined upper and lower bounds, which agrees well with the definition of the proposed interval field model. 5. The interval finite element methods The finite element method (FEM) [82–84] has been widely applied in various fields of engineering. By combining the FEM and the proposed interval field model, the interval finite element analysis is carried out for structural spatial uncertainty analysis. The first order and the second order perturbation-based interval finite element methods (P-IFEMs) are developed, by which response intervals of the spatially uncertain structures can be obtained. As a more generic method, the Monte Carlo simulation (MCS) method is also employed for evaluation of the upper and lower response bounds, which can be referred to for verification of the other IFEMs. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

16

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 9. Sample realizations of the two-dimensional homogeneous interval field with different correlation lengths.

5.1. Formulation of the interval equilibrium equations In the linear FEM, the whole domain A of a structure is partitioned into finite elements. The stiffness equations of an element are constructed as: Ke de = pe

(50)

e

e

where d is the vector of element’s nodal displacements, p is the vector of equivalent nodal forces that applied to the element. The element stiffness matrix Ke is derived by the following integration: ∫ Ke = BT DBdA (51) Ae

where A denotes the whole domain of the structure, Ae is the integral domain of the element, B is the strain– displacement matrix, D is the elasticity matrix. Assembling all sets of element stiffness equations into a global system of equations, the global equilibrium equations can be then formulated. When spatial uncertainties such as the uncertain material properties are involved in the finite element analysis, the equilibrium equations are no longer deterministic and cannot be solved directly. Denote α(u) as a spatially uncertain parameter, the variation of which does not lead to the change of the shape of the element, thus the element stiffness matrix Ke becomes non-deterministic and can be expressed as: ∫ e K (α(u)) = BT D(α(u))BdA (52) Ae

where D is non-deterministic due the involvement of the spatial uncertain parameter α(u). Consequently, the global stiffness matrix K that is assembled from the element stiffness matrices is also a non-deterministic matrix that Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

17

relates to α(u): K(α(u)) =

Ne ∑

GiT Kei (α(u))Gi

(53)

i=1

where Σ stands for the matrix assembling rather than the sum operation, Ne is the number of the elements, Kei is the stiffness matrix of the ith element, and Gi is the transition matrix for the assembling corresponding to Kei . Thus the resulting equilibrium equations finally turn to be: K(α)d = p

(54)

in which K(α) is the non-deterministic global stiffness matrix regarding the spatial uncertain parameter α(u), d is the vector of global nodal displacements, p is the vector of global nodal forces. In this paper, we only deal with the cases where the material property is spatially uncertain and quantified by the proposed interval field model, while the external loads p are deterministic. The adoption of interval field model makes the stiffness matrix K(α) an interval matrix, consequently the displacement response d also forms an interval vector. The above non-deterministic equilibrium equations in (54) are a type of well-known interval linear equations [15,16], the solution set of which generally has a very complicated form and is very difficult to achieve analytically. In most cases, the solution of interest is to estimate the upper and lower bounds for responses such as the displacement or stress at each node of the structure. In the following, based on the perturbation theory [85,86], the first order and the second order P-IFEMs are developed for structural response bounds analysis. 5.2. Perturbation-based interval finite element methods 5.2.1. The first order P-IFEM (a) Displacement analysis Using the truncated interval K–L expansion, the spatial uncertain parameter α(u) represented by the interval field M model can be substituted by M uncorrelated interval variables ζ j , j = 1, 2, . . . , M which satisfy ζ ∈ Ωunit , namely: α(u) = α m (u) +

M ∑

√ αr (u) λ j ϕ j (u)ζ j

(55)

j=1

where α m (u) and αr (u) are respectively the midpoint function and the radius function of α(u). The first order Taylor series expansion of the global stiffness matrix K(α) and the resulting displacement vector d around their midpoint values can be expressed as: ⏐ M ∑ ∂K(α) ⏐⏐ K(α) ≈ K0 + ζj (56) ∂ζ j ⏐ζ=0 j=1 d(α) ≈ d0 +

⏐ M ∑ ∂d(α) ⏐⏐ ζ j , j = 1, 2, . . . , M ∂ζ j ⏐ζ=0 j=1

(57)

where K0 and d0 are respectively the global stiffness matrix and ⏐ ⏐ displacement vector with the uncertain parameter ∂K(α) ⏐ ∂d(α) ⏐ α(u) taking as its midpoint value function. ∂ζ j ⏐ and ∂ζ j ⏐ represent the derivatives of K(α) and d(α) with ζ=0

ζ=0

respect to ζ j evaluated at ζ = 0, respectively. Substituting Eqs. (56) and (57) into Eq. (54) yields: ⎞⎛ ⎞ ⎛ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂K(α) ⏐ ∂d(α) ⏐ ⎝K0 + ⎠ ⎝d0 + ζj⎠ = p ⏐ ζj ∂ζ ∂ζ j ⏐ζ=0 j ζ=0 j=1 j=1

(58)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

18

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

namely, ⎛ ⎞ ⎛ ⎞⎛ ⎞ ⏐ ⏐ ⏐ M M M ∑ ∑ ∑ ⏐ ⏐ ⏐ ∂K(α) ⏐ ∂K(α) ⏐ ∂d(α) ⏐ ⎠ ⎝ ⎠⎝ ζj⎠ K0 d0 + ⎝ ⏐ ζ j d0 + ⏐ ζj ∂ζ ∂ζ ∂ζ j ⏐ζ=0 j j ζ=0 ζ=0 j=1 j=1 j=1 ⎛ ⎞ ⏐ M ∑ ⏐ ∂d(α) ⏐ = p − K0 ⎝ ζj⎠ ∂ζ j ⏐ζ=0 j=1

(59)

By ignoring the second order terms and identifying the coefficients of ζ j on both sides, one can obtain: d0 = K−1 0 p

(60)

⏐ ⏐ ⏐ ∂d(α) ⏐⏐ −1 ∂K(α) ⏐ = −K0 d0 ∂ζ j ⏐ζ=0 ∂ζ j ⏐ζ=0

(61)

Substituting Eq. (61) into Eq. (57), the displacement vector d(α) can be approximated as: ( ) ⏐ M ∑ ⏐ ∂K(α) ⏐ d0 ζ j K−1 d(α) ≈ d0 − 0 ⏐ ∂ζ j ζ=0 j=1

(62)

[ ]T M of the non-deterministic stiffness where ζ = ζ1 , ζ2 , . . . , ζ M satisfies that ζ ∈ Ωunit . The gradients ∂K(α) ∂ζ j matrix K(α) with respect to uncorrelated interval variables ζ j , j = 1, 2, . . . , M can be computed by the following equations: N

e ∂Ke (α) ∂K(α) ∑ G iT i = Gi ∂ζ j ∂ζ j i=1 ∫ ∂Kie (α) ∂D(α(u)) BT = BdA e ∂ζ j ∂ζ j A

(63) (64)

∂D(α(u)) ∂α(u) ∂D(α(u)) = ∂ζ j ∂α ∂ζ j √ ∂α(u) = αr (u) λ j ϕ j (u) ∂ζ j

(65) (66) ∂Kie (α) is ∂ζ j ∂D(α(u)) of the ∂ζ j

where Σ represents the matrix assembling, and

the gradient of the ith element stiffness matrix with respect

elasticity matrix. As given in Eq. (66), the gradient ∂α(u) of to ζ j , which is determined by the gradient ∂ζ j the interval field α(u) can be easily derived from the truncated interval K–L expansion in Eq. (55). For convenience, we denote: ⏐ ∂K(α) ⏐⏐ q j = K−1 d0 (67) 0 ∂ζ j ⏐ζ=0 which is a deterministic quantity, and the displacement vector in Eq. (62) can be then rewritten as: d(α) ≈ d0 −

M ∑

qjζj

(68)

j=1

Now, we denote n dof as the number of degrees of freedom in the finite element model, then the upper bound dkU and the lower bound dkL of the kth (k = 1, 2, . . . , n dof ) degree of freedom of the displacement response can be obtained by solving the following set of convex optimization problems: ⎧ M ⎪ ⎨d U = max d − ∑ q ζ 0,k j,k j k (69) ζ j=1 ⎪ ⎩ M s.t. ζ ∈ Ωunit Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

19

and ⎧ M ⎪ ⎨d L = min d − ∑ q ζ 0,k j,k j k ζ j=1 ⎪ ⎩ M s.t. ζ ∈ Ωunit

(70)

in which d0,k is the value of the kth component in d0 , and q j,k is the value of the kth component in q j . The above convex optimization problems can be analytically solved by using Lagrange multipliers [87,88], which leads to:   M ∑ U (71) dk = d0,k + √ q 2j,k j=1

dkL = d0,k

  M ∑ −√ q2

(72)

j,k

j=1

With the displacement response intervals [dkL , dkU ], k = 1, 2, . . . , n dof of all degrees of freedom computed, the upper and lower bounds of the spatially uncertain response of the structure are finally obtained. (b) Stress analysis The stress vector σe of an element can be obtained by: σe = D(α(u))Bde

(73)

in which σe = [σx x , σ yy , τx y ]T (for plane-stress/strain problems) or σe = [σx x , σ yy , σzz , τx y , τ yz , τzx ]T (for space problems) is the stress vector of the element e. The vector de of element nodal displacements can also be denoted by: de = Ie d

(74)

where Ie is a diagonal matrix with only the quantities corresponding to the nodes of the element e equal to 1 and others 0. Then the stress vector σe of the element can be rewritten as: σe = D(α(u))BIe d

(75)

As discussed previously, the nodal displacement vector d can be approximated by Eq. (62). Hence the stress vector of each element can be denoted as: ⎛ ⎞ ⏐ M ∑ ⏐ ∂K(α) ⏐ d0 ζ j ⎠ ⎝p − σe = D(α(u))BIe K−1 (76) 0 ⏐ ∂ζ j ζ=0 j=1 Now, using a first order Taylor expansion of the elasticity matrix D(α(u)) we have: ⎛ ⎞ ⎛ ⎞ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂D(α) ∂K(α) ⏐ ζ j ⎠ BIe K−1 ⎝p − ⏐ d0 ζ j ⎠ σe ≈ ⎝D0 + 0 ∂ζ ⏐ ∂ζ ⏐ j=1

⎛

j

ζ=0

j=1

j

ζ=0

⎞ ⏐ ⏐ M ∑ ⏐ ⏐ ∂D(α) ∂K(α) ⏐ ζ j ⎠ BIe ⎝d0 − K−1 ⏐ d0 ζ j ⎠ = ⎝D0 + 0 ⏐ ⏐ ∂ζ ∂ζ j j ζ=0 ζ=0 j=1 j=1 ⎛ ⎞ ⎞ ⎛ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂D(α) ⏐ ∂K(α) ⏐ e −1 ⎝ ⎠ e d0 ζ j ⎠ − ∆ζ 2 = D0 BIe d0 + ⎝ ⏐ ζ j BI d0 − D0 BI K0 ∂ζ ∂ζ j ⏐ζ=0 j ζ=0 j=1 j=1 M ∑

⎞

⎛

(77)

where ∆ζ 2

⎞ ⎞ ⎛ ⎛ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂D(α) ⏐ ∂K(α) ⏐ ⎠ e −1 ⎝ =⎝ d0 ζ j ⎠ ⏐ ζ j BI K0 ∂ζ ∂ζ j ⏐ζ=0 j ζ=0 j=1 j=1

(78)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

20

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

thus formulating an analytical expression of the interval variables ζ j , j = 1, 2, . . . , M. If the second order terms, namely ∆ζ 2 were furtherly neglected, the stress vector then degenerates into a group of linear functions of ζ j . Similar to the convex optimization problems for nodal displacements, the upper and lower bounds of the stress vector approximated by Eq. (77) can be conveniently solved by the method of Lagrange multipliers. 5.2.2. The second order P-IFEM (a) Displacement analysis The second order Taylor series expansion of the global stiffness matrix K(α) and the resulting displacement vector d can be expressed as: ⏐ ⏐ M M M ∑ 1 ∑ ∑ ∂K2 (α) ⏐⏐ ∂K(α) ⏐⏐ ζj + ζi ζ j (79) K(α) ≈ K0 + ∂ζ j ⏐ζ=0 2 i=1 j=1 ∂ζi ∂ζ j ⏐ζ=0 j=1 and ⏐ ⏐ M M M ∑ ∂d(α) ⏐⏐ 1 ∑ ∑ ∂d2 (α) ⏐⏐ d(α) ≈ d0 + ζj + ζi ζ j ∂ζ j ⏐ζ=0 2 i=1 j=1 ∂ζi ∂ζ j ⏐ζ=0 j=1 where

⏐

∂K2 (α) ⏐ ∂ζi ∂ζ j ⏐ζ=0

and

⏐

∂d2 (α) ⏐ ∂ζi ∂ζ j ⏐ζ=0

represent the second derivatives of K(α) and d(α) with respect to ζi and ζ j evaluated

at ζ = 0, respectively. Substituting Eq. (79) into Eq. (54) yields: ⎛ ⎞ ⏐ ⏐ M M ∑ M 2 ∑ ∑ ⏐ ⏐ ∂K(α) ⏐ ∂K (α) ⏐ 1 ⎝K0 + ζi ζ j ⎠ ⏐ ζj + 2 ∂ζ ∂ζi ∂ζ j ⏐ζ=0 j ζ=0 j=1 i=1 j=1 ⎞ ⎛ ⏐ ⏐ M M ∑ M 2 ∑ ∑ ⏐ ⏐ ∂d (α) ∂d(α) 1 ⏐ ζj + ⏐ ζi ζ j ⎠ = p · ⎝d0 + ∂ζ ⏐ 2 ∂ζ ∂ζ ⏐ j=1

j

(80)

ζ=0

i=1 j=1

i

(81)

j ζ=0

By ignoring the third and higher order terms in the expansion formulation of the above equation and identifying the coefficients of ζ j and ζi ζ j on both sides, Eqs. (60) and (61) as well as the following equation can be obtained: ( ⏐ ⏐ ⏐ ⏐ ) ∂K2 (α) ⏐⏐ ∂K(α) ∂d(α) ⏐⏐ ∂K(α) ∂d(α) ⏐⏐ ∂d2 (α) ⏐⏐ −1 = −K0 d0 + + (82) ∂ζi ∂ζ j ⏐ζ=0 ∂ζi ∂ζ j ⏐ζ=0 ∂ζi ∂ζ j ⏐ζ=0 ∂ζ j ∂ζi ⏐ζ=0 in which the Hessian matrix

∂K2 (α) ∂ζi ∂ζ j

of K(α) with respect to the uncorrelated interval variables ζi and ζ j , i, j =

1, 2, . . . , M can be computed by the following equations: ( e ) ∂Kk (α) Ne 2 ∂ ∑ ∂ζi ∂K (α) T = Gk Gk ∂ζi ∂ζ j ∂ζ j k=1 ( e ) ∂Kk (α) ∫ ∂ ∂ζi ∂D2 (α(u)) = BT BdA ∂ζ j ∂ζi ∂ζ j Ae ∂D2 (α(u)) ∂D2 (α(u)) ∂α(u) ∂α(u) = ∂ζi ∂ζ j ∂α 2 ∂ζi ∂ζ j

(83)

(84) (85)

By substituting Eqs. (60), (61) and (82) into Eq. (80), the second order approximation of the displacement response vector can be then determined. For each degree of freedom, the displacement is expressed as a second-order M polynomial of the interval variables ζ j , j = 1, 2, . . . , M that satisfy ζ ∈ Ωunit , thus the upper and lower bounds can be conveniently solved by optimization algorithms, such as the SQP algorithm [89]. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

21

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

(b) Stress analysis From Eqs. (75) and (80), the stress vector of each element can be expressed as: ⎛ ⎞ ⏐ ⏐ M M ∑ M 2 ∑ ∑ ⏐ ⏐ 1 ∂d(α) ∂d (α) ⏐ ζj + ⏐ ζi ζ j ⎠ σe = D(α(u))BIe ⎝d0 + ∂ζ ⏐ 2 ∂ζ ∂ζ ⏐ j

j=1

ζ=0

i=1 j=1

i

By second order Taylor expansion of the elasticity matrix D(α(u)) we have: ⎛ ⎞ ⏐ ⏐ M M ∑ M 2 ∑ ∑ ⏐ ⏐ 1 ∂D(α) ∂D (α) e ⏐ ζj+ ⏐ ζi ζ j ⎠ BIe σ ≈ ⎝D0 + ⏐ ⏐ ∂ζ 2 ∂ζ ∂ζ j i j ζ=0 ζ=0 j=1 i=1 j=1 ⎛⎛ ⎞⎞ ⏐ ⏐ M M ∑ M 2 ∑ ∑ ⏐ ⏐ 1 ∂d(α) ∂d (α) ⏐ ζj + ⏐ ζi ζ j ⎠⎠ · ⎝⎝d0 + ∂ζ j ⏐ζ=0 2 i=1 j=1 ∂ζi ∂ζ j ⏐ζ=0 j=1 ⎛ ⎞ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂D(α) ∂d(α) ⏐ BIe d0 + D0 BIe ⏐ ⎠ ζj = D0 BIe d0 + ⎝ ∂ζ ⏐ ∂ζ ⏐ j=1

j

ζ=0

j

j=1

(86)

j ζ=0

(87)

ζ=0

⏐ ⏐ ⏐ ⏐ ) M ∑ M 2 ∑ ⏐ ⏐ ⏐ ∂D(α) 1 1 ∂D2 (α) ⏐⏐ ∂d(α) ∂d (α) ⏐ BIe ⏐ ⏐ BIe d0 + D0 BIe ζi ζ j + ∆ζ 3 ,ζ 4 + ⏐ ⏐ ⏐ ⏐ 2 ∂ζ ∂ζ ∂ζ ∂ζ 2 ∂ζ ∂ζ i j i j i j ζ=0 ζ=0 ζ=0 ζ=0 i=1 j=1 (

where ∆ζ 3 ,ζ 4 denotes the term containing the third and fourth order terms of ζ j , j = 1, 2, . . . , M. By substituting Eqs. (61) and (82) into Eq. (87) and ignoring ∆ζ 3 ,ζ 4 , we finally obtain the second order approximation of the stress vector σe of each element: ⎛ ⎞ ⏐ ⏐ M M ∑ ∑ ⏐ ⏐ ∂D(α) ∂K(α) ⏐ BIe d0 − D0 BIe K−1 ⏐ d0 ⎠ ζ j σe ≈ D0 BIe d0 + ⎝ 0 ∂ζ ⏐ ∂ζ ⏐ j=1

j

ζ=0

j=1

j

ζ=0

) ⏐ ⏐ ⏐ M ∑ M ∑ ⏐ ∂D(α) ⏐⏐ 1 ∂D2 (α) ⏐⏐ −1 ∂K(α) ⏐ e e BI d0 + d0 ζi ζ j BI − K0 + 2 ∂ζi ∂ζ j ⏐ζ=0 ∂ζi ⏐ζ=0 ∂ζ j ⏐ζ=0 i=1 j=1 (( ⏐ ⏐ ⏐ ) ) M ∑ M ∑ 1 ∂K2 (α) ⏐⏐ ∂K(α) −1 ∂K(α) ⏐⏐ ∂K(α) −1 ∂K(α) ⏐⏐ e −1 − D0 BI K0 − K0 − K0 d0 ζi ζ j 2 ∂ζi ∂ζ j ⏐ζ=0 ∂ζi ∂ζ j ⏐ζ=0 ∂ζ j ∂ζi ⏐ζ=0 i=1 j=1 (

(88) Analogously to the computation of the upper and lower bounds for the displacement response, the stress response bounds of each element can be approximated using the nonlinear optimization methods. Based on the constitutive relation, the strain uncertainty analysis can also be carried out in a similar way, through which the upper and lower bounds of the structural strain can be obtained. 5.3. Monte Carlo simulation method for interval finite element analysis Based on the simulation method of the interval field, the Monte Carlo simulation (MCS) method can be employed to evaluate the upper and lower response bounds of structures under spatial uncertainties. The MCS method can generate accurate response bounds when the number of input samples is enough. Therefore, it is presented in the following to provide a reference solution for the P-IFEMs. The MCS method can be implemented in three steps. Firstly, simulate a population of sample realizations of the spatially uncertain parameter based on the interval field model. Secondly, perform the finite element analysis associated with each member of that population. Thirdly, generate a set of structural responses. Finally, compute the Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

22

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

upper and lower bounds of the responses. The detailed implementation procedure of the MCS method for response bounds analysis is given in the following: Step 1: Generate a population of Ns sample realizations α (s) (u), s = 1, 2, . . . , Ns of the interval field parameter α(u) through the simulation method in Section 4.4. Step 2: {For each sample α (s) (u), perform a deterministic finite element analysis, to obtain a population of response } vectors d(s) (u), s = 1, 2, . . . , Ns . { Step 3: For each degree of freedom, compute the maximum and the minimum values of the responses d(s) (u), s = 1, 2, . . . , Ns }, and obtain the upper bound dU (u) and the lower bound d L (u) of the spatial response. 6. Numerical examples In this section, two numerical examples are investigated. First, the displacement response analysis of a square plate with spatially uncertain Young’s modulus is carried out by the P-IFEMs, and the accuracy of the results is verified by the MCS method. Second, the first order P-IFEM is applied to stress response bounds analysis of a plate with a hole. 6.1. A square thin plate A square thin plate subjected to distributed forces at the two sides is shown in Fig. 10(a). The side length of the quadrate plate is a = 1 m, and the thickness is h = 0.0254 m. The sides along the thickness direction of vertices A and B are fixed. The distributed forces pleft = pright = 16 MPa are equal in value but opposite in direction. The Poisson’s ratio of the concrete plate is ν = 0.3. The Young’s modulus is spatially uncertain, which is described as an interval field E(u) ∈ E I (u). It is treated as a plane stress problem. The finite element model is employed with 100 four-node elements, as shown in Fig. 10(b). The displacement of the plate will also be spatially uncertain due to the uncertainty of the Young’s modulus. In this example, the displacement analysis is carried out by the proposed P-IFEMs, and verified by the MCS method. First, the upper and lower bounds of displacement of the plate is given by the first order P-IFEM. Second, the influence of the number of truncated terms in interval K–L expansion on the response bounds is investigated. Finally, the accuracy of the P-IFEMs is verified by the MCS method. (a) Displacement response bounds analysis by the first order P-IFEM The Young’s modulus is represented by a two-dimensional homogeneous interval field with the constant midpoint function E m (u) = 32.5 GPa and radius function E r (u) = 0.2E m (u). The correlation coefficient function R E I (u1 , u2 ) is: ( ) |u 1 − u 2 | |v1 − v2 | − (89) R E I (u1 , u2 ) = exp − lu lv in which lu = 2.4 m and lv = 3.2 m are the correlation lengths. By using the truncated interval K–L expansion, the spatially uncertain Young’s modulus E I (u) is approximated by 24 uncorrelated interval variables, namely: E(u) ≈ E m (u) + E r (u)

24 ∑ √

λ j ϕ j (u)ζ j

(90)

j=1

} { ⏐∑ ⏐ 24 2 with ζ ∈ Ωunit = ζ ⏐ 24 ζ ≤ 1 . The approximation degree κ reaches 98.18%. j=1 j The response bounds of displacement in the directions of coordinates u and v given by the first order P-IFEM are respectively shown in Figs. 11 and 12. From the results we can see that due to the uncertainty of the Young’s modulus E(u), the displacements are also spatially uncertain. The region enveloped by the upper bound and the lower bound indicates the variation range of displacement of the plate. For example, at the node in A′ , the displacement du in the direction of coordinate u belongs to [2.892, 4.384]×10−4 m, and the displacement dv in the direction of coordinate of v belongs to [0.747, 1.362] × 10−4 m. Because the vertices A and B are fixed, their displacement bounds are both degenerated into the value zero. (b) Influence of the number of truncated terms on response bounds Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

23

Fig. 10. A concrete plate subjected to distributed forces.

Since a larger number M of the truncated terms in the interval K–L expansion means a better approximation of the spatial uncertainty of the Young’s modulus, thus the computed response bounds are more accurate. To grasp the influences of the number of truncated terms on the responses of the plate, the displacement response bounds at locations of vertex A′ and the center O are provided with increasing M. Fig. 13 shows the displacement response bounds of A′ and O in the direction of u, respectively, and Fig. 14 shows their response bounds in the direction of v. The first order P-IFEM and the MCS method are separately analyzed. From the results it can be observed that both the response bounds given by the first order P-IFEM and the MCS method tend to be steady with increasing M. From Fig. 13 it can be found that in this case the approximation degree of E(u) ∈ E I (u) has little impact on the displacement response bound of A′ in the direction of u; when M is greater than 6, the response bounds by both the P-IFEM and the MCS method become stable. While for the response bounds at the center point O, it changes greatly when M is smaller than 6 and keeps steady when M is greater than 6. Similarly, from Fig. 14 it can be observed that the displacement response bounds of both the vertex A′ and the center O in the direction of v keep steady when M reaches 5. Therefore, it can be approximately recognized that only the most prominent 6 terms are necessary in the interval K–L expansion for evaluation of the displacement response bounds for this plate with uncertain Young’s modulus. If more terms are retained, the Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

24

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 11. Response bounds of displacement du by the first order P-IFEM.

Fig. 12. Response bounds of displacement dv by the first order P-IFEM.

approximation of the spatially uncertain Young’s modulus that quantified by the interval field model will be more accurate; however, this might have little effects on the precision of the response bounds. (c) Results comparison of the P-IFEMs and the MCS method With the number M of truncated terms taken as 6, the displacement response bounds of the plate with spatially uncertain Young’s modulus computed by the first order P-IFEM and the MCS method are compared. Here 1 million Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

25

Fig. 13. Displacement response bounds of A′ and O in the direction of u.

Fig. 14. Displacement response bounds of A′ and O in the direction of v.

samples of the interval field E(u) ∈ E I (u) are used in the MCS method. Figs. 15 and 16 present respectively the displacement response bounds in the directions of u and v. And the response intervals at nodes marked in circles in Fig. 10(b) are also given in Tables 1 and 2, in which the response intervals given by the first order P-IFEM and the MCS method are compared. The relative errors are also provided in Figs. 17 and 18. From the relative errors of the upper bounds and the lower bounds, it can be observed that the first order P-IFEM presents with good precision for evaluation of the upper and lower bounds of displacement du . All the relative errors are controlled within 9%, as shown in Table 1 and Fig. 17. From Table 2 and Fig. 18 we can find that at most nodes the displacement bounds of dv given by the first order P-IFEM show satisfactory precision, while for some nodes the relative errors exceed 10%, including the lower bounds at nodes 9, 109 and 118, the upper bound at nodes 29 and 89, the relative errors of which are 14.57%, 11.08%, 14.16%, 12.47% and 10.41%, respectively. To improve the accuracy, the second order P-IFEM is also employed for displacement response bounds analysis, the results and corresponding relative errors compared to the MCS method are given in Tables 3 and 4. The relative errors are also provided in Figs. 19 and 20. From the results it can be found that the second order P-IFEM shows better accuracy than the first order P-IFEM for this case. The relative errors of the displacement bounds of du at marked nodes are all controlled within 6%, and those of dv are all less than 9%. However, compared to the first Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

26

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 15. Comparison of displacement (du ) bounds by the first order P-IFEM and the MCS method.

Fig. 16. Comparison of displacement (dv ) bounds by the first order P-IFEM and the MCS method.

order P-IFEM, the second order P-IFEM is more complicated in computational aspect, since that it requires more information of the system. 6.2. A plate with a hole A plate with a circular hole is shown in Fig. 21(a). The width and the length of the plate are 15 mm and 20 mm respectively, and the radius of the hole is 3.45 mm. A uniform tension p = 6 MPa is applied on the short edges of the plate. The Poisson’s ratio is ν = 0.3. The finite element mesh shown in Fig. 21(b) includes a total of 424 triangular elements and 738 nodes. The Young’s modulus E(x) ∈ E I (x) of the plate is modeled as a two-dimensional homogeneous interval field with the midpoint function E m (x) = 210 GPa and radius function E r (x) = 0.1 E m (x). The correlation coefficient function R E I (x1 , x2 ) is: ( ) |y1 − y2 |2 |x1 − x2 |2 − (91) R E I (x1 , x2 ) = exp − l x2 l y2 Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

27

Fig. 17. Relative error of displacement bounds (du ) by the first order P-IFEM.

Fig. 18. Relative error of displacement bounds (dv ) by the first order P-IFEM.

Fig. 19. Relative error of displacement bounds (du ) by the second order P-IFEM.

in which the correlation lengths are respectively l x = 0.06 m and l y = 0.08 m. The truncated interval K–L expansion is employed for representation of the Young’s modulus. The Nystr¨om method [78] is utilized to solve the Fredholm integral equation for derivation of the eigenvalues and eigenfunctions in the interval K–L expansion. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

28

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Table 1 Displacement bounds (du ) at the marked nodes by the first order P-IFEM and the MCS method. Node no.

1

3

9

101

First order P-IFEM (10−4 m) MCS (10−4 m)

[0.000, 0.000] [0.000, 0.000]

[1.562, 2.291] [1.677, 2.279]

[1.992, 2.908] [2.136, 2.857]

[2.911, 4.273] [3.125, 4.185]

Relative error

LB –

LB 6.85%

LB 6.77%

LB 6.85%

Node no.

21

23

29

103

First order P-IFEM (10−4 m) MCS (10−4 m)

[0.395, 0.604] [0.429, 0.596]

[0.874, 1.283] [0.942, 1.273]

[1.238, 1.814] [1.325, 1.778]

[1.516, 2.233] [1.627, 2.189]

Relative error

LB 7.90%

LB 7.24%

LB 6.61%

LB 6.81%

Node no. First order P-IFEM MCS (10−4 m)

(10−4

m)

Relative error Node no. First order P-IFEM MCS (10−4 m)

(10−4

m)

Relative error

UB –

UB 1.47%

UB 0.55%

UB 0.82%

UB 1.79%

UB 2.01%

UB 2.11%

UB 2.03%

81

83

89

109

[−0.604, − 0.395] [−0.598, − 0.430]

[−1.283, − 0.874] [−1.258, − 0.942]

[−1.814, − 1.238] [−1.779, − 1.320]

[−2.233, − 1.516] [−2.193, − 1.617]

LB 1.00%

LB 2.02%

LB 1.98%

LB 1.83%

UB 8.08%

UB 7.24%

UB 6.26%

UB 6.26%

121

112

118

120

[0.000, 0.000] [0.000, 0.000]

[−2.291, − 1.562] [−2.257, − 1.677]

[−2.908, − 1.992] [−2.848, − 2.120]

[−4.273, − 2.911] [−4.217, − 3.105]

LB –

LB 1.49%

LB 2.13%

LB 1.34%

UB –

UB 6.86%

UB 6.07%

UB 6.25%

UB: upper bound, LB: lower bound.

Table 2 Displacement bounds (dv ) at the marked nodes by the first order P-IFEM and the MCS method. Node no.

1

3

9

101

First order P-IFEM (10−4 m) MCS (10−4 m)

[0.000, 0.000] [0.000, 0.000]

[−0.398, − 0.259] [−0.392, − 0.283]

[0.168, 0.372] [0.196, 0.364]

[0.681, 1.092] [0.745, 1.074]

Relative error

LB –

LB 1.49%

LB 14.57%

LB 8.49%

Node no. First order P-IFEM MCS (10−4 m)

(10−4

m)

Relative error Node no.

UB –

UB 8.44%

UB 2.39%

UB 1.68%

21

23

29

103

[−1.361, − 0.931] [−1.351, − 1.005]

[−1.107, − 0.757] [−1.098, − 0.818]

[−0.143, − 0.075] [−0.138, − 0.085]

[0.135, 0.251] [0.150, 0.245]

LB 0.73%

LB 0.90%

LB 3.59%

LB 9.98%

UB 7.29%

UB 7.46%

UB 12.47%

UB 2.23%

81

83

89

109

[−1.361, − 0.931] [−1.332, − 1.003]

[−1.107, − 0.757] [−1.087, − 0.816]

[−0.143, − 0.075] [−0.139, − 0.083]

[0.135, 0.251] [0.152, 0.253]

Relative error

LB 2.16%

LB 1.90%

LB 3.11%

LB 11.08%

Node no.

121

112

118

120

First order P-IFEM (10−4 m) MCS (10−4 m)

[0.000, 0.000] [0.000, 0.000]

[−0.398, − 0.259] [−0.395, − 0.281]

[0.168, 0.372] [0.195, 0.367]

[0.681, 1.092] [0.749, 1.077]

Relative error

LB –

LB 0.76%

LB 14.16%

LB 8.97%

First order P-IFEM MCS (10−4 m)

(10−4

m)

UB 7.13%

UB –

UB 7.28%

UB 7.89%

UB 10.41%

UB 1.66%

UB 1.11%

UB 1.39%

UB: upper bound, LB: lower bound.

The existence of the hole in the plate will cause stress concentration nearby, especially at the location around the element of A, as marked in Fig. 21(b). To perform the stress analysis of the plate with a hole, firstly we focus Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

29

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 20. Relative error of displacement bounds (dv ) by the second order P-IFEM. Table 3 Displacement bounds (du ) at the marked nodes by the second order P-IFEM. Node no. Second order P-IFEM

1 (10−4

m)

3

9

101

[0.000, 0.000]

[1.631, 2.360]

[2.077, 2.994]

[3.040, 4.402]

Relative error

LB –

LB 2.74%

LB 2.76%

LB 2.73%

Node no.

21

23

29

103

Second order P-IFEM (10−4 m)

[0.416, 0.626]

[0.912, 1.322]

[1.292, 1.868]

[1.584, 2.302]

Relative error

LB 3.01%

LB 3.12%

LB 2.51%

LB 2.60%

Node no.

81

83

89

109

Second order P-IFEM (10−4 m)

[−0.625, − 0.416]

[−1.322, − 0.912]

[−1.868, − 1.292]

[−2.302, − 1.584]

Relative error

LB 4.52%

LB 5.10%

LB 5.04%

LB 4.95%

Node no. Second order P-IFEM

UB –

UB 5.01%

UB 3.19%

121 (10−4

m)

Relative error

UB 3.57%

UB 3.87%

UB 3.12%

112

UB 4.79%

UB 5.07%

UB 2.14%

118

UB 5.19%

UB 5.16%

UB 2.03%

120

[0.000, 0.000]

[−2.360, − 1.631]

[−2.994, − 2.077]

[−4.402, − 3.040]

LB –

LB 4.54%

LB 5.14%

LB 4.40%

UB –

UB 2.75%

UB 2.03%

UB 2.10%

UB: upper bound, LB: lower bound.

on evaluation of stress response bounds for the element of A, by means of the first order P-IFEM. With different number of terms M retained in the truncated interval K–L expansion, the response intervals of the normal stresses σx , σ y and the shear stress τx y of the element A are obtained and plotted in Fig. 22. It can be observed that the response bounds of all the normal stresses and the shear stress converge with M. When M reaches 19, the response bounds keep steady. Therefore, it can be recognized that the most principal 19 terms in the interval K–L expansion of the interval field E I (x) are necessary in this case for providing a reliable evaluation of the response bounds. When M = 19, the stress response intervals σxI of the normal stress σx of the whole plate are presented in Fig. 23. The response interval σxI at an arbitrary location in the plate are represented by a solid line. From Fig. 23 we can see, the stress response intervals are relatively large at locations around the hole, while those at other locations are much smaller. The response bounds nearby the hole are replotted in Fig. 24. The result indicates that the response interval of the normal stress σx at location of (x, y) = (0.00474, 0.0154) mm is σxI = [−21.75, −17.88] MPa, in which the lower bound 21.745 MPa is the maximum value of all possible normal stress responses over the whole plate. The response bounds of the normal stress σ y are also provided. In this case, the plate is stretched along the coordinate of y, therefore the maximum principal stress is the normal stress σ y . By the first order P-IFEM, the Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

30

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Table 4 Displacement bounds (dv ) at the marked nodes by the second order P-IFEM. Node no.

1

3

9

101

Second order P-IFEM (10−4 m)

[0.000, 0.000]

[−0.412, − 0.273]

[0.179, 0.391]

[0.722, 1.134]

Relative error

LB –

LB 5.10%

LB 8.51%

LB 3.04%

Node no. Second order P-IFEM

21 (10−4

m)

Relative error Node no. Second order P-IFEM

UB –

23

m)

UB 7.47%

29

UB 5.60%

103

[−1.401, − 0.972]

[−1.140, − 0.790]

[−0.150, − 0.080]

[0.145, 0.262]

LB 3.71%

LB 3.90%

LB 8.23%

LB 3.48%

UB 3.28%

81 (10−4

UB 3.46%

UB 3.44%

83

UB 6.91%

89

UB 7.01%

109

[−1.401, − 0.972]

[−1.140, − 0.790]

[−0.150, − 0.080]

[0.145, 0.262]

Relative error

LB 5.18%

LB 4.94%

LB 7.73%

LB 4.66%

Node no.

121

112

118

120

Second order P-IFEM (10−4 m)

[0.000, 0.000]

[−0.412, − 0.273]

[0.179, 0.391]

[0.722, 1.134]

Relative error

LB –

LB 4.34%

LB 8.08%

LB 3.55%

UB 3.12%

UB –

UB 3.24%

UB 2.88%

UB 4.71%

UB 6.69%

UB 3.52%

UB 5.29%

UB: upper bound, LB: lower bound.

Fig. 21. The plate with a hole under tension.

response bounds of σ y of the plate are given in Fig. 25. Due to stress concentration, Fig. 25 shows that the response intervals of σ y around the hole are obviously larger than those at other locations. The response bounds of σ y nearby the hole are replotted in Fig. 26, from which we know that the maximum stress σ y happens at location of (x, y) = (0.00143, 0.0121) mm. The response interval of σ y at this location is σ yI = [44.28, 53.78] MPa, which indicates the possible variation range of the normal stress σ y under the spatial uncertainty of the Young’s modulus. Indeed, this location is exactly the location of the element A, which from both stress concentration principle and engineering experience is recognized as the most dangerous position that may cause crack and failure of the structure. Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

31

Fig. 22. Effects of M on response intervals of stresses of element A.

Fig. 23. Response intervals of the normal stress σx .

7. Conclusions In this paper, an interval field model for quantification of spatially uncertain parameters is proposed. The fundamental framework and the basic characteristics of the interval field model are given. Specifically, an interval K–L expansion is presented as an effective representation of the interval field, thus greatly facilitating the application of the interval field model in subsequent uncertainty analysis. By retaining the most prominent terms in this expansion, a truncated interval K–L expansion and the error analysis are presented, thus realizing the quantification of the interval field considerably through very limited interval variables. By introducing the proposed interval field model to the finite element analysis of structures with spatially uncertain parameters, the non-deterministic equilibrium equations with interval factors are then formulated. The first order and the second order P-IFEMs are then proposed for this type of problems, in which the stiffness matrix is non-deterministic due to the spatially uncertain parameters. The feasibility and validity of the proposed interval field model and the P-IFEMs are verified. Compared with the first order P-IFEM, the second order P-IFEM can Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

32

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 24. Response intervals of the normal stress σx nearby the hole.

Fig. 25. Response intervals of the normal stress σ y .

provide more accurate solutions. However, the second order P-IFEM is relatively complex to perform, and it requires much more information of the system. The interval field is proposed as a supplement rather than a substitution for random field. It is applicable to problems when the sample information is insufficient for constructing a random field while the bounds of uncertainty Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

33

Fig. 26. Response intervals of the normal stress σ y nearby the hole.

can be given, especially in early design stages when objective probabilistic information often is not available. As a general spatial uncertainty quantification model, the interval field can also be combined with other numerical methods such as the finite difference method, the boundary finite element method, etc., thus formulating effective non-deterministic numerical methods for varieties of engineering problems. Acknowledgments This work is supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51621004), the National Natural Science Foundation of China (Grant No. 11802089), the National Science Fund for Distinguished Young Scholars (Grant No. 51725502), the Open Project Program of Key Laboratory for Precision & Non-traditional Machining of Ministry of Education, Dalian University of Technology (Grant No. JMTZ201701). Appendix A Proof. Denote HiKL (u) as the formulation on the right side of Eq. (27), namely: HiKL (u) = H m (u) +

∞ ∑

√ H r (u) λ j ϕ j (u)ζ j

(A.1)

j=1 I Then it is to prove that the variation range HiKL (u), the covariance function C H I (u1 , u2 ), and the correlation coeffiiKL { } cient function R H I (u1 , u2 ) of HiKL (u) are exactly accordant to those of the interval field H (u) ∈ H I (u), u ∈ D . iKL ∑ (1) Firstly, prove that C H I (u1 , u2 ) = C H I (u1 , u2 ) = H r (u1 )H r (u2 ) ∞ j=1 λ j ϕ j (u1 )ϕ j (u2 ). iKL ˆ The proof contains two parts, including (a) Derive the covariance C I (u1 , u2 ) of Hˆ iKL (u1 ) and Hˆ iKL (u2 ), where HiKL

Hˆ iKL (u1 ) and Hˆ iKL (u2 ) are respectively approximations of HiKL (u1 ) and HiKL (u2 ) by truncating the infinite series in Eq. (A.1) with the first M terms; and (b) prove that Cˆ H I (u1 , u2 ) converges to C H I (u1 , u2 ) as M → ∞. iKL

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

34

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

(a) Derive the covariance Cˆ H I (u1 , u2 ) of Hˆ iKL (u1 ) and Hˆ iKL (u2 ) which are approximations of HiKL (u1 ) and iKL HiKL (u2 ), respectively. From Eq. (A.1) we have: HiKL (u1 ) = H m (u1 ) +

∞ ∑

√ H r (u1 ) λ j ϕ j (u1 )ζ j

j=1 m

HiKL (u2 ) = H (u2 ) +

∞ ∑

(A.2) √ H (u2 ) λ j ϕ j (u2 )ζ j r

j=1

For arbitrary M > 0, M ∈ N = {1, 2, . . .}, the approximations of HiKL (u1 ) and HiKL (u2 ) can be expressed as: Hˆ iKL (u1 ) = H m (u1 ) +

M ∑

√ H r (u1 ) λ j ϕ j (u1 )ζ j

j=1 m

Hˆ iKL (u2 ) = H (u2 ) +

M ∑

(A.3) √ H (u2 ) λ j ϕ j (u2 )ζ j r

j=1

which can also be written as: ⎤ ζ1 ]⎢ζ ⎥ √ 2⎥ √λ M ϕ M (u1 ) ⎢ ⎢ ⎥ (A.4) λ M ϕ M (u2 ) ⎣ ... ⎦ ⎡

[

] [ m ] [ r Hˆ iKL (u1 ) H (u1 ) H (u1 ) = + H m (u2 ) Hˆ iKL (u2 )

] [√ √ λ1 ϕ1 (u1 ) √λ2 ϕ2 (u1 ) · · · √ H r (u2 ) λ1 ϕ1 (u2 ) λ2 ϕ2 (u2 ) · · ·

ζM

Denote: [ ] [ ] [ m ] hˆ iKL (u1 ) Hˆ iKL (u1 ) H (u1 ) = − H m (u2 ) hˆ iKL (u2 ) Hˆ iKL (u2 )

(A.5)

Eq. (A.4) equals: ⎤ ζ1 ]⎢ζ ⎥ √ 2⎥ √λ M ϕ M (u1 ) ⎢ ⎢ .. ⎥ λ M ϕ M (u2 ) ⎣ . ⎦ ⎡

[

] [ r hˆ iKL (u1 ) H (u1 ) = hˆ iKL (u2 )

] [√ √λ1 ϕ1 (u1 ) H r (u2 ) λ1 ϕ1 (u2 )

√ √λ2 ϕ2 (u1 ) λ2 ϕ2 (u2 )

··· ···

(A.6)

ζM

or in brief: hˆ = Tζ ζ

(A.7)

in which: [ ]T hˆ = hˆ iKL (u1 ) hˆ iKL (u2 ) ] [ r ] [√ √ √ H (u1 ) λ1 ϕ1 (u1 ) λ2 ϕ2 (u1 ) · · · λ M ϕ M (u1 ) Tζ = √ √ √ H r (u2 ) λ1 ϕ1 (u2 ) λ2 ϕ2 (u2 ) · · · λ M ϕ M (u2 ) [ ]T ζ = ζ1 ζ2 · · · ζ M

(A.8)

Assume that there is: ζ = Th hˆ

(A.9)

then we have: Th Tζ = I in which I is an identity matrix. Because there exists ζT ζ ≤ 1, which by Eq. (A.9) yields: hˆ T TTh Th hˆ ≤ 1

(A.10) ∑∞

j=1

ζ j2 ≤ 1, it can also be derived that

∑M

j=1

ζ j2 ≤ 1, namely (A.11)

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

35

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

With Eq. (A.10), it can be deduced that: TTh Th = (Tζ TTζ )−1

(A.12)

with which Eq. (A.11) can be rewritten as: hˆ T (Tζ TTζ )−1 hˆ ≤ 1

(A.13)

By comparing Eq. (A.13) with Eq. (6) that determines the uncertainty domain of correlated interval variables as shown in Fig. 2, the covariance Cˆ H I (u1 , u2 ) of interval variables Hˆ iKL (u1 ) and Hˆ iKL (u2 ) can be obtained as: iKL

Cˆ H I (u1 , u2 ) = H r (u1 )H r (u2 )

M ∑

iKL

λ j ϕ j (u1 )ϕ j (u2 )

(A.14)

j=1

(b) Prove that Cˆ H I (u1 , u2 ) converges to C H I (u1 , u2 ) as M → ∞. iKL ∑ Denote the sequence of Cn = H r (u1 )H r (u2 ) nj=1 λ j ϕ j (u1 )ϕ j (u2 ) as {Cn }. The element Cn is called the nth term of the sequence. Therefore, it is to prove that the limit of the sequence {Cn } is C H I (u1 , u2 ), namely, limn→∞ Cn = C H I (u1 , u2 ). ⏐∑ ⏐ ⏐ ⏐ From Eq. (29) we know that for every ελ > 0, there exists an index Nλ such that ⏐ ∞ λ j=n+1 j ⏐ < ελ for all n > Nλ . Denote Φ¯ j = maxu∈D ϕ 2j (u), then for ελ = H r (u )Hεr (u )·Φ¯ in which ε > 0 is arbitrary, there exists N such j 1 2 ⏐ ⏐∑ ⏐ ⏐ that ⏐ ∞ j=n+1 λ j ⏐ < ελ for all n > N . Subsequently, there is: ⏐ ⏐ ⏐Cn − C H I (u1 , u2 )⏐ ⏐ ⏐ ⏐∑ ⏐ ∞ ⏐ ⏐ r r ⏐ = H (u1 )H (u2 ) · ⏐ λ j ϕ j (u1 )ϕ j (u2 )⏐⏐ ⏐ j=n+1 ⏐ ⏐ ⏐ ⏐∑ ⏐ ⏐ ∞ ⏐ ≤ H r (u1 )H r (u2 ) · ⏐⏐ λ j ⏐⏐ · max ϕ 2j (u) ⏐ j=n+1 ⏐ u∈D

(A.15)

< ελ · H r (u1 )H r (u2 ) · Φ¯ j =ε ⏐ ⏐ Namely, for arbitrary ε > 0 there exists N such that ⏐Cn − C H I (u1 , u2 )⏐ < ε for all n > N . In other words, limn→∞ Cn = C H I (u1 , u2 ), which is equivalent to: lim Cˆ H I (u1 , u2 ) = C H I (u1 , u2 ) = H r (u1 )H r (u2 )

M→∞

∞ ∑

iKL

λ j ϕ j (u1 )ϕ j (u2 ) □

(A.16)

j=1

I (2) Secondly, prove that HiKL (u) ∈ HiKL (u) = H I (u) = [H m (u) − H r (u), H m (u) + H r (u)]. Denote: ∞ ∑ √ ∆HiKL (u) = HiKL (u) − H m (u) = H r (u) λ j ϕ j (u)ζ j

(A.17)

j=1

in which ζ j ∈ ζ I = [−1, 1], j = 1, 2, . . . are standard uncorrelated interval variables and satisfy √ ∑ r Because that for arbitrary ∆HiKL (u) = ∞ j=1 H (u) λ j ϕ j (u)ζ j , there exists: ∗ ∆HiKL (u)

=

∞ ∑

√ H r (u) λ j ϕ j (u)(−ζ j ) = −∆H (u)

∑∞

j=1

ζ j2 ≤ 1.

(A.18)

j=1

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

36

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

the midpoint function of ∆HiKL (u) is constantly equal to 0. Therefore, the midpoint function of the interval field H (u) represented by the interval K–L expansion in Eq. (27) equals H m (u), namely: m HiKL (u) = H m (u)

(A.19)

Based on Eq. (8), the radius function can be obtained as:   ∞ ∑ √ √  r HiKL (u) = C H I (u, u) = √(H r (u))2 λ j ϕ 2j (u) = H r (u) R H I (u, u) = H r (u)

(A.20)

j=1

[ ] I Therefore, there is HiKL (u) ∈ HiKL (u) = H I (u) = H m (u) − H r (u), H m (u)+Hr (u) if H (u) is represented by the interval K–L expansion in Eq. (27). □ Appendix B Proof. From the second part of Appendix A, it is not difficult to have H˜ m (u) = H m (u), namely, εm (u) = H m (u) − H˜ m (u) ≡ 0. Furthermore, we can also have: (

H˜ r (u)

)2

M ( )2 ∑ = H r (u) λ j ϕ 2j (u)

(B.1)

j=1

which yields: )2 ( (H r (u))2 − H˜ r (u) =1−

(H r (u))2

M ∑

λ j ϕ 2j (u)

(B.2)

j=1

Therefore there is: εr 2 (u) = =

=

=

)2 ( ∫ (H r (u))2 − H˜ r (u) 1 du |D| D (H r (u))2 ∫ M ∑ 1 λ j ϕ 2j (u))du (1 − |D| D j=1 ∫ ∫ ∑ M 1 1 1du − λ j ϕ 2j (u)du |D| D |D| D j=1 ∫ M ∑ 1 1− λj ϕ 2j (u)du |D| D

(B.3)

j=1 M

1 ∑ = 1− λj |D|

□

j=1

Appendix C Proof. By Eqs. (27), (32), (36) and (37) there are: h(u) =

∞ ∑ √

λ j ϕ j (u)ζ j

(C.1)

λ j ϕ j (u)ζ j

(C.2)

j=1

and ˜ h(u) =

M ∑ √ j=1

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

37

which yields: ∞ ∑ √

˜ h(u) − h(u) =

λ j ϕ j (u)ζ j

(C.3)

j=M+1

˜ Therefore the spatial error variance ε 2h (u) of h(u) is: ∫ ( ) 2 1 ˜ h(u) − h(u) du ε 2h (u) = |D| D ⎛ ⎞2 ∫ ∞ ∑ √ 1 ⎝ = λ j ϕ j (u)ζ j ⎠ du |D| D j=M+1 ∫ ∑ ∫ ∑ ∞ ∞ √ 1 1 = λ j ϕ 2j (u)ζ j2 du + λ j λk ϕ j (u)ϕk (u)ζ j ζk du |D| D |D| D j=M+1

(C.4)

j=M+1

j̸=k

1 = |D|

∞ ∑

λ j ζ j2

j=M+1

∫

∫ ∞ 1 ∑ √ 2 ϕ j (u)du + λ j λk ζ j ζk ϕ j (u)ϕk (u)du |D| j=M+1 D D j̸=k

With

ϕi (u)ϕ j (u)du = δi j we have: ∫ ∞ 1 ∑ 2 2 λjζj ε h (u) = ϕ 2j (u)du |D| D

∫

D

j=M+1

∞ 1 ∑ λ j ζ j2 = |D| j=M+1

∞ 1 ∑ = λj |D| j=M+1 ⎛ ⎞ M ∑ 1 ⎝ |D| − λj⎠ = |D|

(C.5)

j=1

Namely, ε 2h (u) ≤ 1 −

1 |D|

∑M

j=1

λj. □

Appendix D M T {The⏐∑sampling of } the vector ζ = [ζ1 ζ2 · · · ζ M ] from the M-dimensional unit sphere Ωunit ⏐ M 2 = ζ ⏐ j=1 ζ j ≤ 1 can be achieved by the following three steps:

Step 1: Sampling of the unit vector ζu . Firstly, generate M values c j , j = 1, 2, . . . , M by following the uniform distribution within [−a, a], where a is an arbitrary positive real value. Secondly, formulate the unit vector of ζu as ζu = ∥c∥c , where c = [c1 c2 · · · c M ]T and ∥c∥2 is the Euclidean norm of the vector c. 2

Step 2: Sampling of the radius value r, r ∈ [0, 1]. Generate the value of r by following the probability distribution of F(r ) = r M . Step 3: Obtain the sample of the vector ζ, namely, ζ = r · ζu . Note that the samples of ζ generated by the above specific steps will distribute uniformly within the MM dimensional unit sphere Ωunit . Indeed one can also use any other distribution forms in sampling of the unit vector ζu or the radius value r to obtain the samples of ζ. The only requirement is to ensure that any point in this sphere M Ωunit is possible to be obtained. References [1] E. Vanmarcke, Random Fields: Analysis and Synthesis, The MIT Press, Cambridge, 1983. [2] W.K. Liu, T. Belytschko, A. Mani, Random field finite elements, Internat. J. Numer. Methods Engrg. 23 (1986) 1831–1845.

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

38 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx R. Adler, The Geometry of Random Fields, Wiley, Chichester, 1981. S.D. Pietra, V.D. Pietra, J.D. Lafferty, Inducing features of random fields, IEEE Trans. Pattern Anal. Mach. Intell. 19 (1997) 380–393. R.G. Ghanem, P.D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer, New York, 1991. E. Vanmarcke, M. Shinozuka, S. Nakagiri, G.I. Schuëller, M. Grigoriu, Random fields and stochastic finite elements, Struct. Saf. 3 (1986) 143–166. A. Der Kiureghian, J.-B. Ke, The stochastic finite element method in structural reliability, Probab. Eng. Mech. 3 (1988) 83–91. G. Stefanou, The stochastic finite element method: Past, present and future, Comput. Methods Appl. Mech. Engrg. 198 (2009) 1031–1051. M. Kleiber, T.D. Hien, The Stochastic Finite Element Method: Basic Perturbation Technique and Computer Implementation, John Wiley & Sons, New York, 1992. A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc. 39 (1977) 1–38. G.A. Shafer, A Mathematical Theory of Evidence, Princeton University Press, Princeton, 1976. R.R. Yager, J. Kacprzyk, M. Fedrizzi, Advances in the Dempster-Shafer Theory of Evidence, John Wiley & Sons, New York, 1994. L.A. Zadeh, Fuzzy sets, Inf. Control 8 (1965) 338–353. L.A. Zadeh, Fuzzy sets as a basis for a theory of possibility, Fuzzy Sets Syst. 1 (1978) 3–28. R.E. Moore, Interval Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1966. E.R. Hansen, Topics in Interval Analysis, Oxford University Press, 1969. S. McWilliam, Anti-optimisation of uncertain structures using interval analysis, Comput. Struct. 79 (2001) 421–430. J.L. Wu, Z. Luo, Y.Q. Zhang, N. Zhang, L.P. Chen, Interval uncertain method for multibody mechanical systems using Chebyshev inclusion functions, Internat. J. Numer. Methods Engrg. 95 (2013) 608–630. Y. Ben-Haim, I. Elishakoff, Convex Models of Uncertainty in Applied Mechanics, Elsevier, Amsterdam, 1990. I. Elishakoff, P.G. Elisseeff, S.A.L. Glegg, Non-probabilistic, Convex-theoretic modeling of scatter in material properties, AIAA J. 32 (1994) 843–849. Z.P. Qiu, Convex models and interval analysis method to predict the effect of uncertain-but-bounded parameters on the buckling of composite structures, Comput. Methods Appl. Mech. Engrg. 194 (2005) 2175–2189. Z. Kang, Y.J. Luo, Non-probabilistic reliability-based topology optimization of geometrically nonlinear structures using convex models, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3228–3238. C. Jiang, X. Han, G.Y. Lu, J. Liu, Z. Zhang, Y.C. Bai, Correlation analysis of non-probabilistic convex model and corresponding structural reliability technique, Comput. Methods Appl. Mech. Engrg. 200 (2011) 2528–2546. H.-R. Bae, R.V. Grandhi, R.A. Canfield, An approximation approach for uncertainty quantification using evidence theory, Reliab. Eng. Syst. Saf. 86 (2004) 215–225. Z.P. Mourelatos, J. Zhou, A design optimization method using evidence theory, ASME J. Mech. Des. 128 (2005) 901–908. F. Massa, T. Tison, B. Lallemand, A fuzzy procedure for the static design of imprecise structures, Comput. Methods Appl. Mech. Engrg. 195 (2006) 925–941. Z. Zhang, X.X. Ruan, M.F. Duan, C. Jiang, An efficient epistemic uncertainty analysis method using evidence theory, Comput. Methods Appl. Mech. Engrg. 339 (2018) 443–466. C. Jiang, Q.F. Zhang, X. Han, J. Liu, D.A. Hu, Multidimensional parallelepiped model—a new type of non-probabilistic convex model for structural uncertainty analysis, Internat. J. Numer. Methods Engrg. 103 (2015) 31–59. R.L. Muhanna, H. Zhang, R.L. Mullen, Interval finite elements as a basis for generalized models of uncertainty in engineering mechanics, Reliab. Comput. 13 (2007) 173–194. H.U. Köylüo˘glu, A.S. ¸ Çakmak, S.R.K. Nielsen, Interval algebra to deal with pattern loading and structural uncertainties, J. Eng. Mech. 121 (1995) 1149–1157. Z.P. Qiu, S.H. Chen, H.B. Jia, The Rayleigh quotient iteration method for computing eigenvalue bounds of structures with bounded uncertain parameters, Comput. Struct. 55 (1995) 221–227. E.R. Hansen, On solving systems of equations using interval arithmetic, Math. Comp. 22 (1968) 374–384. G. Alefeld, J. Herzberger, Introduction to interval computations, SIAM Rev. 27 (1985) 296–297. A. Neumaier, Interval Methods for Systems of Equations, Cambridge University Press, New York, 1990. S.M. Rump, On the solution of interval linear systems, Computing 47 (1992) 337–353. R.L. Muhanna, R.L. Mullen, Formulation of fuzzy finite-element methods for solid mechanics problems, Comput.-Aided Civ. Infrastruct. Eng. 14 (1999) 107–117. D. Moens, D. Vandepitte, Fuzzy finite element method for frequency response function analysis of uncertain structures, AIAA J. 40 (2002) 126–136. M. Hanss, K. Willner, A fuzzy arithmetical approach to the solution of finite element problems with uncertain parameters, Mech. Res. Commun. 27 (2000) 257–272. S. Rao, J.P. Sawyer, Fuzzy finite element approach for the analysis of imprecisely-defined systems, in: 37th Structure, Structural Dynamics and Materials Conference, American Institute of Aeronautics and Astronautics, Salt Lake City, UT, 1996. R. Muhanna, R. Mullen, Uncertainty in mechanics problems—interval–based approach, J. Eng. Mech. 127 (2001) 557–566. M.V.R. Rao, R.L. Mullen, R.L. Muhanna, A new interval finite element formulation with the same accuracy in primary and derived variables, Int. J. Reliab. Saf. 5 (2011) 336–357. D. Moens, M. Hanss, Non-probabilistic finite element analysis for parametric uncertainty treatment in applied mechanics: recent advances, Finite Elem. Anal. Des. 47 (2011) 4–16.

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

39

[43] D. Moens, D. Vandepitte, A survey of non-probabilistic uncertainty treatment in finite element analysis, Comput. Methods Appl. Mech. Engrg. 194 (2005) 1527–1555. [44] S.S. Rao, L. Chen, Numerical solution of fuzzy linear equations in engineering analysis, Internat. J. Numer. Methods Engrg. 43 (1998) 391–408. [45] A. Pownuk, Numerical Solutions of Fuzzy Partial Differential Equations and Its Applications in Computational Mechanics, Springer, Berlin Heidelberg, 2004. [46] A. Neumaier, A. Pownuk, Linear systems with large uncertainties, with applications to truss structures, Reliab. Comput. 13 (2007) 149–172. [47] N. Xiao, R.L. Muhanna, F. Fedele, R.L. Mullen, Uncertainty analysis of static plane problems by intervals, Sae Int. J. Mater. Manuf. 8 (2015) 374–381. [48] N. Xiao, R.L. Muhanna, F. Fedele, R.L. Mullen, Interval finite element analysis of structural dynamic problems, Sae Int. J. Mater. Manuf. 8 (2015) 382–389. [49] Z.P. Qiu, I. Elishakoff, Antioptimization of structures with large uncertain-but-non-random parameters via interval analysis, Comput. Methods Appl. Mech. Engrg. 152 (1998) 361–372. [50] Z.P. Qiu, S.H. Chen, I. Elishakoff, Bounds of eigenvalues for structures with an interval description of uncertain-but-non-random parameters, Chaos Solitons Fractals 7 (1996) 425–434. [51] Z.P. Qiu, X.J. Wang, J.Y. Chen, Exact bounds for the static response set of structures with uncertain-but-bounded parameters, Int. J. Solids Struct. 43 (2006) 6574–6593. [52] S.H. Chen, X.W. Yang, Interval finite element method for beam structures, Finite Elem. Anal. Des. 34 (2000) 75–88. [53] D. Moens, M. De Munck, W. Desmet, D. Vandepitte, Numerical dynamic analysis of uncertain mechanical structures based on interval fields, in: Proc. IUTAM Symposium on the Vibration Analysis of Structures with Uncertainties, St. Petersburg, Russia, 2011, pp. 71-83. [54] W. Verhaeghe, W. Desmet, D. Vandepitte, D. Moens, Interval fields to represent uncertainty on the output side of a static FE analysis, Comput. Methods Appl. Mech. Engrg. 260 (2013) 50–62. [55] G. Muscolino, A. Sofi, M. Zingales, One-dimensional heterogeneous solids with uncertain elastic modulus in presence of long-range interactions: interval versus stochastic analysis, Comput. Struct. 122 (2013) 217–229. [56] A. Sofi, G. Muscolino, Static analysis of Euler–Bernoulli beams with interval Young’s modulus, Comput. Struct. 156 (2015) 72–82. [57] A. Sofi, G. Muscolino, I. Elishakoff, Static response bounds of Timoshenko beams with spatially varying interval uncertainties, Acta Mech. 226 (2015) 3737–3748. [58] A. Sofi, Structural response variability under spatially dependent uncertainty: Stochastic versus interval model, Probab. Eng. Mech. 42 (2015) 78–86. [59] D. Wu, W. Gao, Uncertain static plane stress analysis with interval fields, Internat. J. Numer. Methods Engrg. 110 (2017) 1272–1300. [60] L.P. Zhu, I. Elishakoff, J.H. Starnes, Derivation of multi-dimensional ellipsoidal convex model for experimental data, Math. Comput. Modelling 24 (1996) 103–114. [61] L. Fryba, N. Yoshikawa, Bounds analysis of a beam based on the convex model of uncertain foundation, J. Sound Vib. 212 (1998) 547–557. [62] X.J. Wang, I. Elishakoff, Z.P. Qiu, Experimental data have to decide which of the nonprobabilistic uncertainty descriptions—convex modeling or interval analysis—to utilize, J. Appl. Mech. 75 (2008) 041018. [63] Y.J. Luo, Z. Kang, Z. Luo, A. Li, Continuum topology optimization with non-probabilistic reliability constraints based on multi-ellipsoid convex model, Struct. Multidiscip. Optim. 39 (2009) 297–310. [64] I. Elishakoff, Y. Bekel, Application of Lamé’s super ellipsoids to model initial imperfections, J. Appl. Mech. 80 (2013) 061006–061009. [65] B.Y. Ni, C. Jiang, Z.L. Huang, Discussions on non-probabilistic convex modelling for uncertain problems, Appl. Math. Model. 59 (2018) 54–85. [66] C. Jiang, B.Y. Ni, X. Han, Y.R. Tao, Non-probabilistic convex model process: a new method of time-variant uncertainty analysis and its application to structural dynamic reliability problems, Comput. Methods Appl. Mech. Engrg. 268 (2014) 656–676. [67] C. Jiang, B.Y. Ni, N.Y. Liu, X. Han, J. Liu, Interval process model and non-random vibration analysis, J. Sound Vib. 373 (2016) 104–131. [68] C. Jiang, N.Y. Liu, B.Y. Ni, X. Han, Giving dynamic response bounds under uncertain excitations - a non-random vibration analysis method, Chin. J. Theor. Appl. Mech. 48 (2016) 447–463. [69] L. Wang, X.J. Wang, R.X. Wang, X. Chen, Time-dependent reliability modeling and analysis method for mechanics based on convex process, Math. Probl. Eng. 2015 (2015) 16. [70] L. Wang, X.J. Wang, X. Chen, R.X. Wang, Time-variant reliability model and its measure index of structures based on a non-probabilistic interval process, Acta Mech. 226 (2015) 3221–3241. [71] K. Karhunen, On Linear Methods in Probability Theory (English Translation By I. Selin), Annales Academiae Scientiarum Fennicae, Helsinki, 1947. [72] M. Loève, Probability Theory I, Springer-Verlag, New York, 1978. [73] R. Gutiérrez, J.C. Ruiz, M.J. Valderrama, On the numerical expansion of a second order stochastic process, Appl. Stoch. Models Data Anal. 8 (1992) 67–77. [74] K.K. Phoon, S.P. Huang, S.T. Quek, Simulation of second-order processes using Karhunen–Loeve expansion, Comput. Struct. 80 (2002) 1049–1060. [75] H.L. Van Trees, Detection, Estimation and Modulation Theory, Part 1, Wiley, New York, 1968. [76] W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Fredholm equations of the second kind, in: Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, New York, 2007.

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.

40

B.Y. Ni and C. Jiang / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

[77] K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. [78] W. Betz, I. Papaioannou, D. Straub, Numerical methods for the discretization of random fields by means of the Karhunen–Loève expansion, Comput. Methods Appl. Mech. Engrg. 271 (2014) 109–129. [79] R. Ghanem, The nonlinear Gaussian spectrum of log-normal stochastic processes and variables, J. Appl. Mech. 66 (1999) 964–973. [80] R.G. Ghanem, P.D. Spanos, Spectral stochastic finite-element formulation for reliability analysis, J. Eng. Mech. 117 (1991) 2351–2372. [81] J. Parvizian, A. Düster, E. Rank, Finite cell method, Comput. Mech. 41 (2007) 121–133. [82] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann, 2013. [83] T. Belytschko, W.K. Liu, B. Moran, K. Elkhodary, Nonlinear Finite Elements for Continua and Structures, Wiley, 2013. [84] D.L. Logan, A First Course in the Finite Element Method, Sixth ed., Cengage Learning, Boston, 2016. [85] T. Kato, Perturbation Theory for Linear Operators, second ed., Springer-Verlag Berlin Heidelberg, New York, 1976. [86] R. Bellman, Perturbation techniques in mathematics, physics, and engineering, Amer. Math. Monthly 72 (1965) 217. [87] L.W. Cohen, Lagrange multipliers for functions of infinitely many variables, Bull. Amer. Math. Soc. 40 (1934) 267–270. [88] M.S. Bazaraa, H.D. Sherali, C.M. Shetty, Nonlinear Programming: Theory and Algorithms, John Wiley & Sons, New York, 1993. [89] P.T. Boggs, J.W. Tolle, Sequential quadratic programming, Acta Numer. 4 (2008) 1–51.

Please cite this article as: B.Y. Ni and C. Jiang, Interval field model and interval finite element analysis, Computer Methods in Applied Mechanics and Engineering (2019) 112713, https://doi.org/10.1016/j.cma.2019.112713.