An improved interval finite element method based on the element-by-element technique for large truss system and plane problems

Abstract

The interval finite element method based on the element-by-element technique is proved to be a rigorous and efficient method for considering interval uncertainties. The structural internal force of the same accuracy as the interval response displacement can be obtained by combining the variational method and the Lagrange multiplier method. However, the increase of the Lagrange multiplier reduces the efficiency of solving the equation, especially for a large truss system. In addition, the method is mostly used to analyse truss system, and the applications in plane problems are still rare. In this article, an improved method involving modification of the Lagrange multiplier interval equations is proposed to calculate interval displacements and interval internal forces at the same time. In this proposed method, the increase of the unknown Lagrange interval does not affect the solution of the equations. For plane problems, a new method for solving the stress interval of the element is proposed. This proposed method can effectively reduce the expansion of the stress interval of the element. Finally, the correctness and rationality of the proposed method are verified by numerical examples.

Introduction

By combining the interval analysis method ¹ with the finite element method, a new method for interval uncertainties, the interval finite element method,^2,3 has been proposed. The main idea of the interval finite element method is to study the influence of the structural system parameters on the structural displacement and stress by considering the structural parameters as interval variables; its core content is to solve the interval governing equations. ⁴ Based on the monotonicity of the function, the combination monotone method was proposed by Rao and Berke. ² However, this form of monotonicity is not always satisfied, and with more uncertain parameters, the calculation of this method is more complicated. A perturbation approach of the interval finite element method based on the perturbation theory was proposed by Qiu and Elishakoff ⁵ and Wang et al.; ⁶ however, their method can solve only small perturbation problems. Subsequently, a sub-interval perturbation method for large perturbation cases was proposed; ⁵ however, although the increase of the sub-interval improves the accuracy of the results, it also requires a large amount of calculation. The interval finite element method has a wide application prospect in engineering, but the interval extension caused by the interval correlation problem restricts the application of the method. 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. In the work of Muhanna et al.,^8–10 based on the element-by-element (EBE) technique, the structural stiffness matrix is assembled, and the correlation problem caused by element stiffness coupling is fully considered by the interval factor separation method. The EBE-based interval finite element method was also extended to calculate the envelope frequency response functions with uncertain parameters by Yaowen et al.,^11,12 and the results indicated the effectiveness of the EBE model to frame structures. To obtain the high-precision axial force of a truss, a Lagrange multiplier method was proposed by Muhanna et al., ⁹ the method provides the same accuracy axial force of truss as the interval displacement. However, the increased number of unknown quantities decreases the effectiveness of the method, especially for a large truss system. So far, the interval finite element method, based on the EBE technique, is mostly used to analyse a frame structure, whereas for plane problems, the use of the EBE technique is seldom mentioned. Moreover, compared with the displacement, the element stress in plane problems is of more concern. The method of element interval stress in plane problems was given by Zhang. ¹³ The method is effective for small uncertainty, but the degree of the extension is larger when the uncertainty is larger.

In this article, an improved interval finite element method based on the EBE technique is proposed based on the full study of the finite element method based on the EBE technique. After improvement, interval internal forces can be calculated with the same accuracy as the interval displacements, eliminating the effect of variation of the unknown quantities to the calculation process. A new formulation to obtain the element stress is proposed for the condition in which the accuracy of stress is low when the parameter uncertainties are fairly large.

Interval conceptions

The interval number is defined by real and limitary set $[\underset{¯}{x},\overline{x}]=\{x|\le x\le \overline{x}\}$ . As a rational number $a/b$ or complex number is regarded as a pair of in-sequence integer, the interval number is regarded as a pair of in-sequence real numbers $\underset{¯}{x}$ and $\overline{x}$ . The interval number is denoted by $X^I$ , and $\underset{¯}{x}$ and $\overline{x}$ are its two limits, that is, $X^I=[\underset{¯}{x},\overline{x}]$ . For two interval numbers $X^I=[\underset{¯}{x},\overline{x}]$ and $Y^I=[\underset{¯}{y},\overline{y}]$ , the operation laws are defined as

$$\{\begin{array}{l}{X}^I+Y^I=\left[\underset{\_}{x},\overline{x}\right]+\left[\underset{\_}{y},\overline{y}\right]=\left[\underset{\_}{x}+\underset{\_}{y},\overline{x}+\overline{y}\right]\\ X^I-Y^I=\left[\underset{\_}{x},\overline{x}\right]-\left[\underset{\_}{y},\overline{y}\right]=\left[\underset{\_}{x}-\underset{\_}{y},\overline{x}-\overline{y}\right]\\ X^I\times Y^I=\left[\underset{\_}{x},\overline{x}\right]\times \left[\underset{\_}{y},\overline{y}\right]=\left[\min \left(\underset{\_}{x}\underset{\_}{y},\overline{x}\overline{y},\underset{\_}{x}\overline{y},\overline{x}\underset{\_}{y}\right),\max \left(\underset{\_}{x}\underset{\_}{y},\overline{x}\overline{y},\underset{\_}{x}\overline{y},\overline{x}\underset{\_}{y}\right)\right]\\ \frac{X^I}{Y^I}=\frac{\left[\underset{\_}{x},\overline{x}\right]}{\left[\underset{\_}{y},\overline{y}\right]}=\left[\underset{\_}{x},\overline{x}\right]\times \left[\frac{1}{\overline{y}},\frac{1}{\underset{\_}{y}}\right],0\notin \left[\underset{\_}{y},\overline{y}\right]\end{array}$$

(1)

The operation laws for a real number, such as commutative law, associative law and identity law, are only partly tenable for an interval number. Other operation laws represent infirm inclusive form, such as distributive law and quits law, for example, $X^I\times (Y^I\pm Z^I)\subseteq X^I\times Y^I\pm X^I\times Z^I$ . Thus, the operations of interval number may be induced in the extension of interval. The interval extension is described by the relations of real-valued function f and interval-valued function F

$$\{f(x_1,x_2,\dots ,x_n)x_i\in X_i^I,i=1,2,\dots ,n\}\subseteq \mathrm{F}(X_1^I,X_2^I,\dots ,X_n^I)$$

(2)

For interval extension, some methods have been proposed, such as the combination monotone method, ² the sub-interval perturbation method ⁵ and the truncation method. ¹⁴

Interval finite element method based on EBE

Generally, the interval finite element equation for the static boundary value problems is described as

$${\mathit{K}}^I{\mathit{U}}^I={\mathit{R}}^I\mathrm{or}\mathit{K}({\mathit{a}}^I){\mathit{U}}^I=\mathit{R}{\left(\mathit{a}\right)}^I$$

(3)

where $\mathit{K}$ is the stiffness matrix, $\mathit{R}$ is the load vector, ${\mathit{a}}^I$ is the interval parameter vector and ${\mathit{U}}^I$ is the static interval displacement vector.

For the interval finite element method based on EBE, ¹⁵ it is assumed that the interval parameters between each element are independent, and the interval parameters in the element are interrelated. By separating the elements and using the EBE technique, the structural global stiffness matrix is assembled to eliminate the correlation of adjacent units during the coupling process in the common junction. Supposing the uncertainty of the elastic modulus is ${\delta }^I$ , that is, ${\mathit{E}}^I=E_c+E_c{\delta }^I$ , the interval element stiffness matrix is

$${\mathit{K}}^I={\mathit{K}}_c+{\mathit{K}}_c{\mathit{D}}^I$$

(4)

where ${\mathit{K}}_c$ has the form of ${\mathit{K}}_c=\left(\begin{array}{ccc}{\mathit{K}}_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\mathit{K}}_{\mathrm{Ne}}\hfill \end{array}\right)$ , and ${\mathit{D}}^I$ is an interval diagonal matrix containing $d_i^I$ , having the form of ${\mathit{D}}^I=\left(\begin{array}{ccc}{\mathit{d}}_1^I\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& \ddots \hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\mathit{d}}_{\mathrm{Ne}}^I\hfill \end{array}\right)$ . ${\mathit{d}}_i^I$ is an interval diagonal matrix containing the interval quantity ${\delta }_i^I$ , having the form of ${\mathit{d}}_{\mathrm{ie}}^I=\mathrm{diag}({\delta }_1^I\cdots {\delta }_n^I)$ , $i=1,2,\dots ,\mathrm{Ne}$ , and $\mathrm{Ne}$ is the total number of the elements, and n is the number of degrees of freedom of each element.

The assembled stiffness matrix is singular because the mutual restraints between elements are released by the separation of elements. To remedy the displacement compatibility condition after the separation of the EBE system, additional constraints should be imposed. Suppose the compatibility of displacements is as follows

$$\mathit{c}\mathit{u}=0$$

(5)

where $\mathit{c}$ is a node-selection matrix with only 0, 1 and −1. For non-interval conditions, the total potential energy functional of the discrete model structure constructed by the Lagrange multiplier method is described as

$${\mathrm{\Pi }}_p^{*}=\frac{1}{2}{\mathit{u}}^T\mathit{Ku}-{\mathit{u}}^T\mathit{P}+\mathbf{\lambda }\mathit{cu}$$

(6)

where ${\mathrm{\Pi }}_p^{*}$ is the total potential energy functional, $\mathit{K}$ is the global stiffness matrix, and $\mathit{P}$ is the structural load array. $\mathbf{\lambda }$ is the Lagrange multiplier vector, having the form of $\mathbf{\lambda }=({\mathrm{\lambda }}_1\dots {\mathrm{\lambda }}_{\mathbf{m}}{)}^{\mathbf{T}}$ , and $\mathrm{m}$ is the number of additional constraint equations. Through the stationary value of the potential energy function, that is, $\partial {\mathrm{\Pi }}^{*}/\partial \mathit{u}$ , the governing equations ⁹ are obtained

$$\left(\begin{array}{cc}\mathit{K}\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& 0\hfill \end{array}\right)\left(\begin{array}{c}\mathit{u}\hfill \\ \multicolumn{1}{c}{\mathbf{\lambda }}\end{array}\right)=\left(\begin{array}{c}\mathit{p}\hfill \\ \multicolumn{1}{c}{0}\end{array}\right)$$

(7)

Considering the uncertainty of the parameters, we have

$$\left(\begin{array}{cc}{\mathit{K}}^I\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& 0\hfill \end{array}\right)\left(\begin{array}{c}{\mathit{u}}^I\hfill \\ \multicolumn{1}{c}{{\mathbf{\lambda }}^I}\end{array}\right)=\left(\begin{array}{c}{\mathit{p}}^I\hfill \\ \multicolumn{1}{c}{0}\end{array}\right)$$

(8)

Define

$$\begin{array}{c}{\mathit{A}}^I=\left(\begin{array}{cc}{\mathit{K}}^I\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& 0\hfill \end{array}\right)=\left(\begin{array}{cc}{\mathit{K}}_c\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& 0\hfill \end{array}\right)\hfill \\ \multicolumn{1}{c}{+\left(\begin{array}{cc}{\mathit{K}}_c\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right)\left(\begin{array}{cc}{\mathit{D}}^I\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right)={\mathit{A}}_c+{\mathit{S}}_c{\mathit{D}}^I}\\ {\mathit{x}}^I=\left(\begin{array}{c}{\mathit{u}}^I\hfill \\ \multicolumn{1}{c}{{\mathbf{\lambda }}^I}\end{array}\right),{\mathit{b}}^I=\left(\begin{array}{c}{\mathit{p}}^I\hfill \\ \multicolumn{1}{c}{0}\end{array}\right)\hfill \end{array}$$

Thus, equation (8) is transformed into

$${\mathit{A}}^{\mathit{I}}{x}^{\mathit{I}}=b^{\mathit{I}}$$

(9)

According to the fixed-point iterative method,^16,17 the solution of the equations is

$${\mathit{x}}^{\mathit{I}}=\sum (A^{\mathit{I}},b^{\mathit{I}})=x_{\mathcal{c}}+{\left(x^{*\mathit{I}}\right)}^{\mathcal{l}+1}$$

(10)

$${\left({\mathit{x}}^{*I}\right)}^{l+1}=\mathit{R}{\mathit{b}}^{\mathit{I}}-{\mathit{x}}_c-\mathit{R}{\mathit{S}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(11)

where ${\mathit{x}}_c=\mathit{Rb}$ , $\mathit{R}=({\mathit{A}}_c{)}^{-1}$ , ${\mathit{A}}_c$ is the midpoint of the interval stiffness matrix ${\mathit{A}}^{\mathit{I}}$ , and ${\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}={\mathit{D}}^{\mathit{I}}({\mathit{u}}_c+({\epsilon }^{\mathit{I}}\Delta {\mathit{u}}^{\mathit{I}}{)}^l)$ . ${\epsilon }^{\mathit{I}}=[1-\beta ,1+\beta ]$ is the interval expansion factor, which is used to improve the efficiency of interval iteration, and $\beta$ is usually between 0.001 and 0.1.

The response displacement interval and the Lagrange multiplier interval of the structure can be separated from ${\mathit{x}}^{\mathit{I}}$ using the Boolean matrix $\mathit{L}$

$${\mathit{u}}^{\mathrm{I}}={\mathit{L}}_1{\mathit{x}}^{\mathrm{I}}={\mathit{L}}_1\mathit{R}{\mathit{b}}^{\mathrm{I}}-{\mathbf{L}}_1\mathit{R}{\mathit{S}}_{\mathrm{c}}{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathrm{I}}$$

(12)

$${\mathbf{\lambda }}^{\mathit{I}}={\mathit{L}}_2\mathit{R}{\mathit{b}}^{\mathit{I}}-{\mathit{L}}_2\mathit{R}{\mathit{S}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(13)

where the Boolean matrix $\mathit{L}$ is either 0 or 1, indicating the relationship between the global stiffness matrix and the element stiffness matrix.

Improved interval finite element method based on EBE

Note that the structural stiffness matrix ${\mathit{K}}_c$ given by equation (8) is singular, and ${\mathit{A}}_c$ contains the items with the principal element of 0. Obviously, ${\mathit{A}}_c$ is also a singular matrix, and it does not exist for the inverse matrix ${\mathit{A}}_c^{-1}$ . In addition, the original solution method of equation (8) not only increases the number of degrees of freedom but also expands the size of the matrix, thereby increasing the storage burden of the computer and greatly reducing the computational efficiency, especially when the scale of the problem is large. In view of the above shortcomings, an improved solution method of the Lagrange multiplier interval equations is proposed in this article.

The displacement interval compatibility condition is transformed as

$$\mathit{c}{\mathit{u}}^{\mathit{I}}=0\Rightarrow \mathit{c}{\mathit{u}}^{\mathit{I}}-\frac{I}{\alpha }{\mathbf{\lambda }}^{\mathit{I}}=-\frac{I}{\alpha }{\mathbf{\lambda }}^{\mathit{I}}$$

(14)

Substituting equation (14) into equation (8), we obtain

$$\left(\begin{array}{cc}{\mathit{K}}^{\mathit{I}}\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& -\frac{\mathit{I}}{\alpha }\hfill \end{array}\right)\left(\begin{array}{c}{\mathit{u}}^{\mathit{I}}\hfill \\ \multicolumn{1}{c}{{\mathbf{\lambda }}^{\mathit{I}}}\end{array}\right)=\left(\begin{array}{c}{\mathit{p}}^{\mathit{I}}\hfill \\ \multicolumn{1}{c}{-\frac{I}{\alpha }{\mathbf{\lambda }}^{\mathit{I}}}\end{array}\right)$$

(15)

where $\alpha$ is a specified constant, which is advisable to be $\mathrm{\eta }{k}_{\max }$ . $k_{\max }$ is the maximum principal component in the stiffness matrix, and $\mathrm{\eta }$ is the penalty factor. At this point, both sides of equation (15) have the term ${\mathbf{\lambda }}^{\mathit{I}}$ , which can be solved via the iterative method, that is

$$\left(\begin{array}{cc}{\mathit{K}}^{\mathit{I}}\hfill & {\mathit{c}}^T\hfill \\ \multicolumn{1}{c}{\mathit{c}}& -\frac{I}{\alpha }\hfill \end{array}\right)\left(\begin{array}{c}{\left({\mathit{u}}^{\mathit{I}}\right)}^{n+1}\hfill \\ \multicolumn{1}{c}{{\left({\mathbf{\lambda }}^{\mathit{I}}\right)}^{n+1}}\end{array}\right)=\left(\begin{array}{c}{\mathit{p}}^{\mathit{I}}\hfill \\ \multicolumn{1}{c}{-\frac{I}{\alpha }{\left({\mathbf{\lambda }}^{\mathit{I}}\right)}^n}\end{array}\right)$$

(16)

The iterative process is expressed as follows

$$({\mathit{K}}^{\mathit{I}}+\mathrm{\alpha }{c}^T\mathit{c}){\left({\mathit{u}}^{\mathit{I}}\right)}^{n+1}={\mathit{p}}^{\mathit{I}}-{\mathit{c}}^T{\left({\mathbf{\lambda }}^{\mathit{I}}\right)}^n$$

(17)

$${\left({\mathbf{\lambda }}^{\mathit{I}}\right)}^{n+1}={\left({\mathbf{\lambda }}^{\mathit{I}}\right)}^n+\alpha \mathit{c}{\left({\mathit{u}}^{\mathit{I}}\right)}^{n+1}$$

(18)

Assume $\mathit{Q}=\mathrm{\alpha }{\mathit{c}}^T\mathit{c}$ , ${\mathit{A}}^{\mathit{I}}=({\mathit{K}}^{\mathit{I}}+\mathit{Q})$ , ${\mathit{b}}^{\mathit{I}}={\mathit{p}}^{\mathit{I}}-{\mathit{c}}^T({\mathbf{\lambda }}^{\mathit{I}}{)}^n$ , ${\mathit{u}}^{\mathit{I}}={\mathit{u}}_c+{\mathit{u}}^{*I}$ , ${\mathbf{\lambda }}^{\mathit{I}}={\mathbf{\lambda }}_c+{\mathbf{\lambda }}^{*I}$ and ${\mathit{b}}^{\mathit{I}}={\mathit{b}}_c+{\mathit{b}}^{*I}$ . Using the fixed-point iterative method for equations (17) and (18), we obtain

$${\left({\mathit{u}}_c\right)}^{n+1}=\mathit{R}{\mathit{b}}_c=\mathit{R}({\mathit{p}}_c-c^T{\left({\mathbf{\lambda }}_c\right)}^n)$$

(19)

$${\left({\mathbf{\lambda }}_c\right)}^{n+1}={\left({\mathbf{\lambda }}_c\right)}^n+\alpha \mathit{c}{\left({\mathit{u}}_c\right)}^{n+1}$$

(20)

$$\mathit{R}{\mathit{b}}^{\mathit{I}}-\mathit{R}{\mathit{A}}^{\mathit{I}}{\left({\mathit{u}}_c\right)}^{n+1}+(\mathit{I}-\mathit{R}{\mathit{A}}^{\mathit{I}}){\left(\epsilon {\left({\mathit{u}}^{*I}\right)}^{n+1}\right)}^l\subseteq {\left({\left({\mathit{u}}^{*I}\right)}^{n+1}\right)}^{l+1}$$

(21)

$${\left({\mathbf{\lambda }}^{*I}\right)}^{n+1}={\left({\mathbf{\lambda }}^{*I}\right)}^n+\mathit{a}\mathit{c}{\left({\left({\mathit{u}}^{*I}\right)}^{n+1}\right)}^{l+1}$$

(22)

where ${\mathit{A}}_c$ is not a singular matrix, $\mathit{R}={\mathit{A}}_c^{-1}=({\mathit{K}}_c+\alpha {\mathit{c}}^T\mathit{c}{)}^{-1}$ .

If only the uncertainty of the elastic modulus of the material is taken into account and load $\mathit{p}$ is the determined value, we obtain ${\mathit{b}}^{*I}=-{\mathit{c}}^T({\mathbf{\lambda }}^{*I}{)}^n$ . Thus, the interval displacement in equations (21) and (22) can be expressed by

$${\left({\left({\mathit{u}}^{*I}\right)}^{n+1}\right)}^{l+1}=-\mathit{R}{\mathit{c}}^T{\left({\mathbf{\lambda }}^{*I}\right)}^n-\mathit{R}{\mathit{K}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(23)

The initial conditions for the above iterative computations are

$$n=0:{\mathbf{\lambda }}_c=0,{\mathbf{\lambda }}^{*I}=0$$

(24)

The converge condition is $\Vert {\left({\mathbf{\lambda }}_{\mathrm{c}}\right)}^{\mathrm{n}+1}-{\left({\mathbf{\lambda }}_{\mathrm{c}}\right)}^{\mathrm{n}}\Vert \le \Vert \mathrm{Toler}{\left({\mathbf{\lambda }}_{\mathrm{c}}\right)}^{\mathrm{n}}\Vert$ . Where ‖.‖ is a matrix norm, and $\mathrm{Toler}$ is the convergence tolerance, which can be chosen as 1.0 × 10⁻⁶.

According to numerical experiences, the iteration equations (19)–(22) converge after only two steps; at this time, n equals to 1. According to equations (19) and (20), the Lagrange multiplier solution of the equations is $({\mathbf{\lambda }}_c{)}^n$ , and the displacement solution is $({\mathit{u}}_c{)}^{n+1}$ . When n is equal to 1

$${\left({\mathit{u}}_c\right)}^1=\mathit{R}{\mathit{b}}_c=\mathit{R}({\mathit{p}}_c-{\mathit{c}}^T{\left({\mathbf{\lambda }}_c\right)}^0)$$

(25)

$${\left({\mathbf{\lambda }}_c\right)}^1={\left({\mathbf{\lambda }}_c\right)}^0+\alpha \mathit{c}{\left({\mathit{u}}_c\right)}^1$$

(26)

$${\left({\mathit{u}}_c\right)}^2=\mathit{R}({\mathit{p}}_c-{\mathit{c}}^T{\left({\mathbf{\lambda }}_c\right)}^1)$$

(27)

The initial multiplier vector $({\mathbf{\lambda }}_c{)}^0=0$ ; thus, equations (26) and (27) can be transformed into

$${\mathbf{\lambda }}_c=\alpha \mathit{c}\mathit{R}{\mathit{p}}_c$$

(28)

$${\mathit{u}}_c=\mathit{R}{\mathit{p}}_c-{\mathit{c}}^T{\mathbf{\lambda }}_c$$

(29)

In the same manner, the interval part equations (22) and (23) can be changed to

$${\mathbf{\lambda }}^{*I}=-\alpha \mathit{c}\mathit{R}{\mathit{K}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(30)

$${\left({\mathit{u}}^{*I}\right)}^{l+1}=-\mathit{R}{\mathit{c}}^T{\mathbf{\lambda }}^{*I}-\mathit{R}{\mathit{K}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(31)

When the interval Lagrange multiplier $-\mathit{R}{\mathit{c}}^T{\mathbf{\lambda }}^{*I}$ has little effect on the response displacement interval, the interval displacement can be approximately taken as

$${\left({\mathit{u}}^{*I}\right)}^{l+1}=-\mathit{R}{\mathit{K}}_c{\mathit{M}}^{\mathrm{Il}}{\Delta }^{\mathit{I}}$$

(32)

The response displacement interval and the Lagrange multiplier interval of the structure can be obtained easily and efficiently using the improved interval finite element method based on EBE. The Lagrange multiplier interval can be converted for the element nodal force interval. In the displacement interval finite element, the nodal force interval of element is described as

$${\mathit{F}}_e^{\mathit{I}}={\mathit{L}}_e({\mathit{K}}^{\mathit{I}}{\mathit{u}}^{\mathit{I}}-{\mathit{p}}_b^{\mathit{I}})$$

(33)

where ${\mathit{L}}_e$ is the element selection matrix, ${\mathit{F}}_e^{\mathit{I}}$ is the nodal force interval of element, and ${\mathit{p}}_b^{\mathit{I}}$ is a structural boundary support reaction interval; for a truss structure, ${\mathit{p}}_b^{\mathit{I}}$ is the constraint reaction interval. According to equation (15), we have

$${\mathit{K}}^{\mathit{I}}{\mathit{u}}^{\mathit{I}}={\mathit{P}}^{\mathit{I}}-{\mathit{c}}^T{\mathbf{\lambda }}^{\mathit{I}}$$

(34)

Substituting equation (34) into equation (33) yields

$${\mathit{F}}_e^{\mathit{I}}={\mathit{L}}_e({\mathit{P}}^{\mathit{I}}-{\mathit{c}}^T{\mathbf{\lambda }}^{\mathit{I}}-{\mathit{p}}_b^{\mathit{I}})$$

(35)

In system of bars, equation (35) provides a method for effectively calculating the axial force of a member bar.

A new interval stress solution method

In the traditional finite element method, the nodal force of element can be obtained by the element stress

$${\mathit{F}}^e=\underset{\mathrm{\Omega }}{\int }{\mathit{B}}^T{\mathbf{\sigma }}_{e}d\mathrm{\Omega }=\left(\underset{\mathrm{\Omega }}{\int }{\mathit{B}}^T\mathit{D}\mathit{B}d\mathrm{\Omega }\right){\mathit{u}}_e$$

(36)

where $\mathit{B}$ is the element strain matrix, $\mathit{D}$ is the elastic matrix, ${\mathbf{\sigma }}_e$ is the element stress array, and ${\mathit{u}}_e$ is the element displacement array.

For a three-node triangular element, the strain matrix $\mathit{B}$ and the element stress ${\mathbf{\sigma }}_e$ are constant. Thus, equation (36) can be provided as

$${\mathit{F}}^e={\mathit{B}}^T{\mathbf{\sigma }}_e\mathrm{St}$$

(37)

where $S$ is the plane element area, and $t$ is the plane element thickness.

For a quadrilateral element, the element integral form of the natural coordinate system is obtained via equation (36) as

$${\mathit{F}}^e=\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathit{B}}^T{\mathbf{\sigma }}_e\left|\mathit{J}\right|\mathrm{td}\epsilon d\eta$$

(38)

For rectangular or parallelogram elements, the Jacobian matrix $\left|\mathit{J}\right|$ is a constant; for general four-node elements, the Jacobian matrix $\left|\mathit{J}\right|$ can be approximated as a constant when the element size is small. Suppose $\mathit{\phi }(\mathbf{\epsilon },\mathbf{\eta })={\mathit{B}}^{\mathbf{T}}{\mathbf{\sigma }}_e$ and the equation (38) can be transformed as

$${\mathit{F}}^e=\left|\mathit{J}\right|t\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathit{B}}^T{\mathbf{\sigma }}_{e}d\epsilon d\eta$$

(39)

Because the stress inside the element is continuous, that is, $\mathit{\phi }(\mathbf{\epsilon },\mathbf{\eta })$ is continuous in the element, according to the double integral mean theorem, it is obtained the following equation

$${\mathit{F}}^e=\mathit{\phi }(\alpha ,\beta )\left|\mathit{J}\right|t\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}d\epsilon d\eta =4t\mathit{\phi }(\alpha ,\beta )\left|\mathit{J}\right|$$

(40)

where $-1\le \alpha \le 1$ , $-1\le \beta \le 1$ and $\mathit{\phi }(\mathbf{\alpha },\mathbf{\beta })$ are the stress at coordinate $(\alpha ,\beta )$ within the element. The approximate and exact solution of stress is equal at the optimal stress point of the plane four-node isoparametric element. ¹⁸ The optimal stress point of the four-node rectangular or parallelogram element is located at the centroid of the element, whose natural coordinates are $\epsilon =0$ and $\eta =0$ . Gaussian numerical integration is typically used for the four-node isoparametric element, and the accuracy of the integral depends on the order of the Gaussian integral. The exact integral scheme of the four-node parallelogram element is $2\times 2$ order Gaussian integral, that is, the element nodal force expression given by equation (27) is precisely established. Therefore, for parallelogram or rectangular elements, using the $2\times 2$ order Gaussian integral scheme, the following equation is obtained

$${\mathit{F}}^e=4t\mathit{\phi }(0,0)\left|\mathit{J}\right|={\left({\mathit{B}}^{\mathbf{T}}{\mathbf{\sigma }}_e\right)}_{(0,0)}4\mathbf{t}\left|\mathit{J}\right|$$

(41)

Equations (37) and (40) can be uniformly expressed as

$${\mathit{F}}^e={\mathit{B}}^T{\mathbf{\sigma }}_e\mathit{V}$$

(42)

where

$$\{\begin{array}{c}\Delta :V=\mathrm{St}\hfill \\ \multicolumn{1}{c}{\square :V=4\left|\mathit{J}\right|t=2\mathrm{St}}\end{array}$$

(43)

Because $\mathit{B}$ is an irregular matrix, there is no inverse matrix for $\mathit{B}$ ; however, the generalized inverse matrix ${\mathit{B}}^{+}$ exists and is unique. If the element nodal force is known, then equation (42) can be regarded as an overdetermined equation of the element stress, and the least square solution of the element stress is

$$\mathbf{\sigma }={\left({\mathit{B}}^{\mathbf{T}}\right)}^{+}{\mathit{F}}^e/\mathbf{V}$$

(44)

Considering the uncertainty of the structural parameters, the interval stress of the element is

$${\mathbf{\sigma }}_e^{\mathit{I}}={\left({\mathit{B}}^{\mathbf{T}}\right)}^{+}{\mathit{F}}^{\mathbf{Ie}}/\mathbf{V}$$

(45)

Substituting equation (35) into equation (45) yields

$$\begin{array}{c}{\mathbf{\sigma }}_e^{\mathit{I}}={\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e({\mathit{P}}^{\mathit{I}}-{\mathit{c}}^T{\mathbf{\lambda }}^{\mathit{I}}-{\mathit{p}}_b^{\mathit{I}})/V\hfill \\ =\{{\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e{\mathit{P}}^{\mathit{I}}-{\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e{\mathit{c}}^T{\mathbf{\lambda }}^{\mathit{I}}-{\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e{\mathit{p}}_b^{\mathit{I}}\}/V\hfill \end{array}$$

(46)

From the properties of the generalized inverse matrix, we obtain $({\mathit{B}}^T{)}^{+}=(B^{+}{)}^T$ . Because $\mathrm{rank}(\mathit{B})=m=3$ , that is, $\mathit{B}$ is a row-full-rank matrix, we obtain $({\mathit{B}}^T{)}^{+}=({\mathit{B}}^T(\mathit{B}{\mathit{B}}^T{)}^{-}{)}^T$ . If the Lagrange multiplier method is also used for model constraints, then ${\mathit{p}}_b^{\mathit{I}}$ can be given by the Lagrange multiplier. Substituting the Lagrange multiplier interval into equation (46), the element stress interval can be obtained as

$${\mathbf{\sigma }}_e^{\mathit{I}}=\{{\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e({\mathit{P}}^{\mathit{I}}-{\mathit{c}}^T{\mathbf{\lambda }}^{\mathit{I}})-{\left({\mathit{B}}^T\right)}^{+}{\mathit{L}}_e\mathit{Q}\mathit{R}{\mathit{K}}_c{\mathit{M}}^{\mathit{I}}{\Delta }^{\mathit{I}}\}/V$$

(47)

A reasonable element stress interval can be obtained by delaying the interval multiplication operation in the calculation process as long as possible. The method is strict and effective for triangular and rectangular or parallelogram elements, and the stress interval for the general quadrilateral element is an approximation.

Numerical examples

Plane truss example

Here, a plane truss ¹³ is considered; the truss is composed of m-bar –n-layer members. The left node of the truss is subjected to a horizontal concentrated load, and the upper node is subjected to a vertical concentrated load, as shown in Figure 1, where the load is 10 kN, elastic modulus of bars E = 210 GPa, the cross-sectional area is 0.0025 m², and L is 1 m. Considering the uncertainty of the cross-sectional area of the element, the uncertainty is 1%. The different m and n values are calculated. The computer is running on the Windows 7 operating system using a processor operating at 1.6 GHz, and the memory is 1.5 GB.

Figure 1.

m-bar –n-layer truss structure.

Table 1 presents the comparison between the results obtained by the original method ¹³ and the improved method of the horizontal x-displacements at point D, and Table 2 presents the corresponding comparison of the axial forces at element 1. Tables 1 and 2 show that the results obtained by the improved method are basically the same as the results of the original method, ¹³ and the error is very small, demonstrating that the improved method is reasonable and feasible.

Table 1.

Comparison of the x-directional displacement at point D (unit: 10⁻⁴ m).

m×n	Number of elements	Mean value			Interval radius
		Original method 13	Improved method	Error	Original method 13	Improved method	Error
2 × 2	18	1.679577	1.679577	0	0.009132	0.009132	0
2 × 5	42	10.18939	10.18939	0	0.054453	0.054453	0
3 × 5	63	6.165798	6.165802	6.49E–07	0.033661	0.033660	1.58E–05
3 × 8	99	22.67490	22.67489	4.41E–07	0.121946	0.121948	1.64E–05
3 × 10	123	45.80690	45.80679	2.40E–06	0.255417	0.255414	1.30E–05

Table 2.

Comparison of the axial force of element 1 (unit: kN).

m×n	Number of elements	Mean value			Interval radius
		Original method 13	Improved method	Error	Original method 13	Improved method	Error
2 × 2	18	5.597979	5.597979	0	0.065747	0.065747	0
2 × 5	45	34.191621	34.191581	1.17E–06	0.197073	0.197069	2.03E–05
3 × 5	65	12.234446	12.234432	1.14E–06	0.152604	0.152604	1.97E–05
3 × 8	104	57.461161	57.461153	1.39E–07	0.371831	0.371829	5.38E–06
3 × 10	130	100.06620	100.06620	0	0.684653	0.684653	0

Table 3 shows the comparison between the original method ¹³ and the improved method of the time consumption of the matrix inversion and interval iteration process, and Figure 2 shows the relationship curve of the central processing unit (CPU) time and the number of interval variables. It may be seen from Table 1 and Figure 2 that the computational time increases with the increase in the number of interval variables and the matrix inverse process takes up more computation time. The improved method reduces the storage scale of the matrix, thereby obviously improving the calculation efficiency.

Table 3.

Comparison of the CPU time of the computational process (unit: second).

m×n	Number of elements	Matrix inverse calculation		Interval iteration calculation
		Original method 13	Improved method	Original method 13	Improved method
2 × 2	18	0.094	0.031	0.016	0.016
2 × 5	42	1.047	0.281	0.250	0.125
3 × 5	63	3.250	0.656	0.797	0.328
3 × 8	99	14.594	1.750	3.125	0.859
3 × 10	123	27.031	2.766	5.641	1.375

Figure 2.

CPU time versus problem scale for the original method ¹³ and the improved method.

Plane quadrilateral element example

Consider the plate with a quarter-circle cutout ¹³ shown in Figure 3. The plate has dimensions of 0.1 m × 0.05 m with a thickness of 0.005 m, and the radius of the circular cutout is 0.02 m. Let Poisson’s ratio be v = 0.3 and the elastic modulus be E = 200 GPa. A uniformly distributed load of 100 kN/m is applied to the right edge, as shown in Figure 3. The plate is analysed as a plane-stress problem.

Figure 3.

Quadrilateral mesh model of plate with quarter-circle cutout.

Considering that the uncertainties of the elastic modulus of the material are 1%, 5% and 10%, Table 4 gives the comparison of the results obtained by the original method ¹³ and the improved method of the x-directional displacement at point A, Table 5 gives the comparison results of the y-directional displacement at point E, Table 6 gives the comparison results of the x-directional nodal force at node 14, and Table 7 gives the comparison results of the y-directional nodal force at node 14. It may be seen from Tables 4–7 that the results obtained by the improved method are basically the same as the results of the original method, ¹³ and the error is very small, demonstrating that the improved method is reasonable and feasible for the planar solid elements.

Table 4.

Comparison of the x-directional displacement at point A (unit: 10⁻⁵ m).

Method	Mean value	Interval radius
		1%	5%	10%
Original method 13	1.15398680	0.00597850	0.03401255	0.08188526
Proposed method	1.15398680	0.00597687	0.03401255	0.08188526
Error (%)	0	0.027	0	0

Table 5.

Comparison of the y-directional displacement at point E (unit: 10⁻⁵ m).

Method	Mean value	Interval radius
		1%	5%	10%
Original method 13	−0.37431038	0.00240586	0.01441726	0.03750738
Proposed method	−0.37431037	0.00240495	0.01441726	0.03750738
Error (%)	0	0.038	0	0

Table 6.

Comparison of the x-directional nodal force at node 14 (unit: kN).

Method	Mean value	Interval radius
		1%	5%	10%
Original method 13	−2.82719000	0.00184940	0.01227474	0.03485942
Proposed method	−2.82718990	0.00184939	0.01227474	0.03485941
Error (%)	0	0	0	0

Table 7.

Comparison of the y-directional nodal force at node 14 (unit: kN).

Method	Mean value	Interval radius
		1%	5%	10%
Original method 13	0.11838433	0.00096324	0.00697892	0.02179134
Proposed method	0.11838431	0.00096324	0.00697892	0.02179134
Error (%)	0	0	0	0

Considering that the uncertainties of the elastic modulus are 1%, 5%, 10% and 15%, Table 8 gives the comparison results obtained by the endpoint combination method, the original method ¹³ and the improved method of the centroid stress at the element 6. Figures 3 and 4 show the change curves of the stress interval width in the x-direction and the y-direction at element 6 as a function of the uncertainty of the elastic modulus for the three methods. Table 8 reveals the following: (1) the element interval stresses obtained by the original method ¹³ and the proposed method are less than those of the endpoint combination method, in agreement with the inherent characteristics of the two methods; (2) the results of the proposed method are superior to those of the original stress interval method; ¹³ and (3) with the increase of interval uncertainty, the result of the original method ¹³ is increasingly deviated from the reasonable range; on the contrary, the proposed method can give more reasonable element stress interval, as confirmed in Figures 4 and 5. Figures 4 and 5 also reveal that the improvement effect of the x-direction stress interval is more obvious than the y-direction stress interval because the x-direction is the direction of the load, and the stress value is greater than that in the y-direction; as a result, the magnitude of the medium value in the interval iterative method has a significant effect on the interval width.

Table 8.

Comparison of the element centroid stress at element 6 (unit: MPa).

Uncertainty (%)	Method	${\underset{¯}{\sigma }}_x$	${\overline{\sigma }}_x$	$\mathrm{wid}{\sigma }_x/{\sigma }_{\mathrm{xc}}$ (%)	${\underset{¯}{\sigma }}_y$	${\overline{\sigma }}_y$	$\mathrm{wid}{\sigma }_x/{\sigma }_{\mathrm{xc}}$ (%)
1	Combination method	17.404	17.516	0.64	1.569	1.642	4.57
	Original method 13	17.256	17.669	2.36	1.549	1.662	7.05
	Proposed method	17.402	17.522	0.69	1.566	1.644	4.90
5	Combination method	17.176	17.734	3.20	1.423	1.790	22.84
	Original method 13	16.308	18.653	13.42	1.264	1.961	43.21
	Proposed method	17.099	17.825	4.16	1.361	1.850	30.47
10	Combination method	16.890	18.010	6.42	1.237	1.970	45.71
	Original method 13	14.742	20.376	32.08	0.723	2.572	112.23
	Proposed method	16.527	18.397	10.71	0.950	2.260	81.66
15	Combination method	16.606	18.289	9.65	1.048	2.147	68.78
	Original method 13	12.489	23.028	59.35	0.203	3.637	223.64
	Proposed method	15.600	19.324	21.33	0.241	2.969	170.01

Figure 4.

Stress interval width of the x-direction versus the uncertainty of the elastic modulus.

Figure 5.

Stress interval width of the y-direction versus the uncertainty of the elastic modulus.

Conclusion

In this article, an improved interval finite element method based on the EBE technique was proposed. In the proposed method, the original governing equations of finite element method (FEM) are divided into two parts, namely, the displacement equations and the Lagrange multiplier equations, and the iterative process is used to solve the corresponding equations. With an increasing number of the unknown Lagrange multiplier equations, the process of solution of structural response displacement is not affected, and it can obtain the structural response interval displacement and the interval force with the same precision at the same time. The numerical examples showed that the improved method has the same accuracy as the original method, but it greatly improves the efficiency of the solution of the large-scale bar system.

A new method of solving the element stress was proposed by applying the improved method to the plane problem. In the proposed method, the element centroid stress interval can be solved by the element node force interval, and high accuracy of the element interval stress can be guaranteed. A numerical example showed that the improved method can maintain the interval accuracy, in agreement with the original method, ¹³ and even when the uncertainty of the structure is large, reasonable interval results can be obtained using the proposed method.