Modified finite element analysis for exterior boundary problems in infinite medium

Modified algorithm of the 3d finite element analysis developed for solution of external boundary problem is considered. The method in question is based on incorporating the FEM and Somigliana’s integral formula. Efficiency of software implementations of the algorithm has been tested. A stress-strain analysis of inhomogeneous medium with a cavity has been carried out to display the approach.


Introduction
The underground non-reinforced caverns subjected to internal pressure or any other loads are often used in construction and mining as the gas or oil storages or for the waste dumping (see Fig. 1, a).Large-sized cavern is often located in an inhomogeneous medium because of complicated structure of the rock formation.Analysis of stability of such caverns should be made with account for specifics of the layered rock structure [1][2][3][4][5].Also for the land surface subsidence prediction and risk analysis of buildings it is necessary to estimate deformation of the rock mass caused by the presence of underground karst caverns (see Fig. 1, b).Consequently, the ground subsidence and horizontal displacement under building foundation are the important topics for the analysis [6][7][8][9][10][11].
a Corresponding author : chernat0000@mail.ruIn many cases of practical importance the analysis of structures considering their interaction with infinite foundation or the surrounding elastic medium is based on numerical solution of the threedimensional external boundary problem [12][13][14].Application of algorithms of the finite element method (FEM) or the boundary element method (BEM) reveals some peculiar properties of solution in infinite domain [15][16][17].The FEM analysis is performed assuming a reduced domain of finite size, or by application of special "infinite" elements [18,19].The problems inevitably arise in using this approach, such as validation of the reduced domain size and selection of boundary conditions on the external boundary in the first case, or complications encountered in solution of the system of equations in the second case.When the BEM approach application is considered the analysis of interaction of a structure and of inhomogeneous medium occurs inconvenient, as well as the case of domains with several cavities, meanwhile these problems are not rare in practice.
These specifics of the boundary problem in the case of infinite domain cause some computational complications and increase the computing cost [20,21].Various known algorithms offer to avoid such problems at the expense of combining different methods and equations [22][23][24][25][26][27].In the following, the algorithm based on combining the FEM facilities and Somigliana's integral formula, combined method (CM), is considered and its efficiency displayed [28][29][30].

Somigliana's integral formula
The Somigliana's integral formula allows determining components of the displacement vector T at arbitrary point [ in a space : bounded by a surface S (see Fig. 2); in the absence of the volume forces:

Algorithm of combined method (CM)
In the following, the algorithm of the 3d analysis for external boundary problem by combined method (CM) is presented.Let the displacement vector u([) should found in elastic space : bounded by a surface S, satisfying at [: the Navier's equilibrium equations and boundary conditions.The boundary conditions are where n i ([) are components of the unit outer normal at the surface S, V ij ([) are the components of stress tensor, p i ([) and u iS ([) are given components of force and displacement vectors.
A sub-space : 0 bounded by a surface S 0 is selected in the semi-infinite space : (see Fig. 3).The iteration procedure can be performed as follows.The first step in process is establishing of boundary condition u i (0) at the surface S 0 (assumed usually u i (0) = 0).Then the boundary problem in the subspace : 0 can be solved by FEM to complete boundary condition (Eq.2) on the surface S.This allows using the Somigliana's integral formula (Eq. 1) to determine the first approximation, u i (1) .Generally, having (k-1)-th approximation the k-th approximation u i (k) ([), [ S 0 can be found by applying the Somigliana's integral formula (Eq.1), using u i (k-1) ([), [ S 1 and t i (k-1) ([), [ S 2 determined by the FEM analysis of boundary problem in the sub-space : 0 with boundary conditions (Eq.2,3).The process can be terminated when the difference between the successive approximations u i The following approach is used for integrating the first term in Eq. 1 which contains t i (K), KS.The nodal forces P i in the nodes [ l of a finite element 'S in the right-hand side of the FEM equations can be treated as some integral characteristics of surface forces t i (K): An approximation is assumed: Strictly, this formula (Eq.6) is not a quadrature formula but it provides reasonable accuracy as an alternative of numerical differentiation in the FEM solution.
The software developed based on implementation of this algorithm has been tested in a series of the model problems.These were external boundary problems with self-balanced and non-selfbalanced external loads applied to a sphere-shaped concavity, such as uniform pressure, two selfbalanced point forces and one point force applied symmetrically or non-symmetrically.These problems were also analyzed by the FEM for the purpose of comparing the results.The FE meshes were different by the fineness and the number of finite element rows.In general, the accuracy of the CM results even at the number of elements far less that in the FEM modeling was higher than for the FEM.Apart from that, the displacements obtained by the CM for the case of non-self-balanced external loads include the rigid-body components unlike those determined by the FEM.This feature of the CM may be considered as providing better accuracy of results in such cases.The algorithm of combined method was modified for the case of inhomogeneous medium on the assumption that the inhomogeneous medium may have a substantial influence on the strain-stress state of the structure within a limited surface L. The remaining space : L can be considered as a homogeneous one.A surface S 0 is selected in the space : L (see Fig. 4) and a closed sub-space : 0 bounded by surfaces S and S 0 is considered.The first step in the process is establishing of boundary condition u i (0) at the surface S 0 (usually u i (0) = 0).Then the boundary problem in the sub-space : 0 can be solved by FEM to define boundary condition (Eq.2) on the surface L. The domain of integration in Eq. 1 is the boundary of homogeneous region, consequently it will be the surface L instead of the surface S. Then the boundary problem in the sub-space : 0 is solved by FEM to define the second approximation u i (2) on the surface L and so far.It is a double-boundary algorithm of combined method (DCM).
The values of stress and strain in Eq. 1 in the modified algorithm in question are calculated by the FEM approximately, while in the original algorithm part of them is defined by Eq. 2. By this reason designing of the FE mesh has to be more careful to decrease the error.
The developed software based on implementation of the DCM algorithm has been tested in stressstrain analysis of sphere-shaped cavity surrounded by a spherical layer with physical properties different from those in the surrounding space under the internal uniform pressure.The dimension of the analyzed domain for the CM analysis was 3 times less than that for the FEM analysis.The accuracy of the CM and the FEM results was comparable.The both methods resulted in partly underestimated stresses and strains compared to the results of exact solution due to insufficient dimensions of the considered domain.

Strain analysis of rock massif with a cavity
For example selected is an underground closed cylindrical cavern subjected to internal pressure on the boundary of two rock layers; a complete cohesion of the layers on the boundary was assumed (see Fig. 5).The layers' properties are characterized by the Young's modules, E 1 and E 2 , with E 2 varied from E 1 to 10 E 1 .The variety may be caused by variability of rock crumbling even in a homogeneous massif.The influence of relationship of Young's modules E 1 e E 2 on the stress-strain state of massif was evaluated.Calculations were carried out by applying the DCM as well as FEM to verify results.The modeled area for the DCM was two times smaller than that for FEM.The DCM FE model had two layers where the first one : L consisted of two different materials and the second one (: 0 _: L ) had assumed average Young's modulus E 3 .The FEM model had four layers consisting of two different materials.
The intrinsic stress state of massif caused by the weight of upper rock was accounted for by the method of deleting loads.For this purpose the FEM analysis of massif without cavern under the action of vertical and horizontal mass forces corresponding to depth of location was performed.By this way the loads on the surface S to be deleted were determined as nodal forces and then equivalent forces opposite in sign to above deleted forces were used to define the boundary conditions.
Safety factors K V = V u /V were calculated in the inner layers of finite elements around the cavern.There V u * = 0.75V u is the weakened rock compressive strength diminished with respect to the standard rock compressive strength V u .As might be expected, the tangential stress V T was the highest at the side surface of cavern (KV T = 1.09 y 1.43) as well as the vertical stress V z at the upper and the bottom surfaces (KV z = 1.15 y 1.71 for the upper surface and KV z = 1.05 y 1.15 for the bottom one).Values of the safety factor for the sound layer with modulus E 2 were smaller than those for the weak one; they decreased accordingly the ratio E 1 e E 2 .Therefore, it is the reasonable conclusion that the volumes where KV 1.2 should be reinforced.
Fig. 6, 7 illustrate displacements u z at the points A, C and u U at the point B compared with corresponding displacements in homogeneous massif with Young's modulus E 1 for point A or with modulus E 2 for point C as functions of the ratio E 1 e E 2 .Fig. 4 shows that displacement u zC of the point C in the weak rock volume is almost unaffected by the increasing modulus E 2 , while displacement u zA of the point A in the sound rock volume decreases irrespective the changes of modulus E 2 .The radial displacement u UB of the point B on the contact surface is smaller than that in homogeneous massif with Young's modulus E 1 .The displacement u UB is larger than that in homogeneous massif with Young's modulus E 2 if E 1 e E 2 !0.3 because of the weak rock dominant influence.In case of E 1 e E 2 0.3 displacement u UB decreases because of the strong rock influence.
As is also evident from the above (see Fig. 6, 7) all the results of the FEM and DCM analyses are mostly the same.So far, the stress-strain state of the massif can be estimated by applying the DCM.Besides, the considered domain and the number of finite elements necessary for the DCM analysis were far smaller than those for the FEM analysis.

Figure1.7KLV
Figure1.(a) Underground oil storage in a cavern (b) Underground karst caverns (i.e.solution cavities caused by ablation with underground water) where [:, KS, i, j = 1, 2, 3 are indexes of Cartesian coordinate axes x i , t i (K) is a component of the actual force vector and u i (K) is a component of the actual displacement vector at the point K at the surface S; G ij (K,[) and F ij (K,[) are the components of vectors of force and displacement caused by the unit force acting in direction x j at the point [ (the fundamental solution of Navier's equilibrium equations).

Figure 2 .Figure 3 .
Figure 2. Exterior domain bounded by a surface S Figure 3. Selection of analyzed domain in a space :

Figure 4 .
Figure 4. Scheme of the inhomogeneous space selection Figure 5. Underground cavern for analysis in inhomogeneous rock massif

Figure 6 . 7 .
Figure 6.The vertical displacements of the upper and bottom Figure 7.The radial displacement of point surfaces of the cavern on the cohesion surface