Direct methods for solving the Variation problem for multicriteria estimation of the bearing capacity of geomaterials Текст научной статьи по специальности «Математика»

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

UDC 539.3

DIRECT METHODS FOR SOLVING THE VARIATION PROBLEM

FOR MULTICRITERIA ESTIMATION OF THE BEARING CAPACITY

OF GEOMATERIALS

Igor A. BRIGADNOV

Saint-Petersburg Mining University, Saint-Petersburg, Russia

The article deals with direct methods for solving the variational problem in stresses for multicriteria estimation of the bearing capacity of a geomaterial sample in the current configuration, which can be both reference (unde-formed) and actual (deformed). The problem is to minimize the integral quadratic functional from the various stress components in the selected control subdomain on a set of stress fields statically balanced with external influences. For the simplest configurations of the sample, it is proposed to use the method of generalized Fourier series in Hilbert spaces. For complex configurations of a sample with stress concentrators, it is suggested to use finite element approximation with the subsequent minimization of a finite-dimensional quadratic function with linear constraints of equalities. A substantial numerical example is given for estimating the bearing capacity of a sample from a geoma-terial under pure compression.

Key words: geomaterial; variation problem in stresses; bearing capacity; multicriteria estimation; generalized Fourier series; finite-element approximation

How to cite this article: Brigadnov I.A. Direct Methods for Solving the Variation Problem for Multicriteria Estimation of the Bearing Capacity of Geomaterials. Journal of Mining Institute. 2018. Vol. 232, p. 368-374. DOI: 10.31897/PMI.2018.4.368

Introduction. Rocks and concrete are some of the basic construction materials, and therefore the assessment of their bearing capacity is a very urgent scientific and technical problem [3, 13]. In [1] the author proposed an original approach, based on a multicriteria estimation of the bearing capacity of a final sample from a geomaterial in a reference (undeformed) or actual (deformed) configuration. Within the framework of this approach, a variational problem for stresses in a given subdomain was set, where, depending on engineering concepts, the root-mean-square values of various stress components are estimated, then based on their combination the bearing capacity of the current body configuration relative to the specified external influences is assessed. The proposed approach can be attributed to nonlocal methods, in which the destruction of a geomaterial is considered as a physical process occurring near the stress concentrator [7], and not at the material point where the maximum stress value is reached. With the help of the new method, it is possible to construct comprehensive lower estimation impact of the external forces that guarantee the destruction of the final geomaterial sample, which is extremely important for the rock crushing processes [13].

Formulation of the problem. Let us consider a final geomaterial sample as a deformable solid that in the current configuration occupies a bounded Lipschitz domain QeR3. The body is in an equilibrium stress-strain state under the action of external stationary effects: a fixed displacement U is given on a non-empty section of a boundary r1 with a Lebesgue measure I T1! > 0, at Q the volume force with a density f is applied, and at the section of a boundary T2 = 5Q\GT1 a surface force is applied with a density g. Then, according to the principle of virtual power, the basic equation of statics is valid [2, 5]:

n(Q,T0, v) = A(Q, f, g, v), Vv eF (Q), (1)

where n(Q,T0,v) = JT0 --(V® v)TdQ and A(Q, f, g, v) = J f ■ v dQ + J g ■ v dy - powers of internal

Q Q G2

stress T0 and external forces (f, g) at permissible speeds of displacements v e F(Q) = {v: Q ^ R3; v(x) = 0, x e T1} respectively. Here and below, the dot denotes the scalar product of vectors, and two consecutive dots - the double scalar product (convolution); V = d/dx is the symbolic Hamiltonian operator; the symbol ® is a tensor (dyadic) product; the superscript T denotes the transpose operation [4]. In reference configuration T0 is the first (asymmetrical) tensor of

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

nominal stresses of Piola - Kirchhoff, and in actual configuration T0 is the symmetric tensor of true Cauchy stresses [2, 5, 9]. Note that the actual configuration of the sample can be obtained by solving, for example, the general problem of elastic-viscoplastic deformation [5].

When assessing the bearing capacity of a deformable solid with respect to specified external influences, it is necessary to take into account the features of the topology of the current configuration, since the most critical are the stress concentration zones (holes, grooves, cuts, etc.). Therefore, we single out a bounded Lipschitz subdomain ra e Q, containing a stress concentrator. The influence of the remaining part of the body on the selected (further control) subdomain is replaced by the following boundary conditions: a section of the subdomain boundary y1 e r1 is fixed, the rest of the boundary y2 = dra/y1 has a preset external force nT0, where n is the external in relation to a control domain normal to the boundary y2. Note that the fixed section of the boundary y1 may be absent.

Let us introduce into consideration a set of permissible velocities of displacements in ra:

V(ra) = {v e W 1,2(ra,R3): v(x) = 0, x e y1} and a quadratic functional on the Hilbert space of stresses L2(ra, M3) in the form of

K(ra,r) = -"^| fT--B-TdQ, (2)

2 ra J

I I ra

where B - a tensor of rank 4, depending only on the coordinate x. Here and below M3- space of real matrices 3 x 3.

Depending on the choice of the tensor B the functional of K(ra,T) has different physical meaning:

1/2

• if B = E is the first unit tensor of rank 4 [4], then K is the root-mean-square intensity of the stresses in ra;

1 1/2

• if b = E — I ® I, where I is the unit tensor of rank 2 [4], then K is - the root-mean-square

3

intensity of the shearing stresses in ra;

1 1/2

• if b = — I ® I, then K is the root-mean-square intensity of the hydrostatic pressure in ra;

• if b = 1 e —— I ® 11 is the tensor of material hardness, the inverse tensor of elasticity,

2^ 1 + v )

where ^ - shear modulus, and v - Poisson's ratio [2, 3, 5, 9], then K is the average specific internal energy of the deformed solid in ra .

In general, the choice of the control subdomain and the tensor are determined by engineering and technical considerations.

As before, in the control subdomain, the basic equation of the statics is valid n(ra,T 0,v) = A(ra, f, nT0, v) for any admissible displacement velocity field veV(ra), where the functional of n is determined in problem (1), and the functional of A has the form of

A(ra, f, n •T0, v) = J f •vdQ + J nT0 • v dy .

Let us set the following problem: among all statically equilibrium stress fields in the control subdomain ra we need to find one that minimizes the functional of K(ra,T), namely,

T* = arginf{K(p, T) : T e G}, (3)

G = {T e Z2(ra,M3) : n(ra,T,v) = A(ra, f,nT0,v), Vv e V(ra)} .

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

From the mathematical point of view, in this problem we search for the minimum of a quadratic functional on a linear affine manifold in a Hilbert space of stresses L2(ra, M3). Further we will consider only tensors B, satisfying the strict condition of Coleman-Noll, for which there exists a

constant a > 0 such that Q- • B - Q > a Q|2 for any Q e M3 and almost all x era [2, 9]. In this case

the functional K(ra,T) is strictly convex and coercive on L2(ra,M3) and therefore the existence and uniqueness of the solution of the problem (3) is guaranteed [6, 9].

Thus, in the control subdomain ra we consider the weakest stress fields, universally balanced in current configuration of the body Q, i.e. it is assumed that in ra the solid body resists

the external impact in the weakest way. In reality the body can balance much more intense ex-

1/2

ternal forces. If K exceeds the preset limit of stress, then the topology of the subdomain ra cannot provide the bearing capacity of the current body configuration in relation to set external factors, since any material can withstand only the final stress level (tangential, crimp, discontinuous, combined, etc.). Thus, we search for a lower estimate of the bearing capacity of the current configuration of the solid body with respect to the given external influences and the control subdomain. The proposed approach can be attributed to nonlocal methods for assessing the destruction of geomaterials [7].

In [1], the content of the proposed approach is analytically shown on the example of the problem of the spherically symmetric loading of a solid sphere by internal pressure in the reference configuration.

Problem-solving procedures. A peculiarity of the variational problem (3) is that the admissible set G contains the first derivatives of the sought-for functions and does not have the integral quadratic functional K(ra,T). Therefore, it is convenient to use direct methods to solve it. The most common one is the method of generalized Fourier series in Hilbert spaces [10]. However, it is applicable only for the simplest regions for which it is relatively easy to construct a countable basis. For complex regions, the finite element method [5, 8] is the most acceptable, which can be treated as a special case of finite Fourier series. In this case, as a result of finite-element approximation, the initial variational problem (3) is reduced to the simplest finite-dimensional problem of minimizing a quadratic function with linear constraints equalities, that can be easily solved in the MATLAB engineering calculations environment.

The method of generalized Fourier series. Let us choose a countable basis {^k(x)}keN, in Hilbert space L2(ra) and in Hilbert space W1'2(ra) a calculation basis {^m(x)}meN such that 9m(x) = 0 at y1, and in Hilbert space L2(y2) a calculation basis {yr(x)}reN. Then the following expansions in Fourier series are valid [10] (here and below we use the summation rule for repeating indices from 1 to œ):

a) for admissible stresses T(x) = Tk^k(x), where TkeM3 are real matrix coefficients;

6) for admissible displacement velocities v(x) = vm9m (x), where vm e R3 - real vector coefficients;

b) for set external forces f (x) = fk ^k (x) and g(x) = gryr (x), where fkeR3 and greR3 are the

corresponding real vector coefficients.

Because of the arbitrariness of the permissible velocities of displacements, i.e. vector coefficients vm (meN), the main problem (3) has the following form:

T }keN = arginf {iT • A • •Tj : T'cim = bm }, (4)

where the tensors of rank 4 Ay = Aj , satisfying the strict Coleman-Noll condition, are determined by the formula

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

vectors c'm are determined by the formula

cim =J<hV9m dQ and vectors bm are determined by the formula

bm =

m-f |фг фт dQ + gr JV rV

r T m

У2

Thus, the problem has been reduced to minimizing of a strongly convex and coercive quadratic form on a linear affine manifold or hyperplane in infinite-dimensional Euclidean (Hilbert) space, which ensures the existence and uniqueness of the solution [6].

To solve problem (3), we use the method of Lagrange multipliers [6]. Let us introduce the Lagrange function

L(Tk, r) = IT • -A1 • •T1 + Xm (T'c'm - bm),

the necessary condition for the stationarity of which has the form

r)L r)L

*h' •-T1 + Xm ® ckm = O, — = Tkckm - bm = 0, (5)

dTk

where O - null tensor of rank 2.

Each tensor of rank 4 A1, satisfying the strict Coleman-Noll condition, always has the inverse tensor of rank 4 B1, still satisfying the strict Coleman-Noll condition and such that B1 = B1 and B1 • -A1 = A1 • -B1 = E - the first unit tensor of rank 4 [2, 4]. Then from the first condition (5) we find that

Tk = -Bk • -Xm ® c'm = -r • Bk • c'm .

Substituting this expression into the second condition (5), we obtain an infinite system of linear algebraic equations (SLAE) with respect to Lagrange multipliers Xm. From the properties of fourth-rank tensors B1 the non-degeneracy of the second-rank tensorc'm • B1 • cJS follows [2, 4], therefore the obtained SLAE is uniquely solvable:

xm = -bs (c'm • Bj • cjs )-1.

Thus, the required matrix coefficients in the expansion of the stress field in the Fourier series with respect to the basis functions of the Hilbert space L (o) have the form

Tk = -bs(cim • Bj • cjsJ-1 • Bkr • crm . (6)

For example, if © = E and orthonormalize the basis functions in L2(o) by the Schmidt procedure [10], then A4 = В4 = 5^E, where ô^ is Kronecker symbol. In this case, the solution (6) takes a simpler form:

Tk = -(cim • cis)-1 bs ® ckm .

Unfortunately, the constructed analytic solution can be used only for the simplest subdomains o, поскольку since the choice of a countable basis in Hilbert spaces for arbitrary domains is an independent and intractable problem [6, 10].

Variational-difference method. To solve the original variational problem (3), we use the standard finite-element approximation with the subsequent numerical solution of the corresponding finite-dimensional problem.

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

As a rule, the current configuration is either a reference undeformed one, or it is determined numerically because of solving a boundary value problem for the whole domain Q within the framework of a model (for small or finite deformations, for elastic or inelastic medium, for displacements or stresses). In this case the standard finite element approximation is used. For dimension QeRn (n = 1, 2, 3) we create sets Qh and rh = dQh such that |Q \ Qh| ^ 0 u I r \ rh| ^ 0 for h ^ +0 with a regular pattern, where h - the largest step of triangulation and Qh is a set of simplest matched simplexes - segments, triangles or tetrahedra for spaces of dimension n = 1, 2 or 3 respectively [8].

The choice of control subdomian rah for further analysis is narrowed to selection the set of

finite elements in the vicinity of the stress concentrator with the transfer of boundary conditions to the subdomain boundary. Next, we can use a smaller triangulation in the control subdomain. Each triangulation is characterized by a set of nodes {xk}keI, where I - set of numbers of nods in . Let us note set of numbers of nodes at the sections of boundaries yh and y2 as I1 and I2 respectively.

For the admissible velocities of displacements, admissible stresses and external forces, we use the standard piecewise linear continuous approximation of the form

Vh ( x) = vm O m ( x), Th ( x) = Tk O k ( x), A ( x ) = fk Ok ( x ), gh ( x) = gr % ( x),

(7)

Where the lowest indices m e I \ I1, k e I and r e 12, and the values with a superscript denote the corresponding nodal values. Here we use standard continuous scalar functions of the following form [8]:

• linear function on each simplex O k : &h ^ R such that Ok (x1) = 5k., which finite support

consists of adjacent simplexes with common vertex xk;

• linear on each face of the corresponding simplex function : yh ^ R such that (x1) = 5h, which finite support consists of adjacent faces at y h with common vertex xr.

It is known [8], that linear space spanned by a linearly independent basis {Ok }keI, is a finite-dimensional subspace of the Hilbert space L2(®h), which is dense everywhere in L2(®h) with h ^ +0 in a regular pattern. Therefore, the approximations (7) can be considered as finite Fourier series of the corresponding fields, which allows all the preceding arguments to be completely transferred to the finite-dimensional case.

The problem (3) is approximated by the following finite-dimensional problem:

{Tk }kI = argmin {2T • -Aj • -T1 : Tcim = bm }, (8)

where indices i, j e I and m e I \ I1. Here tensors of rank 4 A11 = A11 are calculated according to formula

A =-L f BOOdQ.

ra I J

vectors Cm using formula

'1 i

I h\ rah

Cm = jOiVO m dQ

rah

and vectors bm by formula

bm = fk jOkOmdQ + gr f%r%mdy .

rah yh

372

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

Let us transform problem (8) to a form convenient for numerical solution. To do this, we collect all the matrix variables {Tk }keI in a global vector y with a dimension of Mn2, where M =| 11 - number of nods in &h. Further using the standard method we create global matrix G of square function with a dimension of Mn2 x Mn2, global matrix Q of linear constrictions with a dimension of Kn x Mn2, where K =| I \ Ij | - number of nods in ®h \ y\, and a global vector d of right right-hand sides of linear constraints with a dimension of Kn. Then the finite-dimensional problem takes the simple form

y* = arg min {\yTGy : Qy = d}, (9)

where index T is transposition. Note that the number of variables in problem (9) is always greater than the number of linear constraints, since M > K and n = 1, 2, 3.

For a fourth-rank tensor B in the functional (3), satisfying the strong Coleman-Noll condition, a regular finite-element approximation with subsequent standard assembly guarantees a strong positive definiteness of the matrix G of the quadratic function [8], therefore there always exists a unique solution of the finite-dimensional problem (9). Using the method of Lagrange multipliers, the solution is easily written in an explicit form, convenient for implementation in the MATLAB engineering calculations system,

y* = U (QU)-1 d, (10)

where U = G~lQT .

Computational experiments. The effectiveness of the variational-difference method for solving the main problem (3) can be easily demonstrated by evaluating the bearing capacity of a long square beam with pure compression along the longitudinal side surfaces [12]. In the context of the flat deformed state model, the cross-section of the bar (unit square) was divided into 2N2 equal triangular finite elements, where N is the number of equal fractions of the sides of the square. As a current, the deformed configuration of the beam section was found, which was found as a result of solving the following model problem of the linear theory of elasticity in the variational form for displacements during rigid kinematic deformation [12]:

1 1

u = arginf {J(u): u eV} , J(u) = JJ(I2(V® u) + X/2(V® u))dxdy, (11)

0 0

V = {u e W 12(Q,R2): ux(0,y) = 0,ux(1,y) = Ux,uy(0,y) = uy(1,y) = 0,y e (0,1)},

where u = (ux, uy) - flat displacement vector; I1 (V ® u) = x + uyy - intensity of volume deformation; I2(V® u) = (| (u2xx - uxxuyy + u2yy)+ 2u2xy)12 - intensity of shear deformation (the first and

second invariants of the Cauchy tensor of small deformations, respectively [2,9]). Here the dimen-sionless parameter A, = 3(1 - 2v)/(1 + v) depends on Poisson's ratio ve (0,05) [2, 9]; V - a set of ki-nematically possible displacements - such that the left side of the square is rigidly fixed, and the right is kinematically rigidly moved by Ux without vertical slip. As in [12], the finite-dimensional

approximation analogue of problem (11) was solved in the engineering calculation environment of MATLAB using the standard procedure fmincon for the number of fractions of the sides of the square N = 20 .

It is known that the destruction of a sample from a geomaterial under uniaxial compression (i.e., at Ux < 0) begins on a free lateral surface due to peeling [3, 11, 13]. Therefore, by the symmetry of the problem, the area of the upper free surface of the cross-section of a beam of three finite element layers was selected as the control subdomain, where the minimum root-mean-square values of vertical (flaking) stresses were estimated from formula (10).

For example, for the parameter X = 2, which corresponds to the Poisson's ratio v = 0.125 (granite, marble [3, 13]), and displacement Ux = -0.2 the minimum mean square vertical stresses T* found as a result of solving the problem (3) in the control subdomain, on average, 18 % less

^ Igor A Brigadnov

Direct Methods for Solving the Variation Problem

v 1

0.8

0.6

0.4

0,2

-0.2

i A A A / / / /

A A A A A V ¡3 / / V $ A A A K A > A /

A M A A A A A ■A / / / A A A A / A A /

A A / A / A A A / f / A A A A i A / A ?■■

A A / / A / A A / / / A A A A / A A A

A A / A & A / A / / / A A A A & A A: A A

/ 1 / A % A / A. I j / A A 1 /, A A A / A

/ rd A A A A A / £ / / A A / / 1 A A / A

i 0 / A A / A '£ / / / A A / % t A A / A

/ 0: A A A A / i / / / A A u A 1 A A. / A

t / A jg / / / 1 / / / A A A A / A A / A

A: / / $ A / / / / / / A A A A / A A / A

A A / A A / / / / / / / A A A 4 i A / A

A A / A A / A A

than the same stresses found as a result of solving the complete problem (11). Thus, if the value of K12 s greater than the ultimate brittle fracture stress for the studied material, the sample is guaranteed to collapse [11]. In the figure the red zone indicates areas where there may be a loss of global material stability.

In conclusion, we note once again that on the basis of the proposed multicriteria nonlocal analysis of the bearing capacity of the geomate-rial sample it is possible to construct absolute low estimates for external forces and impacts that are guaranteed to destroy the body, which is extremely important for the process of rock crushing [13].

Conclusions

0.2

0.4

0.6

0.8

The deformed configuration of the beam section. Zones of possible destruction of the sample due to flaking are marked with red.

1. A multicriterial estimation of the bearing

capacity of a geomaterial sample is reduced to the variational problem of minimizing the integral quadratic functional from the various stress components in the selected subdomain on a set of stress fields statically balanced with external influences.

2. A direct method for solving the variational problem on the basis of generalized Fourier series in Hilbert space is applicable only for the simplest domains for which it is relatively easy to construct a countable basis.

3. A direct method for solving the variational problem on the basis of finite-element approximation is applicable to complex dimensions with stress concentrators and is reduced to the simplest finite-dimensional problem of minimizing a quadratic function with linear constraints by the equalities.

0

REFERENCES

1. Brigadnov I.A. Multi-criteria assessment of the bearing capacity of geomaterials. Zapiski Gornogo instituta. 2016. Vol. 218, p. 289-295 (in Russian).

2. Lur'e A.I. Nonlinear theory of elasticity. Moscow: Nauka, 1980, p.512 (in Russian).

3. Nikolaevskii V.N. Geomechanics and Fluid Dynamics. Moscow: Nedra, 1996, p. 447 (in Russian).

4. Pal'mov V.A. Elements of tensor algebra and tensor analysis. St. Petersburg: Izd-vo Politekh. un-ta, 2008, p. 109 (in Russian).

5. Pozdeev A.A., Trusov P.V., Nyashin Yu.I. Large elastic-plastic deformations. Moscow: Nauka, 1986, p. 232 (in Russian).

6. Sea Zh. Optimization. Theory and algorithms. Moscow: Mir, 1973, p. 244 (in Russian).

7. Suknev S.V. Application of nonlocal and gradient criteria for the evaluation of fracture of geomaterials in zones of tension stress concentration. Fizicheskaya mezomekhanika. 2001. Vol. 14(2), p. 67-75 (in Russian).

8. S'yarle F. The finite element method for elliptic problems. Moscow: Mir, 1980, p. 512 (in Russian).

9. S'yarle F. Mathematical theory of elasticity. Moscow: Mir, 1992, p. 472.

10. Trenogin V.A. Functional Analysis. Moscow: Nauka, 1980, p. 496.

11. Cherepanov G.P. Mechanics of brittle failure. Moscow: Nauka, 1974, p. 640.

12. Brigadnov I.A. Regularization of non-convex strain energy function for non-monotonic stress-strain relation in the Hencky elastic-plastic model. ActaMechanica. 2015. Vol. 226. Iss. 8, p.2681-2691.

13. Verruijt A. Computational geomechanics. Dordrecht: Springer Science+Business Media, B.V., 1995, p.384.

Author Igor A. Brigadnov, Doctor of Physics and Mathematics, Professor, [email protected] (Saint-Petersburg Mining University, Saint-Petersburg, Russia).

The paper was received on 3 May, 2017.

The paper was accepted for publication on 2 March, 2018.