Interval element-free Galerkin method for uncertain mechanical problems

Abstract

An interval element-free Galerkin method was proposed to solve some issues in structural design and analysis of structural parameters that have errors or uncertainties caused by manufacture, installation, measurement, or computation. Based on the interval mathematics and perturbation theory, we deduced interval element-free Galerkin method equilibrium equations using the element-free Galerkin method and solved them using the parameter perturbation of interval number. A bar under uniform load, a bi-material cantilever beam subjected to uniform stress, and collinear double edged cracked plate subjected to uniform tensile stress, including uncertainty parameters, were analyzed. Numerical simulations show that interval element-free Galerkin method is accurate and effective in solving the uncertainty problems.

Introduction

Due to different reasons, real engineering structures always have some inevitable uncertainties, such as loads, physical and geometric parameters, boundary conditions, and system failure conditions.^1–5 Many calculation methods have been implemented to solve the uncertainty in solid mechanics,^6–8 structural mechanics,^9,10 fluid mechanics,^11,12 and fluid–structure interaction,^13,14 including stochastic finite element method (FEM),^15–18 fuzzy FEM,^19–21 and interval parameter perturbation method.^22–26 The stochastic FEM requires the statistic characteristics of uncertainty parameters, while the fuzzy FEM is very slow. In comparison, the interval analysis can quantificationally inspect the influence of uncertainty parameters without knowing their probability distribution and is very effective for small-range uncertainty problems. The interval parameter perturbation method, proposed by Wang and Qiu, ²⁷ has been widely applied in the structural response analysis due to its simplicity and efficiency. Wang and Qiu^28,29 developed the interval and subinterval perturbation methods to solve the heat convection–diffusion problem with uncertain-but-bounded parameters. A collocation method was proposed further in the article by Wang et al. ³⁰ to analyze the heat convection–diffusion problems with interval input parameters. Compared with previous interval approaches, the computational cost of the interval parameter perturbation method is smaller, and the convergence condition which is related to the ranges of interval parameters is more easily guaranteed.

Before each computation, FEMs should generate and maintain meshes of sufficient quality, which leads to a big amount of data preparations and requires complicated pre- and post-processes.^31–33 The meshes might suffer from severe distortion in metal stamping forming, high velocity impact and explosion, dynamic crack propagation, fluid–solid coupling, and other problems that involve the extreme large deformation, largely reducing the precision. As these defects of FEMs, the upsurge of researching element-free method^34–38 has begun since 1990s in the field of international computational mechanics.

The element-free method has many advantages such as strong anti-distortion ability, high precision, convenient post-process, elimination of volume self-locking, and quick convergence. This method has been developed fast and applied to solve problems about elastoplasticity, ³⁹ buckling, ⁴⁰ heat conduction, ⁴¹ elastodynamic, ⁴² large deformation, ⁴³ and micromechanical study on the electro-elastic behavior. ⁴⁴ However, its application on research about the uncertainty problems is rare, ⁴⁵ and especially, there is no research about the uncertainty problem of simulating the cracked structure. Ma et al. ⁴⁶ used the interval element-free method just to evaluate the bending defection range of plates on Winkler foundation. Interval element-free method needs to be further developed.

In the present study, we proposed the interval element-free Galerkin method (IEFGM) based on the investigation of the interval mathematics and element-free Galerkin method. This method only requires the node information, without demand for any element connectivity, and resolves the interval equilibrium formula using the interval number decomposing method. We also deduced the interval J-integration formula in detail. Finally, to verify the accuracy of IEFGM, we analyzed the issues of a bar under uniform load, a bi-material cantilever beam, and a collinear double edged cracked finite square plate under uniform tension with uncertainty parameters.

Moving least square approximation

Using the values of several unrelated nodes, moving least square (MLS) can achieve a fitting function that is smooth and has connective derivatives. Considering the uncertainty of material characteristics and loads, EFGM approximates the field function using an MLS-generated smoothing function. The function $u^h(\mathbf{x})$ , which is the approximating function $u(\mathbf{x})$ in the domain $\mathrm{\Omega }\subseteq {\Re }^K$ ( $K=1,2,3$ ), is defined as follows

$$u\left(x\right)\cong u^h\left(x\right)={\displaystyle \sum _{i=1}^m{p}_i\left(x\right)a_i\left(x\right)}=p^{\mathrm{T}}\left(x\right)a\left(x\right)$$ (1)

where ${\mathbf{p}}^{\mathrm{T}}(\mathbf{x})=\{p_1(\mathbf{x}),p_2(\mathbf{x}),\dots ,p_m(\mathbf{x})\}$ is a vector of complete basis functions of order m and $\mathbf{a}(\mathbf{x})=\{a_1(\mathbf{x}),a_2(\mathbf{x}),\dots ,a_m(\mathbf{x})\}$ is a vector of unknown parameters that depend on x. $p_i(\mathbf{x})$ are the basis functions and $a_i(\mathbf{x})$ are the coefficients to be determined.

In a two-dimensional space, a linear basis and a quadratic basis are respectively

$$p^{\mathrm{T}}\left(x\right)=\left\{\begin{array}{ccc}1& x& y\end{array}\right\},m=3\left(2\mathrm{D}\right)$$ (2)

$$p^{\mathrm{T}}\left(x\right)=\left\{\begin{array}{cccccc}1& x& y& x^2& xy& y^2\end{array}\right\},m=6\left(2\mathrm{D}\right)$$ (3)

Because of the singularity of $1/\sqrt{r}$ in the vicinity of the crack tip stress field and to improve the precision, we add a singular item ²¹ using the partial expansion basis

$$p^{\mathrm{T}}\left(x\right)=\left\{\begin{array}{cccc}1& x& y& \sqrt{r}\end{array}\right\}m=4$$ (4)

or the complete expansion basis

$$p^{\mathrm{T}}\left(x\right)=\left\{\begin{array}{ccccccc}1& x& y& \sqrt{r}\cos \left(\theta /2\right)& \sqrt{r}\sin \left(\theta /2\right)& \sqrt{r}\sin \left(\theta /2\right)\sin \theta & \sqrt{r}\cos \left(\theta /2\right)\sin \theta \end{array}\right\}m=7$$ (5)

where $r$ is the distance from one point to the crack tip, and $\theta$ is the angle between the straight line, connecting this point to the crack tip, and the crack line.

Coefficient vector $\mathbf{a}(\mathbf{x})$ is determined according to the weighted least square method, and the error function is defined as

$$\begin{array}{l}J\left(x\right)={\displaystyle \sum _{I=1}^{n}w\left(x-x_I\right){\left[u^h\left(x\right)-u\left(x_I\right)\right]}^2}\\ ={\displaystyle \sum _{I=1}^{n}w\left(x-x_I\right){\left[{\displaystyle \sum _{i=1}^m{p}_i\left(x_I\right)a_i\left(x\right)}-u_I\right]}^2}\end{array}$$ (6)

where n is the number of the nodes corresponding to the weighted function $w(\mathbf{x}-{\mathbf{x}}_I)\ne 0$ of point $\mathbf{x}$ and u_I is the nodal parameter for node I.

In order to minimize the error, we aim to minimize $J(\mathbf{x})$ through its extreme condition $\partial J(\mathbf{x})/\partial \mathbf{a}=0$

$$\frac{\partial J\left(x\right)}{\partial a}=A\left(x\right)a\left(x\right)-B\left(x\right)u=0$$ (7)

where

$$A\left(x\right)=P^{\mathrm{T}}WP={\displaystyle \sum _{I=1}^{n}w\left(x-x_I\right)p\left(x_I\right)p^{\mathrm{T}}\left(x_I\right)}$$ (8)

$$B\left(x\right)=P^{\mathrm{T}}W=\left[B_1\left(x\right),\mathrm{\dots },B_I\left(x\right),\mathrm{\dots },B_n\left(x\right)\right]$$ (9)

$$B_I\left(x\right)=w\left(x-x_I\right)p\left(x_I\right)$$ (10)

$$\mathbf{u}=\{u_1,u_2,\dots ,u_n\}$$ (11)

Hence, we have

$$\mathbf{a}\left(\mathbf{x}\right)={\mathbf{A}}^{-1}\left(\mathbf{x}\right)\mathbf{B}\left(\mathbf{x}\right)\mathbf{u}$$ (12)

Substituting equation (12) into equation (1), we have

$$u^h\left(\mathbf{x}\right)={\mathbf{\Phi }}^{\mathrm{T}}\left(\mathbf{x}\right)\mathbf{u}$$ (13)

where

$${\mathbf{\Phi }}^{\mathrm{T}}\left(\mathbf{x}\right)={\mathbf{p}}^{\mathrm{T}}\left(\mathbf{x}\right){\mathbf{A}}^{-1}\left(\mathbf{x}\right)\mathbf{B}\left(\mathbf{x}\right)$$ (14)

$\mathbf{\Phi }(\mathbf{x})$ as the shape function of EFGM is defined as follows

$${\mathbf{\Phi }}^{\mathrm{T}}\left(\mathbf{x}\right)=\{{\Phi }_1\left(\mathbf{x}\right),\dots ,{\Phi }_I\left(\mathbf{x}\right),\dots ,{\Phi }_n\left(\mathbf{x}\right)\}$$ (15)

$${\Phi }_I\left(\mathbf{x}\right)=\sum _{j=1}^{m}p\left(\mathbf{x}\right){\left[{\mathbf{A}}^{-1}\left(\mathbf{x}\right)\mathbf{B}\left(\mathbf{x}\right)\right]}_{\mathrm{jI}}$$ (16)

The t-distribution weight function is

$$w(\mathbf{x}-{\mathbf{x}}_I)=\{\begin{array}{cc}\frac{{(1+{\beta }^2{r}^2)}^{-(1+\beta /2)}-{(1+{\beta }^2)}^{-(1+\beta /2)}}{1-{(1+{\beta }^2)}^{-(\frac{1+\beta }{2})}},\hfill & \multicolumn{1}{c}{r\le 1}\\ \multicolumn{1}{c}{0,}& \multicolumn{1}{c}{r>1}\end{array}$$ (17)

where $r=d_I/d_{\mathrm{mI}}$ , $d_I=|\mathbf{x}-{\mathbf{x}}_I|$ , $d_{\mathrm{mI}}$ is the supporting radius of node $I$ , $d_{\mathrm{mI}}=\mathrm{scale}\times c_I$ , $c_I$ is the maximum distance between ${\mathbf{x}}_I$ and other nodes, $\mathrm{scale}$ is the scale factor, and β = 2. The t-distribution used in equation (17) represents the probability density function of a standard Gaussian random variable divided by the square root of a chi-squared random variable with β degrees of freedom. This weight function has been successfully used in solving various problems in linear–elastic fracture mechanics. ⁴⁷

Discrete scheme

MLS has no interpolation, so we use the penalty function as the natural boundary condition. Utilizing the theory of minimum potential energy, we get the Galerkin discrete form of the elastic body equilibrium governing equation as follows

$$\mathbf{K}\left(a\right)\tilde{\mathbf{u}}=\mathbf{F}\left(a\right)$$ (18)

where $\tilde{\mathbf{u}}$ is the generalized displacement

and

$$\mathbf{K}\left(a\right)=\left[\begin{array}{cccc}{\mathbf{K}}_{11}\left(a\right)\hfill & {\mathbf{K}}_{12}\left(a\right)\hfill & \cdots \hfill & {\mathbf{K}}_{1\left(2N\right)}\left(a\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{21}\left(a\right)}& {\mathbf{K}}_{22}\left(a\right)\hfill & \cdots \hfill & {\mathbf{K}}_{2\left(2N\right)}\left(a\right)\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{\left(2N\right)1}\left(a\right)}& {\mathbf{K}}_{\left(2N\right)2}\left(a\right)\hfill & \cdots \hfill & {\mathbf{K}}_{\left(2N\right)\left(2N\right)}\left(a\right)\hfill \end{array}\right]$$ (19)

is the stiffness matrix with

$${\mathbf{K}}_{\mathrm{IJ}}\left(a\right)=\underset{\mathrm{\Omega }}{\int }{\mathbf{B}}_I^{\mathrm{T}}\mathbf{D}\left(a\right){\mathbf{B}}_{J}d\mathrm{\Omega }+\alpha \left(a\right)\underset{{\mathrm{\Gamma }}_u}{\int }{\Phi }_I\left(\mathbf{x}\right)\mathbf{S}{\Phi }_J\left(\mathbf{x}\right)d\mathrm{\Gamma }$$ (20)

where

$${\mathbf{B}}_I=\left[\begin{array}{cc}{\Phi }_{I,x}\left(\mathbf{x}\right)\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\Phi }_{I,y}\left(\mathbf{x}\right)\hfill \\ \multicolumn{1}{c}{{\Phi }_{I,y}\left(\mathbf{x}\right)}& {\Phi }_{I,x}\left(\mathbf{x}\right)\hfill \end{array}\right]$$ (21)

where $\mathbf{D}(a)$ is the elastic matrix of the planer problem, ${\mathbf{B}}_I$ is the strain matrix, and N is the total number of nodal points in Ω

and

$$\mathbf{F}\left(a\right)=\left\{\begin{array}{c}{\mathbf{F}}_1\left(a\right)\hfill \\ \multicolumn{1}{c}{{\mathbf{F}}_2\left(a\right)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathbf{F}}_{2N}\left(a\right)}\end{array}\right\}$$ (22)

where

$$F_I\left(a\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }{\Phi }_I\left(x\right)b\left(a\right)d\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\Phi }_I\left(x\right)\overline{t}\left(a\right)d\mathrm{\Gamma }}+\alpha \left(a\right){\displaystyle \underset{{\mathrm{\Gamma }}_u}{\int }{\Phi }_I\left(x\right)Sud\mathrm{\Gamma }}$$ (23)

$$\mathbf{S}=\left[\begin{array}{cc}{s}_x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& s_y\hfill \end{array}\right]$$ (24)

$$s_i=\{\begin{array}{cc}1\hfill & \mathrm{if}{u}_i\mathrm{is\; prescribed\; on}{\mathrm{\Gamma }}_u\hfill \\ \multicolumn{1}{c}{0}& \mathrm{if}{u}_i\mathrm{is\; not\; prescribed\; on}{\mathrm{\Gamma }}_u\hfill \end{array}i=x,y$$ (25)

where $\overline{\mathbf{u}}$ and $\overline{\mathbf{t}}(a)$ are the displacement and surface force on the boundary, respectively; $\mathbf{b}(a)$ is the body force; the penalty function $\alpha (a)$ is generally ${10}^5-{10}^8$ times of the elastic module. To compute equations (20) and (23), we can use the background mesh algorithm to proceed the numerical integration.

While a varies in the a^I, $\mathbf{K}(a)$ and $\mathbf{F}(a)$ also vary in $\mathbf{K}(a^I)$ and $\mathbf{F}(a^I)$ , respectively. Utilizing the natural interval expansion of the function, we have

$$\mathbf{K}\left(a^I\right){\tilde{\mathbf{u}}}^I=\mathbf{F}\left(a^I\right)$$ (26)

Equations (20) and (23) can be written as

$$\mathbf{K}\left(a^I\right)=\underset{\mathrm{\Omega }}{\int }{\mathbf{B}}^{\mathrm{T}}\mathbf{D}\left(a^I\right)\mathbf{B}d\mathrm{\Omega }+\alpha \left(a\right)\underset{{\mathrm{\Gamma }}_u}{\int }{\Phi }_I\left(\mathbf{x}\right)\mathbf{S}{\Phi }_J\left(\mathbf{x}\right)d\mathrm{\Gamma }$$ (27)

$$F_I\left(a^I\right)={\displaystyle \underset{\mathrm{\Omega }}{\int }{\Phi }_I\left(x\right)b\left(a^I\right)d\mathrm{\Omega }}+{\displaystyle \underset{{\mathrm{\Gamma }}_t}{\int }{\Phi }_I\left(x\right)\overline{t}\left(a^I\right)d\mathrm{\Gamma }}+\alpha \left(a\right){\displaystyle \underset{{\mathrm{\Gamma }}_u}{\int }{\Phi }_I\left(x\right)Sud\mathrm{\Gamma }}$$ (28)

Equation (26) is exactly the static force interval element-free problem, which generally aims to determine the displacement set $\mathrm{\Gamma }$ using all possible combinations of the stiffness matrix $\mathbf{K}(a^I)$ and load $\mathbf{F}(a^I)$ in any given interval

$$\mathrm{\Gamma }=\{\tilde{\mathbf{u}}:\mathbf{K}\left(a\right)\tilde{\mathbf{u}}=\mathbf{F}\left(a\right),a\in a^I\}$$ (29)

The interval form of ${\tilde{\mathbf{u}}}^I$ is

$${\tilde{\mathbf{u}}}^I=[\underset{¯}{\tilde{\mathbf{u}}},\overline{\tilde{\mathbf{u}}}]$$ (30)

where

$$\underset{¯}{\tilde{\mathbf{u}}}=min\left\{{\mathbf{K}}^{-1}\left(a^I\right)\mathbf{F}\left(a^I\right)\right\}$$ (31)

$$\overline{\tilde{\mathbf{u}}}=max\left\{{\mathbf{K}}^{-1}\left(a^I\right)\mathbf{F}\left(a^I\right)\right\}$$ (32)

Since the uncertainty parameter a varies in small range, we have

$$\begin{array}{l}{K}^m=K\left(a^m\right),{K^{\prime }}_{,i}\left(a^m\right)=\frac{\partial K}{\partial a_i}{|}_{a=a_i},\\ F^m=F\left(a^m\right),{F^{\prime }}_{,i}\left(a^m\right)=\frac{\partial F}{\partial a_i}{|}_{a=a_i},\left(i=1,2,\mathrm{\dots },n\right)\end{array}$$ (33)

where ${\mathbf{K}}^m$ , ${\mathbf{F}}^m$ are the means of the whole stiffness matrix and the load array, respectively; $a^m$ is the mean of a, and n is the number of interval parameters.

Ignoring the high-order items, we have the Taylor expansions at $a=a^m$ of $\mathbf{K}(a^I)$ and $\mathbf{F}\in \mathbf{F}(a^I)$ that include the interval parameters

$$\mathbf{K}\left(a^I\right)={\mathbf{K}}^m+\sum _{i=1}^n\mathrm{\Delta }{a}_i^I·{\mathbf{K}\prime }_{,i}\left(a^m\right)={\mathbf{K}}^m+\mathrm{\Delta }{\mathbf{K}}^I$$ (34)

$$\mathbf{F}\left(a^I\right)={\mathbf{F}}^m+\sum _{i=1}^n\mathrm{\Delta }{a}_i^I·{\mathbf{F}\prime }_{,i}\left(a^m\right)={\mathbf{F}}^m+\mathrm{\Delta }{\mathbf{F}}^I$$ (35)

where $\mathrm{\Delta }{a}^I$ is the deviation of a.

From the interval element-free governing equation, we have

$$({\mathbf{K}}^m+\mathrm{\Delta }{\mathbf{K}}^I)({\tilde{\mathbf{u}}}^m+\mathrm{\Delta }{\tilde{\mathbf{u}}}^I)=({\mathbf{F}}^m+\mathrm{\Delta }{\mathbf{F}}^I)$$ (36)

Since the decomposed form of the interval is unique, we have

$${\tilde{\mathbf{u}}}^m={\left({\mathbf{K}}^m\right)}^{-1}{\mathbf{F}}^m$$ (37)

$$\mathrm{\Delta }{\tilde{\mathbf{u}}}^I={\left({\mathbf{K}}^m\right)}^{-1}\mathrm{\Delta }{\mathbf{F}}^I-{\left({\mathbf{K}}^m\right)}^{-1}\mathrm{\Delta }{\mathbf{K}}^I{\tilde{\mathbf{u}}}^m$$ (38)

Substituting $\mathrm{\Delta }{\mathbf{F}}^I$ and $\mathrm{\Delta }{\mathbf{K}}^I$ into equations (38), we have

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }{\tilde{\mathbf{u}}}^I={\left({\mathbf{K}}^m\right)}^{-1}\sum _{i=1}^n\mathrm{\Delta }{a}_i^I{\mathbf{F}\prime }_{,i}\left(a^m\right)}\\ \multicolumn{1}{c}{-{\left({\mathbf{K}}^m\right)}^{-1}\sum _{i=1}^n\mathrm{\Delta }{a}_i^I{\mathbf{K}\prime }_{,i}\left(a^m\right){\tilde{\mathbf{u}}}^m}\end{array}$$ (39)

$$\mathrm{\Delta }\tilde{\mathbf{u}}=\mathrm{\Delta }{a}_i|{\left({\mathbf{K}}^m\right)}^{-1}\sum _{i=1}^n{\mathbf{F}\prime }_{,i}\left(a^m\right)-{\left({\mathbf{K}}^m\right)}^{-1}\sum _{i=1}^n{\mathbf{K}\prime }_{,i}\left(a^m\right){\tilde{\mathbf{u}}}^m|$$ (40)

The upper and lower bounds of the structural static displacement are

$${\tilde{\mathbf{u}}}^I=[\underset{¯}{\tilde{\mathbf{u}}},\overline{\tilde{\mathbf{u}}}]=[{\tilde{\mathbf{u}}}^m-\mathrm{\Delta }\tilde{\mathbf{u}},{\tilde{\mathbf{u}}}^m+\mathrm{\Delta }\tilde{\mathbf{u}}]$$ (41)

where ${\tilde{\mathbf{u}}}^m$ and $\mathrm{\Delta }\tilde{\mathbf{u}}$ are the mean and deviation of the generalized displacement, respectively. ${\mathbf{\epsilon }}^I$ and ${\mathbf{\sigma }}^I$ can be determined in the same way. By the perturbation method of linear equation, it can be seen that the convergence condition of the perturbation method cannot be satisfied when variation range of the interval variable is large. In this case, the subinterval perturbation method is applicable. ⁴⁸

Because ${\tilde{\mathbf{u}}}^I$ is the generalized displacement, we have

$${\mathbf{u}}^I={\mathbf{\Phi }}^{\mathrm{T}}\left(\mathbf{x}\right){\tilde{\mathbf{u}}}^I$$ (42)

Interval J-integration theory

J-integration, proposed independently by Rice in 1968, was found to be the first translation integration of the energy momentum tensor. ⁴⁹ As shown on the two-dimensional cracked body in Figure 1, we choose an arbitrary smooth closed curve $\mathrm{\Gamma }$ , then from an arbitrary point on the surface below the crack, go counterclockwise along $\mathrm{\Gamma }$ to surround the crack tip, and stop at an arbitrary point on the upper surface. J-integration is defined as follows

$$J=\underset{\mathrm{\Gamma }}{\int }(\mathbf{W}\mathrm{dy}-\overrightarrow{\mathbf{T}}\frac{\partial \overrightarrow{\mathbf{u}}}{\partial x}\mathrm{d}s)$$ (43)

where $\mathbf{W}$ is the strain energy density of the plate; $\overrightarrow{\mathbf{T}}$ , $\overrightarrow{\mathbf{u}}$ , and $\mathrm{d}s$ are the stress component, displacement component, and elementary arc on $\mathrm{\Gamma }$ , respectively. Due to the conservativeness of J-integration, the value of integral is unrelated to the chosen path.

Figure 1.

J -integration curve surrounding the crack tip.

Ignoring the higher-order item and from the Taylor expansion at $a=a^m$ of $J(a^I)$ that includes the interval parameters, we have

$$J\left(a^I\right)=J^m+\sum _{i=1}^n\mathrm{\Delta }{a}_i^I·{J\prime }_{,i}\left(a^m\right)=J^m+\mathrm{\Delta }{J}^I$$ (44)

$$J{\prime }_{,i}\left(a^m\right)=\frac{\partial J^m}{\partial a_i}{|}_{a=a_i}$$ (45)

A rectangular closed curve was used to compute the J-integration (Figure 2). Due to the symmetry around the x-axis, we have

$$J^m=J_W^m-J_T^m$$ (46)

$$J_W^m=2[\underset{0}{\overset{c}{\int }}{W}_4^m\mathrm{dy}+\underset{c}{\overset{0}{\int }}{W}_6^m\mathrm{dy}]$$ (47)

$$J_T^m=2[\underset{0}{\overset{c}{\int }}{(-{\sigma }_x^m{\epsilon }_x^m-{\tau }_{\mathrm{xy}}^m\frac{\partial v^m}{\partial x})}_6\mathrm{ds}+\underset{\frac{d}{2}}{\overset{-\frac{d}{2}}{\int }}{({\tau }_{\mathrm{xy}}^m{\epsilon }_x^m+{\sigma }_y^m\frac{\partial v^m}{\partial x})}_5\mathrm{ds}+\underset{0}{\overset{c}{\int }}{({\sigma }_x^m{\epsilon }_x^m+{\tau }_{\mathrm{xy}}^m\frac{\partial v^m}{\partial x})}_4\mathrm{ds}]$$ (48)

$$\mathrm{\Delta }{J}^I=\sum \mathrm{\Delta }{a}_i^I{J\prime }_{,i}\left(a^m\right)$$ (49)

$$\mathrm{\Delta }J=\mathrm{\Delta }{a}_i|\sum {J\prime }_{,i}\left(a^m\right)|$$ (50)

$$J^I=[\underset{¯}{J},\overline{J}]=[J^m-\mathrm{\Delta }J,J^m+\mathrm{\Delta }J]$$ (51)

Figure 2.

Rectangular closed curve.

Under the situation of linear elasticity, when the passion ratio $\mu$ is not the interval variable, there is a simple relationship between $K_{\mathrm{I}}^I$ and $J^I$ integration

$$K_I^I=\sqrt{E^I{J}^I}\left(\mathrm{Plane\; stress}\right)$$ (52)

$$K_I^I=\sqrt{\frac{E^I{J}^I}{1-{\mu }^2}}\left(\mathrm{Plane\; strain}\right)$$ (53)

Numerical calculation

Bar under linearly distributed load

One end of the bar was fixed, the other was free, with the linear load P(x) = x applied. The bar parameters including length L = 1.0, cross-sectional area A = 1.0, and elastic modulus E^m = 1000, ΔE = 100 (Figure 3). Single Gaussian point and 11 equidistant nodes were used at scale = 2.1 (Figure 4). The example is only used to verify the correctness of IEFGM.

Figure 3.

Unit-length bar under linearly distributed load: geometric model.

Figure 4.

Node distribution at unit length bar (11 nodes).

Figure 5 shows the displacements at the upper boundary $\overline{u}$ and lower boundary $\underset{¯}{u}$ of the bar that had uncertainty parameters calculated by IEFGM and analytical solution. Figure 6 shows the stresses at the upper boundary $\overline{\sigma }$ and lower boundary $\underset{¯}{\sigma }$ of the bar that had uncertainty parameters calculated by IEFGM and analytical solution. Clearly, the medians of displacement and stress (u^m, σ^m) are exactly the displacement and stress at the median using the EFGM, which show that IEFGM is correct and effective.

Figure 5.

Displacement distribution interval of a one-dimensional bar.

Figure 6.

Stress distribution interval of a one-dimensional bar.

Bi-material cantilever beam

As shown in Figure 7, a bi-material cantilever beam was imposed at the upper part with even-distributed load $P=25\mathrm{MPa}$ . The beam was in size length $L=60\mathrm{mm}$ , height D = 12 mm, and unit thickness. The elastic moduli of the two materials were $E_1^m=3\times {10}^7\mathrm{MPa}$ , $\mathrm{\Delta }{E}_1=3\times {10}^6\mathrm{MPa}$ , $E_2^m=2\times {10}^7\mathrm{MPa}$ , $\mathrm{\Delta }{E}_2=2\times {10}^6\mathrm{MPa}$ , Poisson’s ratio ${\mu }_1={\mu }_2=0.3$ . Figure 8 shows the distribution of 16 × 7 regular nodes, and the interface nodes are shared by materials 1 and 2. The 15 × 6 background grids and 4 × 4 Gaussian nodes were used at the scale = 3.5.

Figure 7.

Cantilever beam under uniform load distribution.

Figure 8.

Node distribution.

Table 1 lists displacement and stress by IEFGM and the Monte Carlo method (MCM) with 10,000 samples along height and at x = 50 mm. At y = 0 mm, ${\mathit{u}}_1^I$ and ${\mathit{u}}_2^I$ of material 1 are equal to those of ${\mathit{u}}_1^I$ and ${\mathit{u}}_2^I$ of material 2, which accord with displacement continuity at the interface. As for material 1, ${\mathbf{\sigma }}_x^I=[1.6290,4.4712]\mathrm{MPa}$ and ${\mathbf{\sigma }}_y^I=[-13.679,-12.829]\mathrm{MPa}$ ; as for Material 2, ${\mathbf{\sigma }}_x^I=[-8.282,-6.274]\mathrm{MPa}$ , ${\mathbf{\sigma }}_y^I=[-13.679,-12.829]\mathrm{MPa}$ , which validate the bilateral positive stress ${\mathbf{\sigma }}_x^I$ is discontinuous and ${\mathbf{\sigma }}_y^I$ is continuous at the interface. These results satisfy the stress interface condition and validate the correctness and effectiveness of the IEFGM in resolving interface mechanical problems.

Table 1.

Displacement and stress along the height and at x = 50 mm.

	y (mm)	${\mathit{u}}_1^I$ (10−4mm)		${\mathit{u}}_2^I$ (10−3mm)		${\mathbf{\sigma }}_x^I$ (MPa)		${\mathbf{\sigma }}_y^I$ (MPa)
		IEFGM	MCM	IEFGM	MCM	IEFGM	MCM	IEFGM	MCM
Material 1	6	[12.504, 15.630]	[12.864, 15.179]	[8.472, 10.588]	[8.596, 10.432]	[46.212, 48.768]	[47.179, 47.796]	[−24.752, −24.716]	[−24.726, −24.733]
	4	[7.8770, 9.8470]	[8.0937, 9.5760]	[8.470, 10.586]	[8.601, 10.421]	[36.619, 37.579]	[36.976, 37.188]	[−22.883, −22.609]	[−22.848, −22.712]
	2	[2.9250, 4.4350]	[3.0076, 4.3096]	[8.468, 10.584]	[8.601, 10.417]	[17.943, 19.335]	[18.608, 18.658]	[−18.164, −17.466]	[−17.981, −17.627]
	0+	[−2.269, −0.717]	[−2.199, −0.739]	[8.467, 10.583]	[8.594, 10.424]	[1.6290, 4.4712]	[1.6722, 4.352]	[−13.679, −12.829]	[−13.320, −13.164]
Material 2	0−	[−2.269, −0.717]	[−2.199, −0.739]	[8.467, 10.583]	[8.594, 10.422]	[−8.282, −6.274]	[−8.0400, −6.4573]	[−13.679, −12.829]	[−13.242, −13.164]
	−2	[−7.498, −5.788]	[−7.258, −5.972]	[8.466, 10.582]	[8.593, 10.422]	[−20.854, −20.056]	[−20.694, −20.190]	[−6.281, −5.741]	[−6.103, −5.903]
	−4	[10.498, 13.122]	[10.790, 12.756]	[8.466, 10.582]	[8.593, 10.423]	[−35.755, −34.789]	[−35.718, −34.799]	[−1.941, −1.747]	[−1.886, −1.795]
	−6	[15.118, 18.898]	[15.613, 18.279]	[8.467, 10.583]	[8.591, 10.428]	[−44.289, −42.281]	[−43.723, −42.778]	[−0.228, −0.192]	[−0.221, −0.197]

IEFGM: interval element-free Galerkin method; MCM: Monte Carlo method.

A collinear double edged cracked finite square plate

Figure 9 shows a collinear double edged cracked finite square plate subjected to uniform tension stress P^I on both up and down boundaries. The crack length is a = 5 m and l = 20 m. Poisson’s ratio µ = 0.3. The elastic module E^I and the load P^I are the interval variables, with mean values P^m = 1.0 MPa, E^m = 1000 MPa, and deviations ΔP = 0.1 MPa, ΔE = 100 MPa. It is supposed to be a plane strain problem. Due to the bidirectional symmetry, we only used a quarter of the model in the up-left part. As shown in Figure 10, the symmetrical boundary conditions were applied on the right and bottom sides. Totally, 472 nodes were arranged on the plate, with the refined nodes on the crack tip. An area surrounded by Q₁Q₂Q₃Q₄ (length c = 4.0 m, width b = 2.0 m) was selected to solve the J -integration. In the integration, 8 × 8 and 4 × 4 Gauss points were used in the inner and other of the area, respectively. Complete expansion basis ³⁰ was used when the distance from the node to the crack tip is less than 0.1a, and otherwise, the quadratic basis was used. The penalty function is $\alpha ={10}^6{E}^m$ at $\mathrm{scale}=3.5$ .

Figure 9.

Double-edged cracked plate under mode I loading.

Figure 10.

Lay-outs of nodes, boundary conditions and the J-integration domain.

The interval of the stress intensity factor is $K_I^I=[4.3949,5.3660]\mathrm{MPa}\sqrt{\mathrm{m}}$ by IEFGM. The MCM with 10000 samples is introduced, and simulation result is $K_I^I=[4.1782,5.0881]\mathrm{MPa}\sqrt{\mathrm{m}}$ . When the deviations of the interval variables are ΔP = 0, ΔE = 0, we have $K_I^m=4.8804\mathrm{MPa}\sqrt{\mathrm{m}}$ with an error of only 1.57% compared to the analytical solution of the definite structure. Figure 11 shows that the radial stress interval $[{\underset{¯}{\sigma }}_{\mathrm{rr}},{\overline{\sigma }}_{\mathrm{rr}}]$ ( $\theta =0$ ) by IEFGM and MCM with 10,000 samples in front of the crack tip and the crack length change along with r/a, which is the ratio of the “distance from the given point to the crack tip” to the crack length. Clearly, the stress interval at r/a = 0.002 is [18.9394, 23.1329] MPa. The MCM with 10,000 samples is introduced, and simulation result is [17.6112, 20.3124] MPa. The curves of ${\underset{¯}{\sigma }}_{\mathrm{rr}}$ and ${\overline{\sigma }}_{\mathrm{rr}}$ show the singularity of $1/\sqrt{r}$ near the crack tip stress field. It can be seen that this method is correct and effective.

Figure 11.

Stress interval ahead of the crack tip.

Conclusion

We proposed the IEFGM after investigating interval mathematics and combining with the inner product space and element-free Galerkin method. The IEFGM only needs the node information, without requirement for element connectivity. The interval number decomposition was used to solve the interval equilibrium equations. The interval J-integration formula was deduced in detail. Several numerical calculations examples including a bar under uniform load, a bi-material cantilever beam subjected to uniform stress, and collinear double edged cracked plate subjected to uniform tensile stress, including uncertainty parameters all show that IEFGM is accurate and effective in solving uncertainty problems.