A substructure-based interval finite element method for problems with uncertain but bounded parameters

Abstract

In this article, the dependency between different elements in solid structures is considered and a substructure-based interval finite element method is used to model the interval properties. The penalty method is applied to impose the necessary constraints for compatibility. In order to obtain the interval stresses, an approximation solution based on the Taylor expansion method is presented. Then, the proposed interval substructure model is expanded to nonlinear problems. In consideration of the nonlinear property of the elasticity modulus, an interval elastoplastic substructure analysis method using constant matrix based on the incremental theory is proposed and the interval expression of the interval stress updated formation is derived. Finally, numerical examples are carried out to demonstrate the reasonability and feasibility of the proposed method and evaluation system.

Introduction

The concept of uncertainty plays an important role in the investigation of various sciences and engineering problems. ¹ Normally, uncertainties in reality are described as uncertain but bounded parameters. It is naturally compatible with the interval approach which only needs limited information about the bounds of variables.

The main shortage of interval approach to model uncertain problems is the dependency phenomenon which may result in some overestimation of the results.^2,3 In recent literature, plenty of effort has been made to overcome the drawback. Dong and Shah ⁴ and Rao and Berke ⁵ utilized the monotonic characters of variables to convert the interval problems to deterministic ones, and the bounds of responses are obtained based on the vertex combinations of the interval parameters. Rump ⁶ and Wang et al. ⁷ used the popular perturbation method to predict the response interval of the coupled structural–acoustic problems with uncertainties and the uncertain heat convection–diffusion problems. A collocation method was proposed further in the article by Wang et al. ⁸ to analyze the heat convection–diffusion problems with interval input parameters. The main point of these methods is to transfer the uncertain problems into determined ones. Unfortunately, the solutions may be trapped in local optimum.

In interval arithmetic, intervals are used directly and the results are guaranteed to include all the exact solutions. Rump ⁶ proposed an interval iteration method based on Brouwer’s fixed point theorem which has been shown to be efficient and accurate for solving linear interval equations. But it may result in meaninglessly wide and even catastrophic results without dependency consideration. An interval-based uncertainty treatment called element-by-element (EBE) technique was proposed by Muhanna and colleagues.^10–12 The global stiffness is an assembled EBE to eliminate the dependency in the element-coupling process. The EBE-based interval finite element method (EBE-IFEM) was also extended to calculate the envelope frequency response functions with uncertain parameters by Yang et al.,^13,14 and the results indicated the effectiveness of the EBE model to frame structures. Since frame structures are assembled by beams and trusses which are independent from each other, it is feasible to assume the independence of all the elements. Different from frame structures, parameters of solid structures are normally confirmed depending on the material partition. Thus, elements belonging to the same material are often considered to vary synchronously. In this situation, independent assumption of the elements is not appropriate. Besides, people seldom do research on nonlinear problems with uncertain parameters in interval arithmetic.

In this article, dependence between elements of solid structures is considered. A substructure-based interval finite element method (Sub-IFEM) is presented, and the interval solutions of element stresses are discussed. Then, researches on nonlinear problems are done based on the proposed interval substructure model (ISM). In section “Interval conceptions,” a brief summary of interval arithmetic and some basic theory are presented. Section “Sub-IFEM” introduces the proposed Sub-IFEM and the interval solutions of stresses. Section “Nonlinear ISM” shows the derivation process of the interval elastoplastic finite element method and the interval stress updated formulation. In section “Numerical examples,” numerical examples are applied to indicate the feasibility of the method. Finally, some conclusions are made in section “Conclusion.”

Interval conceptions

Basic definitions

In this article, matrix will be introduced in bold face and interval quantities (interval number, interval vector, and interval matrix) will be denoted by adding a superscript “I”

$$x^I=[\underset{¯}{x},\overline{x}]=\{x\in R:\underset{¯}{x}\le x\le \overline{x}\}$$

(1)

where $\underset{¯}{x}$ and $\overline{x}$ are, respectively, the lower and upper bounds of interval numbers. Let $x_c=(\overline{x}+\underset{¯}{x})/2$ and $x^{*}=(\overline{x}-\underset{¯}{x})/2$ be the middle point and the radius of $x^I$ , respectively, then $x^I$ can be expressed as

$$x^I=x_c+x^{*I}=[x_c-x^{*},x_c+x^{*}]$$

(2)

Sometimes, it is convenient to write the interval in the midpoint form

$$x^I=x_c(1+{\delta }^I)$$

(3)

in which ${\delta }^I$ is a symmetric interval number named the interval multiplier. For example, when we say x has 2% uncertainty level, it means ${\delta }^I=[-0.01,0.01]$ and $x^I=x_c(1+[-0.01,0.01])$ . ¹⁵

An interval vector consists of interval numbers. An interval matrix ${\mathbf{A}}^I$ is interpreted as a set of real $m\times n$ matrices by the convention

$${\mathbf{A}}^I=\{\mathbf{A}\in {\mathbf{R}}^{m\times n}|A_{\mathrm{ij}}\in A_{\mathrm{ij}}^I\mathrm{for}i=1,\dots ,m;j=1,\dots ,n\}$$

(4)

It can also be expressed as

$${\mathbf{A}}^I={\mathbf{A}}_c+{\mathbf{A}}^{*I}=[{\mathbf{A}}_c-{\mathbf{A}}^{*},{\mathbf{A}}_c+{\mathbf{A}}^{*}]={\mathbf{A}}_c(1+{\mathbf{d}}^I)$$

(5)

where ${\mathbf{A}}_c=(\underset{¯}{\mathbf{A}}+\overline{\mathbf{A}})/2$ and ${\mathbf{A}}^{*}=(\overline{\mathbf{A}}-\underset{¯}{\mathbf{A}})/2$ . ${\mathbf{d}}^I$ represents the interval multiplier matrix and ${\mathbf{d}}^I={\mathbf{A}}^{*I}/{\mathbf{A}}_c$ .

Dependency phenomenon in interval arithmetic

The interval arithmetic implicitly assumes that all intervals are independent even if they represent the same physical quantity. Therefore, $x^I-x^I$ is treated as $x^I-y^I$ and the result is not equal to zero. Actually, all the possible results of difference lying within $x^I$ are calculated instead of dependence between the two identical numbers. Ignorance of the dependency results in overestimation of the solutions.

Reducing occurrences of the same interval parameters in interval expressions can effectively eliminate such situation. Especially when the parameter appears only once, the interval formulation will give enough sharp enclosures of the results.

Penalty method

In steady-state analysis, the variational formulation of the model is given in the following form ¹⁶

$$\Pi =\frac{1}{2}{\mathbf{u}}^T\mathbf{Ku}-{\mathbf{u}}^T\mathbf{P}$$

(6)

with the conditions $\partial \Pi /\partial \mathbf{u}=0$

$$\partial \Pi /\partial \mathbf{u}=0$$

(7)

where $\Pi$ is the total potential energy. Assume that we want to impose the constraint conditions

$$\mathbf{Cu}-\mathbf{q}=0$$

(8)

where $\mathbf{C}$ and $\mathbf{q}$ consist of a series of constants. If $\mathbf{C}=[010-1]$ , it means $u_2-u_4=0$ , that is, $u_2=u_4$ . Define $\mathbf{t}=\mathbf{Cu}-\mathbf{q}$ and augment the potential energy function $\Pi$ by $(1/2){\mathbf{t}}^T\mathbf{\eta }\mathbf{t}$ , where $\mathbf{\eta }$ is a diagonal matrix of “penalty numbers” ${\eta }_i$ . Therefore

$$\tilde{\Pi }=\frac{1}{2}{\mathbf{u}}^T\mathbf{Ku}-{\mathbf{u}}^T\mathbf{P}+\frac{1}{2}{\mathbf{t}}^T\mathbf{\eta }\mathbf{t}$$

(9)

Invoking the stationary of $\tilde{\Pi }$ , that is, $\delta \tilde{\Pi }=0$ , results in

$$(\mathbf{K}+{\mathbf{C}}^T\mathbf{\eta }\mathbf{C})\mathbf{u}=\mathbf{p}+{\mathbf{C}}^T\mathbf{\eta }\mathbf{q}$$

(10)

or

$$(\mathbf{K}+\mathbf{Q})\mathbf{u}=\mathbf{p}+{\mathbf{C}}^T\mathbf{\eta }\mathbf{q}$$

(11)

Normally, the constraint conditions are $\mathbf{Cu}=0$ and $\mathbf{q}=0$ , then equation (11) can be reduced to

$$(\mathbf{K}+\mathbf{Q})\mathbf{u}=\mathbf{p}$$

(12)

where $\mathbf{Q}={\mathbf{C}}^T\mathbf{\eta }\mathbf{C}$ is called the penalty matrix. The element ${\eta }_i$ in $\mathbf{\eta }$ is often chosen from 10⁶ to 10⁸.

Sub-IFEM

ISM

In reality, Young’s modulus is usually uncertain and Poisson’s ratio is treated as a constant. ¹⁷ Assume that Young’s modulus is uncertain and described by an interval $E^I$ , where $E^I$ can be written as $E^I=E_c(1+{\delta }^I)$ .

Figure 1 shows a grid of a plane plate with four elements, which are of the same parameters. According to the article by Wang et al., ⁸ the interval global stiffness matrix of Figure 1 is given by the EBE model as follows

$${\mathbf{K}}^I={\mathbf{K}}_c(\mathbf{I}+{\mathbf{D}}^I),{\mathbf{D}}^I=\left(\begin{array}{cccc}{\mathbf{d}}_1^I\hfill & \hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\mathbf{d}}_2^I\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\mathbf{d}}_3^I\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & {\mathbf{d}}_4^I\hfill \end{array}\right)$$

(13)

where ${\mathbf{K}}^I$ is the interval global matrix, ${\mathbf{D}}^I$ and ${\mathbf{d}}^I$ are, respectively, the interval global multiplier matrix and the interval element multiplier matrix. ${\mathbf{d}}_i^I$ , $i=1,2,3,4$ , has the form of

$${\mathbf{d}}_i^I={\left(\begin{array}{ccc}{\delta }^I\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\delta }^I\hfill \end{array}\right)}_n$$

(14)

in which n stands for the degree of freedom of the elements.

Figure 1.

Grid of a plane plate with four elements.

${\mathbf{d}}_i^I$ in equation (13) is always treated as an independent one in the EBE model. Actually, the interval variables of the four elements in Figure 1 change synchronously, that is, they are dependent with each other. Thus, the use of the EBE model in this situation increases the occurrences of the same parameters, which will overestimate the interval results. Besides, the situation that the number of the same intervals increases with the number of the elements intensifies the expansion of the results. To overcome this drawback, the substructure method is used to decrease the multi-occurrence of the same intervals.

As the four elements are dependent with each other in Figure 1, they can be regarded as one substructure. Then, the global stiffness matrix can be assembled based on the local stiffness of the substructure. That is

$${\mathbf{K}}^I={\mathbf{K}}_c(\mathbf{I}+{\mathbf{D}}^I),{\mathbf{D}}^I={\mathbf{d}}^I$$

(15)

The matrix ${\mathbf{d}}^I$ is no longer the interval element multiplier matrix, and it represents the interval multiplier’s one of the substructure. Since there is only one substructure in Figure 1, ${\mathbf{D}}^I$ in equation (15) is just equal to ${\mathbf{d}}^I$ .

Generally, when invoking the dependency of variables, elements of the same parameters can be treated as a substructure, and the global matrix can be obtained as follows based on different substructures

$${\mathbf{K}}^I={\mathbf{K}}_c(\mathbf{I}+{\mathbf{D}}^I)$$

(16)

where ${\mathbf{K}}_c$ and ${\mathbf{D}}^I$ are written as

$${\mathbf{K}}_c=\left(\begin{array}{ccc}{\mathbf{K}}_{1_c}\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\mathbf{K}}_{N_c}\hfill \end{array}\right),{\mathbf{D}}^I=\left(\begin{array}{ccc}{\mathbf{d}}_1^I\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\mathbf{d}}_N^I\hfill \end{array}\right)$$

(17)

in which ${\mathbf{K}}_{\mathrm{ic}}$ , ${\mathbf{d}}_i^I$ , $i=1,\dots ,N$ , are, respectively, the mean matrix and the interval multiplier matrix of different substructures. N stands for the number of substructures. This substructure-based interval method is referred to as the ISM.

Governing equations for ISM and its solutions

According to equation (12), the governing equations of the ISM should be

$$({\mathbf{K}}^I+\mathbf{Q}){\mathbf{u}}^I={\mathbf{P}}^I$$

(18)

where ${\mathbf{u}}^I$ and ${\mathbf{P}}^I$ are, respectively

$${\mathbf{u}}^I=(\begin{array}{ccc}{\mathbf{u}}_1^I\hfill & \cdots \hfill & {\mathbf{u}}_N^I\hfill \end{array}{)}^T,{\mathbf{P}}^I=(\begin{array}{ccc}{\mathbf{P}}_1^I\hfill & \cdots \hfill & {\mathbf{P}}_N^I\hfill \end{array}{)}^T$$

(19)

where ${\mathbf{u}}_j^I$ and ${\mathbf{P}}_j^I$ , $j=1,\dots ,N$ , are, respectively, the displacement vector and the load vector of different substructures.

Let ${\mathbf{A}}^I=({\mathbf{K}}^I+\mathbf{Q})$ , using the interval iteration method proposed by Rump, ⁶ interval part of the solutions for equation (18) can be obtained as

$$({\mathbf{u}}^{*I}{)}^{n+1}=\mathbf{R}{\mathbf{P}}^I-\mathbf{R}{\mathbf{A}}^I{\mathbf{u}}_c+(\mathbf{I}-\mathbf{R}{\mathbf{A}}^I)({\epsilon }^I({\mathbf{u}}^{*I}{)}^n)$$

(20)

where $\mathbf{R}$ is a preconditioning matrix and ${\mathbf{A}}_c^{-1}$ is recommended to be a good choice of $\mathbf{R}$ . ¹⁸ ${\mathbf{u}}_c$ is calculated by $\mathbf{R}{\mathbf{P}}_c$ , and ${\mathbf{u}}^{*I}$ is the interval part of ${\mathbf{u}}^I$ . ${\epsilon }^I$ is an interval called the “inflation parameter” to accelerate the iteration process. ¹⁹ ${\epsilon }^I=[1-\beta ,1+\beta ]$ , where $\beta$ is a number between 0.001 and 0.1. The iterations stop when $({\mathbf{u}}^{*I}{)}^{n+1}\subseteq \mathrm{int}({\mathbf{u}}^{*I}{)}^n$ , and the algorithm converges if and only if $\rho (|(\mathbf{I}-\mathbf{R}{\mathbf{A}}^I)|)<1$ , that is, the spectral radius of the absolute value of $(\mathbf{I}-\mathbf{R}{\mathbf{A}}^I)$ is less than 1. The smaller the $\rho$ , the faster the convergence and the sharper the solutions.

Substituting ${\mathbf{A}}^I$ for ${\mathbf{A}}_c+{\mathbf{K}}_c{\mathbf{D}}^I$ and ${\mathbf{P}}^I$ for ${\mathbf{P}}_c+{\mathbf{P}}^{*I}$ in equation (20), the result is

$$\begin{array}{c}\multicolumn{1}{c}{({\mathbf{u}}^{*I}{)}^{n+1}=\mathbf{R}({\mathbf{P}}_c+{\mathbf{P}}^{*I})-\mathbf{R}({\mathbf{A}}_c+{\mathbf{K}}_c{\mathbf{D}}^I){\mathbf{u}}_c}\\ +(\mathbf{I}-\mathbf{R}({\mathbf{A}}_c+{\mathbf{K}}_c{\mathbf{D}}^I))({\epsilon }^I{\mathbf{u}}^{*I}{)}^n\hfill \\ \multicolumn{1}{c}{=\mathbf{R}{\mathbf{P}}_c+\mathbf{R}{\mathbf{P}}^{*I}-{\mathbf{u}}_c-\mathbf{R}{\mathbf{K}}_c{\mathbf{D}}^I({\mathbf{u}}_c+({\epsilon }^I({\mathbf{u}}^{*I}{)}^n)}\\ \multicolumn{1}{c}{=\mathbf{R}{\mathbf{P}}^{*I}-\mathbf{R}{\mathbf{K}}_c({\mathbf{M}}^I{)}^n{\mathbf{\Delta }}^I}\end{array}$$

(21)

The concept of the matrix multiplication is used here to represent ${\mathbf{D}}^I({\mathbf{u}}_c+({\epsilon }^I({\mathbf{u}}^{*I}{)}^n)$ as $({\mathbf{M}}^I{)}^n{\mathbf{\Delta }}^I$ , ¹⁵ where $({\mathbf{M}}^I{)}^n$ and ${\mathbf{\Delta }}^I$ in substructure model are given, respectively, by

$${(M^I\mathrm{)}}^l\mathrm{=}\left(\begin{array}{ccc}{u}_1+{\epsilon }^I{(u_1^{*I})}^l& & \\ & \ddots & \\ & & u_N+{\epsilon }^I{(u_N^{*I})}^l\end{array}\right),{\Delta }^I=\left(\begin{array}{c}{\mathrm{\delta }}_1^I\\ ⋮\\ {\mathrm{\delta }}_N^I\end{array}\right)$$

(22)

where ${\mathbf{u}}_j+{\epsilon }^I({\mathbf{u}}_j^{*I}{)}^l$ , for $j=1,\dots ,N$ , is the interval displacement vector for substructures. All deterministic calculations in equation (21) should be done before interval arithmetic to keep a sharp enclosure. If the applied load is determined, that is, ${\mathbf{P}}^I={\mathbf{P}}_c$ , equation (21) can be rewritten as

$$({\mathbf{u}}^{*I}{)}^{l+1}=-\mathbf{R}{\mathbf{K}}_c({\mathbf{M}}^I{)}^l{\mathbf{\Delta }}^I$$

(23)

Finally, interval solutions for the ISM are added up by

$${\mathbf{u}}^I={\mathbf{u}}_c+({\mathbf{u}}^{*I}{)}^{l+1}$$

(24)

Interval stresses

According to the traditional expression, interval stresses in each element can be obtained from

$${\mathbf{\sigma }}^I={\mathbf{D}}^I{\mathbf{Bu}}_e^I$$

(25)

where ${\mathbf{D}}^I$ is the interval element elasticity matrix, $\mathbf{B}$ is the strain matrix, and ${\mathbf{u}}_e^I$ is the element displacement vector calculated by ${\mathbf{u}}_e^I=\mathbf{L}{\mathbf{u}}^I$ , in which $\mathbf{L}$ is an element Boolean connectivity matrix. Rewriting ${\mathbf{D}}^I$ into ${\mathbf{D}}_c(1+{\delta }_i^I)$ and substituting equation (21) in equation (25), the interval stresses are obtained as

$${\mathbf{\sigma }}^I=(1+{\delta }_i^I)({\mathbf{\sigma }}_c+{\mathbf{D}}_c\mathbf{BLR}{\mathbf{P}}^{*I}-{\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c({\mathbf{M}}^I{)}^l{\mathbf{\Delta }}^I)$$

(26)

where ${\mathbf{\sigma }}_c$ is calculated by ${\mathbf{\sigma }}_c={\mathbf{D}}_c\mathbf{BL}{\mathbf{u}}_c$ and i represents the uncertain substructure to which the element belongs to. Due to the dependency between $(1+{\delta }_i^I)$ and $({\mathbf{M}}^I{)}^l{\mathbf{\Delta }}^I$ , the results of equation (26) get worse with respect to the increase in uncertainties. The Taylor expansion method is used to obtain reasonable interval results.

The first Taylor expansion of equation (26) at the midpoint of the interval parameters can be given by

$${\mathbf{\sigma }}^I={\mathbf{\sigma }}_c+\sum _{j=1}^N{{\delta }_j^I\frac{\partial {\mathbf{\sigma }}^I}{\partial {\delta }_j^I}|}_{{\mathbf{\Delta }}^I={\mathbf{\Delta }}_c}+R_2^I({\mathbf{\Delta }}^I)$$

(27)

where $R_2^I({\mathbf{\Delta }}^I)$ is the interval allowance and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{\sigma }}_c={\mathbf{D}}_c\mathbf{BL}{\mathbf{u}}_c}\\ \multicolumn{1}{c}{{\frac{\partial {\mathbf{\sigma }}^I}{\partial {\delta }_j^I}|}_{{\mathbf{\Delta }}^I={\mathbf{\Delta }}_c}=\{\begin{array}{c}({{\mathbf{\sigma }}_c+{\mathbf{\sigma }\prime }^I-{\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c{\mathbf{M}}^I{\mathbf{\Delta }}^I)|}_{{\mathbf{\Delta }}_c}+{(1+{\delta }_i^I)|}_{{\mathbf{\Delta }}_c}{(-{\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c{\mathbf{M}}^I\mathbf{\Delta }\prime )|}_{{\mathbf{\Delta }}_c}\hfill \\ \multicolumn{1}{c}{-{({\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c{\mathbf{M}}^I\mathbf{\Delta }\prime )|}_{{\Delta }_c}}\end{array}}\\ =\{\begin{array}{cc}{\mathbf{\sigma }}_c+{\mathbf{\sigma }\prime }^I-{\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c{\mathbf{M}}_c{\mathbf{\Delta }\prime }_i\hfill & j=i\hfill \\ \multicolumn{1}{c}{-{\mathbf{D}}_c\mathbf{BLR}{\mathbf{K}}_c{\mathbf{M}}_c{\mathbf{\Delta }\prime }_j}& j\ne i\hfill \end{array}\hfill \end{array}$$

(28)

in which ${\mathbf{\sigma }\prime }^I={\mathbf{D}}_c\mathbf{BLR}{\mathbf{P}}^{*I}$ , ${\mathbf{M}}_c=\mathrm{mid}({\mathbf{M}}^I)$ , and ${\mathbf{\Delta }\prime }_m=(0_1,\dots ,1_m,\dots ,0_N{)}^T$ , $m=i,j$ . When the parameters are at a low level of uncertainty, $R_2^I({\mathbf{\Delta }}^I)$ can be omitted and an approximation of the interval stresses is given by

$${\mathbf{\sigma }}^I={\mathbf{\sigma }}_c+\sum _{j=1}^N{{\delta }_j^I\frac{\partial {\mathbf{\sigma }}^I}{\partial {\delta }_j^I}|}_{{\mathbf{\Delta }}^I={\mathbf{\Delta }}_c}$$

(29)

Each interval only appears once in equation (29) and the dependency eliminates.

Nonlinear ISM

Considering the nonlinear material characteristics, the ISM can be extended to solve nonlinear interval problems. Then, the governing equations of the nonlinear ISM is

$$({\mathbf{K}}^I\left({\mathbf{u}}^I\right)+\mathbf{Q}){\mathbf{u}}^I={\mathbf{P}}^I$$

(30)

The main idea to solve equation (30) is to divide the nonlinear procedure into several linear steps. Then, the equation can be written in the incremental form as

$$({\mathbf{K}}^I\left(\Delta {\mathbf{u}}^I\right)+\mathbf{Q})\Delta {\mathbf{u}}^I=\Delta {\mathbf{P}}^I$$

(31)

where $\Delta$ is a sign of increment (the same below). The way to solve equation (31) is to update the interval global stiffness matrix ${\mathbf{K}}^I$ with $\Delta {\mathbf{u}}^I$ and then calculate the linear interval equations according to equation (24). However, the preconditioning matrix that is chosen as the inverse of the coefficient matrix in equation (31) should be given at the beginning of each linear step. This may utilize 78.4%–90.7% of the total CPU time. ¹⁵ The following parts propose an interval constant stiffness iteration method where the inverse matrix needs to be calculated only once during the complete computational process.

Interval constant stiffness iteration method for elastoplastic problems

The incremental form of the elastoplastic constitutive equation is

$$\Delta \mathbf{\sigma }=({\mathbf{D}}_e-(1-r){\mathbf{D}}_p)\Delta \mathbf{\epsilon },\{\begin{array}{cc}r=1\hfill & \mathrm{elastic}\hfill \\ \multicolumn{1}{c}{0<r<1}& \mathrm{transition}\mathrm{field}\hfill \\ \multicolumn{1}{c}{r=0}& \mathrm{plastic}\hfill \end{array}$$

(32)

where $r=-f_{m-1}/(f_m-f_{m-1})$ and f is the yield function. $f_m$ and $f_{m-1}$ , respectively, represent the value of f after and before loading. ${\mathbf{D}}_e$ and ${\mathbf{D}}_p$ are, respectively, the elastic and plastic matrices. Let $\Delta {\mathbf{\sigma }}_0=-(1-r){\mathbf{D}}_p\Delta \mathbf{\epsilon }$ , then equation (32) can be written as

$$\Delta \mathbf{\sigma }={\mathbf{D}}_e\Delta \mathbf{\epsilon }+\Delta {\mathbf{\sigma }}_0$$

(33)

here, $\Delta {\mathbf{\sigma }}_0$ is called the initial stress. According to the principle of virtual work, we have

$$\Delta {\mathbf{u}}_e^T\Delta {\mathbf{R}}_e={\int }_V\Delta {\mathbf{\epsilon }}_e^T\Delta \mathbf{\sigma }\mathrm{dV}$$

(34)

where $\Delta {\mathbf{u}}_e^T$ , $\Delta {\mathbf{\epsilon }}_e^T$ , and $\Delta {\mathbf{R}}_e$ are, respectively, the displacement increment, the stain increment, and the load increment of an element. Substituting $\Delta {\mathbf{\epsilon }}_e=\mathbf{B}\Delta {\mathbf{u}}_e$ and equation (33) in equation (34) results in

$${\int }_V{\mathbf{B}}^T{\mathbf{D}}_e\mathbf{B}\mathrm{dV}\Delta {\mathbf{u}}_e=\Delta {\mathbf{R}}_e-{\int }_V{\mathbf{B}}^T\Delta {\mathbf{\sigma }}_0\mathrm{dV}$$

(35)

Define ${\mathbf{K}}_e={\int }_V{\mathbf{B}}^T{\mathbf{D}}_e\mathbf{B}\mathrm{dV}$ as the element stiffness matrix and $\Delta {\mathbf{R}}_{e0}=-{\int }_V{\mathbf{B}}^T\Delta {\mathbf{\sigma }}_0\mathrm{dV}$ as the element pseudo load vector, then equation (35) changes into the form as

$${\mathbf{K}}_e\Delta {\mathbf{u}}_e=\Delta {\mathbf{R}}_e+\Delta {\mathbf{R}}_{e0}$$

(36)

For ISM, adding up all the elements, the constant stiffness governing equations can be written as

$$({\mathbf{K}}^I+\mathbf{Q})\Delta {\mathbf{u}}^I=\Delta {\mathbf{R}}^I+\Delta {\mathbf{R}}_0^I\mathrm{or}{A}^I\Delta {\mathbf{u}}^I=\Delta {\mathbf{R}}^I+\Delta {\mathbf{R}}_0^I$$

(37)

Then, the equations can be solved by updating the interval pseudo load vector $\Delta {\mathbf{R}}_0^I$ .

Interval stress updating method

The key step for calculating equation (37) is to update $\Delta {\mathbf{R}}_0^I$ . Since it is associated with $\Delta {\mathbf{\sigma }}_0^I$ , the result can be calculated by

$$\Delta {\mathbf{\sigma }}_0^I=-(1-r){\mathbf{D}}_p^I\Delta {\mathbf{\epsilon }}^I$$

(38)

According to chapter 3, we have $\Delta {\mathbf{\epsilon }}^I=\mathbf{BL}\Delta {\mathbf{u}}^I=\mathbf{BL}\Delta {\mathbf{u}}_c-\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}^{\mathrm{Il}}{\mathbf{\Delta }}^I$ , then equation (38) is

$$\Delta {\mathbf{\sigma }}_0^I=-(1-r){\mathbf{D}}_p^I(\mathbf{BL}\Delta {\mathbf{u}}_c-\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}^{\mathrm{Il}}{\mathbf{\Delta }}^I)$$

(39)

The Taylor expansion method is used again similar to equation (27), then we obtain

$$\Delta {\mathbf{\sigma }}_0^I=-(1-r)(\Delta {\mathbf{\sigma }}_{0c}+{\sum _{j=1}^N{\delta }_j^I\frac{\partial \Delta {\mathbf{\sigma }}_0^I}{\partial {\delta }_j^I}|}_{{\Delta }_c})$$

(40)

where ${\delta }_j^I$ is the interval multiplier of the jth substructure, and N stands for the number of substructures

$$\begin{array}{c}\multicolumn{1}{c}{\Delta {\mathbf{\sigma }}_{0c}={{{\mathbf{D}}_p^I|}_{{\mathbf{\Delta }}_c}(\mathbf{BL}\Delta {\mathbf{u}}_c-\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}^{\mathrm{Il}}{\mathbf{\Delta }}^I)|}_{{\mathbf{\Delta }}_c}={\mathbf{D}}_P\mathbf{BL}\Delta {\mathbf{u}}_c}\\ \multicolumn{1}{c}{{\frac{\partial \Delta {\mathbf{\sigma }}_0^I}{\partial {\delta }_j^I}|}_{{\mathbf{\Delta }}_c}={({\mathbf{D}}_p^I)\prime |}_{{\mathbf{\Delta }}_c}{(\mathbf{BL}\Delta {\mathbf{u}}_c-\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}^{\mathrm{Il}}{\mathbf{\Delta }}^I)|}_{{\mathbf{\Delta }}_c}}\\ +{{\mathbf{D}}_p^I|}_{{\mathbf{\Delta }}_c}{(-\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}^{\mathrm{Il}}\mathbf{\Delta }\prime )|}_{{\mathbf{\Delta }}_c}\hfill \\ ={({\mathbf{D}}_p^I)\prime |}_{{\mathbf{\Delta }}_c}\mathbf{BL}\Delta {\mathbf{u}}_c-{\mathbf{D}}_p\mathbf{BLR}{\mathbf{K}}_c\Delta {\mathbf{M}}_c\mathbf{\Delta }\prime \hfill \end{array}$$

(41)

The interval plastic matrix has the form

$${\mathbf{D}}_p^I=\frac{{\mathbf{D}}_e^I\left(\frac{\partial f}{\partial {\mathbf{\sigma }}^I}\right){\left(\frac{\partial f}{\partial {\mathbf{\sigma }}^I}\right)}^T{\mathbf{D}}_e^I}{H+{\left(\frac{\partial f}{\partial {\mathbf{\sigma }}^I}\right)}^T{\mathbf{D}}_e^I\left(\frac{\partial f}{\partial {\mathbf{\sigma }}^I}\right)}$$

(42)

where H is the hardening parameter. For idea elastoplastic problems, H is equal to 0. Defining ${\mathbf{a}}^I=\partial f/\partial {\mathbf{\sigma }}^I$ , $A=({\mathbf{a}}^I{)}^T{\mathbf{D}}_e^I{\mathbf{a}}^I$ , and writing ${\mathbf{D}}_e^I=(1+{\delta }_j^I){\mathbf{D}}_e$ , the plastic matrix changes to

$${\mathbf{D}}_p^I=\frac{1}{A}(1+{\delta }_j^I){\mathbf{D}}_e{\mathbf{a}}^I({\mathbf{a}}^I{)}^T{\mathbf{D}}_e$$

(43)

For ${\mathbf{D}}_p^I$ is a function of ${\mathbf{\sigma }}^I$ and ${\mathbf{\sigma }}^I$ contains the vector ${\mathbf{\Delta }}^I$ , when $j=i$ , ${\mathbf{D}}_p^I$ is not only an explicit function but also an implicit function of ${\delta }_j^I$ ; while $j\ne i$ , ${\mathbf{D}}_p^I$ is merely an explicit function of ${\delta }_j^I$ . Thus the differential of ${\mathbf{D}}_p^I$ is

$${({\mathbf{D}}_p^I)\prime |}_{{\mathbf{\Delta }}_c}=\{\begin{array}{cc}{\frac{\partial {\mathbf{D}}_p^I({\delta }_i^I,{({\mathbf{\sigma }}^I)}^{*})}{\partial {\delta }_i^I}|}_{{\Delta }_c}+{\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\delta }_i^I}|}_{{\Delta }_c}\hfill & j=i\hfill \\ \multicolumn{1}{c}{{\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\delta }_j^I}|}_{{\Delta }_c}}& j\ne i\hfill \end{array}$$

(44)

where the symbol “ $(•{)}^{*}$ ” indicates that the variable is constant to take a derivative.

So, when $j=i$ , the first item is

$${\frac{\partial {\mathbf{D}}_p^I({\delta }_i^I,{({\mathbf{\sigma }}^I)}^{*})}{\partial {\delta }_i^I}|}_{{\Delta }_c}={\frac{{\mathbf{D}}_e({\mathbf{a}}^I){({\mathbf{a}}^I)}^T{\mathbf{D}}_e}{{({\mathbf{a}}^I)}^T{\mathbf{D}}_e{\mathbf{a}}^I}|}_{{\Delta }_c}={\mathbf{D}}_p$$

(45)

The second is

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\delta }_i^I}={\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\mathbf{\sigma }}^I}|}_{{\Delta }_c}{\frac{\partial {\mathbf{\sigma }}^I}{\partial {\delta }_i^I}|}_{{\Delta }_c}}\\ \multicolumn{1}{c}{=\frac{\partial {\mathbf{D}}_p({\mathbf{\sigma }}^h)}{\partial {\mathbf{\sigma }}^h}{(\frac{\partial {({\mathbf{\sigma }}^I)}^{h-1}}{\partial {\delta }_i^I}+\frac{\partial {(\Delta {\mathbf{\sigma }}^I)}^{k-1}}{\partial {\delta }_i^I})|}_{{\Delta }_c}}\end{array}$$

(46)

when $j\ne i$ , the express is similar to equation (46), that is

$${\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\delta }_j^I}|}_{{\Delta }_c}=\frac{\partial {\mathbf{D}}_p({\mathbf{\sigma }}^h)}{\partial {\mathbf{\sigma }}^h}{(\frac{\partial {({\mathbf{\sigma }}^I)}^{h-1}}{\partial {\delta }_j^I}+\frac{\partial {(\Delta {\mathbf{\sigma }}^I)}^{k-1}}{\partial {\delta }_j^I})|}_{{\Delta }_c}$$

(47)

where ${\mathbf{\sigma }}^I=({\mathbf{\sigma }}^I{)}^{h-1}+(\Delta {\mathbf{\sigma }}^I{)}^{k-1}$ , here $({\mathbf{\sigma }}^I{)}^{h-1}$ is the total interval stress after the increment of $h-1$ and $(\Delta {\mathbf{\sigma }}^I{)}^{k-1}$ is the stress increment of the $(k-1)\mathrm{th}$ iteration in the $h\mathrm{th}$ increment step. They can be calculated before the current iteration step.

To solve equations (46) and (47) directly is very complicated, hence the difference method is suggested ²⁰ and the corresponding expression is

$$\frac{\partial {\mathbf{D}}_P({\mathbf{\sigma }}^h)}{\partial {\mathbf{\sigma }}^h}\approx \frac{{\mathbf{D}}_P({\mathbf{\sigma }}^h+\beta {\mathbf{\sigma }}^h)-{\mathbf{D}}_P({\mathbf{\sigma }}^h)}{\beta {\mathbf{\sigma }}^h}$$

(48)

Let $\mathbf{\alpha }={(\frac{\partial {({\mathbf{\sigma }}^I)}^{h-1}}{\partial {\delta }_j^I}+\frac{\partial {({\mathbf{\sigma }}^I)}^{k-1}}{\partial {\delta }_j^I})}_{{\mathbf{\Delta }}_c}$ , here $\mathbf{\alpha }$ is a column vector with four elements, then equations (46) and (47) can be rewritten as

$${\frac{\partial {\mathbf{D}}_p^I(({\delta }_i^I{)}^{*},{\mathbf{\sigma }}^I)}{\partial {\delta }_j^I}|}_{{\mathbf{\Delta }}_c}=\sum _{k=1}^4{\alpha }_k{\mathbf{D}}_{\mathrm{pk}}$$

(49)

where ${\alpha }_k$ is the element of $\alpha$ , and ${\mathbf{D}}_{\mathrm{pk}}$ is the kth difference of ${\mathbf{D}}_P({\mathbf{\sigma }}^h)$ . For plane strain problems, there are a total of four differences for four stress components, that is, $k=1,2,3,4$ . Normally, fine results can be obtained when $\beta$ varies from 0.001 to 0.01. ${\mathbf{D}}_p$ is given by the traditional method.

Finally, the updated stresses can be given by adding up the interval elastic stresses and the interval initial stresses

$$\Delta {\sigma }^I=\Delta {\sigma }_e^I+\Delta {\sigma }_0^I=\Delta {\sigma }_c+{\displaystyle \sum _{j=1}^N{{\delta }_j^I\left(\mathrm{(}r\mathrm{-1)}\frac{\partial \Delta {\sigma }_0^I}{\partial {\delta }_j^I}\right)|}_{{\Delta }_c}+{\frac{\partial \Delta {\sigma }_e^I}{\partial {\delta }_j^I}|}_{{\Delta }_c}}$$

(50)

The pseudo load vector is given by

$$\begin{array}{c}\Delta {\mathbf{R}}_0^I=\Delta {\mathbf{R}}_{0c}+\Delta {\mathbf{R}}_0^{*I}\hfill \\ \Delta {\mathbf{R}}_{0c}=\sum _{p=1}^{\mathrm{np}}\int (1-r){\mathbf{B}}^T\Delta {\mathbf{\sigma }}_{0c}\mathrm{dV},\hfill \\ \Delta {\mathbf{R}}_0^{*I}=\sum _{p=1}^{\mathrm{np}}\int (1-r){\mathbf{B}}^T{\sum _{j=1}^N{\delta }_j^I\frac{\partial \Delta {\mathbf{\sigma }}_0^I}{\partial {\delta }_j^I}|}_{{\Delta }_c}\mathrm{dV}\hfill \end{array}$$

(51)

Solutions of the elastoplastic equations

Substituting equation (51) in equation (37), the interval displacement of each linear step is

$$\begin{array}{c}\Delta {\mathbf{u}}_c=\mathbf{R}(\Delta {\mathbf{R}}_c+\Delta {\mathbf{R}}_{0c})\hfill \\ \Delta {\mathbf{u}}^{*I}=\mathbf{R}(\Delta {\mathbf{P}}^{*I}+\Delta {\mathbf{R}}_0^{*I})-\mathbf{R}{\mathbf{K}}_c{\mathbf{M}}^I{\mathbf{\Delta }}^I\hfill \end{array}$$

(52)

Update the pseudo load vector constantly until it converges.

For interval variables, it is hard to check the iteration method’s convergence. As we can see in equation (52), the interval part is accompanied by the deterministic results. Thus, the iteration criterion and the yield conditions of the interval iteration process can still follow the traditional way. Adding up all the incremental steps, we finally reach the results as

$${\mathbf{u}}^I=\sum _{h=1}^m{(\Delta {\mathbf{u}}_c+\Delta {\mathbf{u}}^{*I})}^h,{\mathbf{\sigma }}^I=\sum _{h=1}^m\Delta {\mathbf{\sigma }}^I$$

(53)

Numerical examples

Linear ISM

The calculation model is shown in Figure 2, where the upper solid is made of concrete with the material constants of $E_{1c}=2.4\times {10}^4\mathrm{MPa}$ , $\mu =0.167$ , and $\gamma =24\mathrm{KN}/{\mathrm{m}}^3$ , parameters of the foundation are $E_{2c}=2.0\times {10}^4\mathrm{MPa}$ and $\mu =0.25$ (regardless of the deadweight of the foundation) and parameters of the soft interlayer located between the upper part and the foundation are $E_{3c}=1.25\times {10}^3\mathrm{MPa}$ , $\mu =0.25$ , $\gamma =21.5\mathrm{KN}/{\mathrm{m}}^3$ , $c=0.1\mathrm{MPa}$ , and $\phi ={32}^{\circ }$ . The thickness of the soft interlayer is 0.02 m and the applied load F is 300 KN/m.

Figure 2.

Two-dimensional (2D) model with a soft interlayer.

As shown in Figure 3, there are 60 elements and 80 nodes altogether, where the soft interlayer is divided into four elements. The left side, the right side, and the bottom of the foundation are under normal constraints. Assuming elements of the same material to be dependent, the solid can be treated as a model with three substructures. Suppose the elasticity moduli of all the substructures are intervals with 8% uncertainties. The load varies in [290, 310] KN/m.

Figure 3.

Finite element analysis (FEA) mesh model of the structures.

To solve this system, the Sub-IFEM and the vertex approach are used. Figure 4 shows the displacement of the nodes on the interface between the upper part and the soft interlayer. The result shows that the Sub-IFEM provides a sharper exclusion of the vertex method compared to the EBE-IFEM. The normal stresses ${\sigma }_y$ and the shearing stresses ${\tau }_{\mathrm{xy}}$ of the interlayer elements are shown, respectively, in Tables 1 and 2.

Figure 4.

X-displacement of nodes on the interface between the upper part and the soft interlayer.

Table 1.

Normal stresses of the interlayer elements ${\sigma }_y$ (unit: 10 KN).

Element ID	EBE-IFEM	Sub-IFEM	Vertex method
57	[6.679, 9.279]	[7.038, 8.892]	[7.128, 8.788]
58	[−10.738, −9.363]	[−10.437, −9.639]	[−10.339, −9.741]
59	[−8.219, −7.170]	[−7.909, −7.357]	[−7.889, −7.485]
60	[−30.749, −26.773]	[−29.91, −27.556]	[−29.467, −27.980]

EBE-IFEM: element-by-element–based interval finite element method; Sub-IFEM: substructure-based interval finite element method.

Table 2.

Sharing stresses of the interlayer elements ${\tau }_{\mathrm{xy}}$ (unit: 10 KN).

Element ID	EBE-IFEM	Sub-IFEM	Vertex method
57	[22.448, 26.477]	[22.720, 26.175]	[23.340, 25.510]
58	[7.964, 9.488]	[8.150, 9.280]	[8.220, 9.205]
59	[8.233, 9.767]	[8.414, 9.564]	[8.547, 9.415]
60	[16.684, 19.473]	[16.926, 19.204]	[17.164, 18.957]

EBE-IFEM: element-by-element–based interval finite element method; Sub-IFEM: substructure-based interval finite element method.

In Tables 1 and 2, the stresses given by EBE-IFEM are calculated by equation (29) and the stresses given by Sub-IFEM are calculated by equation (32). The results show that the Sub-IFEM provides sharper enclosures of the vertex method compared to the EBE-IFEM. Besides, the number of degrees of freedom in the Sub-IFE model is 157 which reduce the calculation cost greatly compared to the value in the EBE-IFE model with 440.

Nonlinear ISM

A classical thick-walled cylinder is taken as an example. Considering its axisymmetric characteristic, a grid of quarter of the cylinder is shown in Figure 5. A gradually increasing pressure P with $P_{max}=200\mathrm{MPa}$ is applied inside the cylinder. Young’s modulus, Poisson’s ratio, and the yield stress of the cylinder are, respectively, E = 2.1 × 10⁵ MPa, µ = 0.3 and σ_s = 260 MPa. Supposing Young’s modulus varies in [1.89, 2.31] × 10⁵ MPa, two different conditions are discussed as follows using the proposed nonlinear ISM.

Figure 5.

Grid of quarter of a cylinder.

Elastic interval analysis

In order to check the applicability of the incremental theory to interval problems, an elastic procedure is first carried out using the proposed nonlinear interval method. The yield stress is set at a high value to avoid the plasticity. Then, the applied load is, respectively, divided into 1, 3, 5, 7, and 10 different increments. Table 3 gives the displacement of node A in both X and Y directions. Table 4 shows the stresses of the first integral point in element 6. The distribution diagram of the Gaussian integral points in element 6 is shown in Figure 6.

Table 3.

Interval displacement at point A (unit: m).

Load increments	X-displacement		Y-displacement
	Middle point	Interval radius	Middle point	Interval radius
1	1.27485E−4	1.41650E−5	1.27485E−4	1.41650E−5
3	1.27485E−4	1.41650E−5	1.27485E−4	1.41650E−5
5	1.27485E−4	1.41650E−5	1.27485E−4	1.41650E−5
7	1.27485E−4	1.41650E−5	1.27485E−4	1.41650E−5
10	1.27485E−4	1.41650E−5	1.27485E−4	1.41650E−5

Table 4.

Interval stresses of the first integral point in element 6 (unit: KPa).

Load increments	Horizontal stresses ${\sigma }_x$		Vertical stresses ${\sigma }_y$
	Middle point	Interval radius	Middle point	Interval radius
1	0.32746E + 05	0.68723E−06	−0.16444E + 05	0.65389E−07
3	0.32746E + 05	0.68723E−06	−0.16444E + 05	0.65389E−07
5	0.32746E + 05	0.68723E−06	−0.16444E + 05	0.65389E−07
7	0.32746E + 05	0.68723E−06	−0.16444E + 05	0.65389E−07
10	0.32746E + 05	0.68723E−06	−0.16444E + 05	0.65389E−07

Figure 6.

Gaussian integral points in element 6.

The results in Tables 3 and 4 show that the interval displacements and the interval stresses obtained based on the incremental theory are not affected by change in the load increments. It indicates that the computation of each linear increment is independent with others and the proposed way to replace the nonlinear interval problems with several linear interval procedures is feasible.

Besides, according to the analytical solutions of the thick-walled cylinder, the stresses can be given as follows

$${\sigma }_r=\frac{a^{2}P}{b^2-a^2}(1-\frac{b^2}{r^2}),{\sigma }_{\theta }=\frac{a^{2}P}{b^2-a^2}(\frac{b^2}{r^2}+1)$$

(54)

where a and b are the inner and the outer diameters, respectively; r is the radius of the stress point; and ${\sigma }_r$ and ${\sigma }_{\theta }$ are, respectively, the radial stress and circumferential stress. Equation (54) shows that there is no effect of the elastic modulus to the final stresses. Thus, the interval radius of the stresses in Table 4 which is at a very low level is theoretically reasonable.

Elastoplastic interval analysis

The yield criterion of the material is Mises and the yield stress is σ_s = 260 MPa. The load P is divided into 20 increments in order to model the evolutionary process of the nonlinear problems. According to the elastoplastic analytical solutions of the cylinder, final results can be given by

$$u_r=\frac{1}{\sqrt{3}}\frac{1.5r_s^2{\sigma }_s}{\mathrm{Er}}$$

(55)

$${\sigma }_r=\frac{2}{\sqrt{3}}{\sigma }_s\ln \frac{r}{a}-P,{\sigma }_{\theta }=\frac{2}{\sqrt{3}}{\sigma }_s(1+\ln \frac{r}{a})-P$$

(56)

Figure 7 shows the load–displacement curve in y-direction of point A given by the proposed nonlinear ISM and the analytical method. Figure 8 shows the load–stress curve of the first integral point in element 6 given by the proposed nonlinear ISM and the analytical method.

Figure 7.

Load–displacement curve of $y_{\mathbf{A}}$ .

Figure 8.

Load–stress curve.

The difference between the two displacements in Figure 7 is mainly caused by the error of the middle-point value given by the two methods. In this case, $\underset{¯}{E}=0.9E_c$ and $\overline{E}=1.1E_c$ , according to equation (55), bounds of the displacements are $\underset{¯}{\mathbf{u}}={\mathbf{u}|}_{E=E_c}/1.1$ and $\overline{\mathbf{u}}={\mathbf{u}|}_{E=E_c}/0.9$ , respectively. Thus, the middle-point value given by the analytical method is ${\mathbf{u}}_c=1.01{\mathbf{u}|}_{E=E_c}$ , while the proposed method is just ${\mathbf{u}|}_{E=E_c}$ . It means the results given by the analytical method are a little bigger than the nonlinear ISM and this situation corresponds to Figure 7.

Since the interval stresses are based on the interval displacements, differences of the displacements result in the error of the stresses provided by different methods in Figure 8. Besides, the interval expansion of interval arithmetic increases the stresses in the plastic stage. Actually, the stresses calculated by equation (56) are not affected by elastic modulus.

Overall, the interval results obtained by the nonlinear ISM are very close to the analytical ones shown in Figures 7 and 8. This indicates that using the proposed nonlinear ISM to evaluate the range of interval elastoplastic problems is quite feasible.

Conclusion

An IFEM using substructure method has been presented to determine the responses for solid structures with uncertain parameters. The dependency of parameters is considered. Elements of the same material are treated as a substructure, and then, the substructure method is used to decrease the dependency phenomenon. The Taylor expansion function is used to provide more narrow stresses. Then, the proposed ISM is extended to nonlinear problems. A constant stiffness method is applied to simplify the elastoplastic iterations. Then, the difference method is adapted to updating the interval stresses. The feasibility and rationality of the proposed method has been demonstrated by numerical examples.

For linear problems, it is found that the results given by the proposed Sub-IFEM is much closer to the vertex solutions compared to results provided by the EBE-IFEM. Furthermore, the Sub-IFEM can greatly reduce the number of equations. When extending to nonlinear problems, the substructure model can be integrated effectively with the incremental theory. The elastoplastic example indicates that the results obtained by the nonlinear ISM are mostly the same with analytical solutions. It can be seen that the nonlinear interval method will still work on other nonlinear problems, such as the nonlinear elastic problems and the viscoelasticity problems. It should be pointed out that in this article, the evaluation of the yield state is totally based on the middle point of the results. For interval problems, it needs deeper researches.