CC BY
7
DOI: 10.17516/1997-1397-2021-14-1-21-27 УДК 539.37
Numerical Modelling of Slow Motion of a Granular Medium
Olga I. Kuzovatova*
Siberian Federal University Krasnoyarsk, Russian Federation
Received 10.08.2020, received in revised form 16.09.2020, accepted 20.10.2020 Abstract. The aim of this work is to find approximate analytic solution of the problem of granular medium motion in a convergent channel, to develop computational algorithm based on the finite element method, and to carry out numerical calculations of the problem.
Keywords: variational inequality, materials with different strengths, strains localization. Citation: O.I. Kuzovatova, Numerical Modelling of Slow Motion of a Granular Medium, J. Sib. Fed. Univ. Math. Phys., 2021, 14(1), 21-27. DOI: 10.17516/1997-1397-2021-14-1-21-27.
Introduction
Many natural and artificial materials have different strengths in tension and compression, such are for example geomaterials constituting Earth's crust that include rock or granular materials, dry or saturated soil etc.
One of the fundamental problems of geomechanics is the problem of geomaterial motion in a convergent channel. Analysis of a motion of fragmented rock or a granular medium in a convergent channel is of prime interest for many mining technological processes (caving, drawing the stopes, motion in orepasses) and for grain storage and processing (storage bunkers or grain tank emptying). Similar processes take place in natural conditions as well, e.g. rock or soil displacement around mine shafts or grooves, shift trough formation above worked out areas or karst caverns and so on. The approximate (engineering) solution of this problem and the results of real experiments can be found in [1, 2].
The aim of this work is to find approximate analytic solution of the problem of loose medium motion in a convergent channel on the base of a model that takes into account different strengths of materials, to develop computational algorithm based on the finite element method, and to carry out numerical calculations of the problem.
1. Mathematical model
For description of the stress-strain state of a granular medium as a material with different strengths in tension and compression we shall use a model of a medium with plastic connections [3, 4]. Under the action of compressive or tensile stresses which are less than the cohesion coefficient (the yield point of plastic connections) such a material is not deformed. Until the yield point is reached the deformation follows the law of linear hardening. The rheological scheme of the model is given on Fig. 1.
*oik17@yandex.ru https://orcid.org/0000-0003-3635-4807 © Siberian Federal University. All rights reserved
"2
c
0 \
K
Fig. 1. Rheological scheme
Fig. 2. Deformation C and stress K cones
According to this scheme, there is an additive decomposition a = ac + a0 + ae, where a is the total strain tensor, ac is the rigid contact component, a0 is the cohesion tensor, ae = E : e is the elastic component, e is the deformation tensor, and E is the symmetric positively defined elastic modulus tensor. The tensor ac satisfies the variational inequality
oc : (e — e) < 0, e,ë G C,
(1)
where C is the cone of admissible deformations (Fig. 2) C = (e) < 6 (e)}, k is the dilatancy parameter, y (e) = %/2e' : e' is the intensity of shear, e' = e — 65/3 is the deviator of deformations, 6 (e) = e : 5 is the volume deformation, and 5 is the Kronecker symbol. In this notation, the inequality (1) take the form
(E : e — a + a0) : (e — e) > 0, e,e G C.
By definition of a projection, this means that
e = [E-1 (a — a0)],
where n is the projection operator onto the cone C with respect to the norm \e\ = %/e : E : e.
Consider an element of a construction from a material with different strengths filling a planar domain Q with the boundary dQ = r that consists of two non-intersecting parts T„ and . On the first part displacements are absent and on the second part the distributed load p is given. There hold equilibrium equations in variational form and boundary conditions:
/ (V • o + f ) (U - u) dil = 0, JQ
u = U = 0 on r„, o • n = p on rCT.
(2) (3)
The problem (2)-(3) reduces to the problem of finding the minimum min J (u) = J (u), where
u£Uc
J (u)
In (
-e:E:e + e: o —f•u Q \ 2
)da — I
p • udr,
Uc = {u G H1 (Q) | u\r
e (u) G C}.
Since C is the cone with the apex at the origin, we have J (u) = min min J (Xu)
ùeuc \>o
J (\u) = J ^e (u) : E : e (u) + Xe (u) : o" — Xf • uj dQ — XJ p • udr.
^ dQ — X J
0
2
A direct computation of minJ (Xu) shows that
\>0
A = — Therefore,
Jo (e (u) : a0 — f ■ u) dQ — fr p ■ udF
fn e (u) : E : e (u) dQ
where 2+
{0:
Z > 0:
z < 0.
J (u) = — max
ueuc u=0
Jo (e (u) : a0 — f ■ w) dQ — fr p ■ udF
2 fn e (u) : E : e (u) dQ
From this it can be shown that the displacement fields vanishes identically if and only if
(e (U) : a0 - f • U) dQ - p • UdF < 0, VU € UC. (4)
JQ JTa
A load (f,p) is called safe if u = 0. Let p = 0, f = m • f0, where m is the loading parameter. Then it follows from (4) the the load is safe for m varying from zero to the limit value (safety factor)
=+= ./o
m = mm o
fo e(u) : a0dQ
eu [Jo f ■ MQ] + ■
(5)
u=0
The proved statement is a formulation of a kinematic theorem on limiting equilibrium from plasticity theory [5].
2. Analytic solution
Consider as an example the problem of planar gravity flow of a granular medium in a convergent channel. Assume that a > ¡3 and consider two cases given on Figs. 3 and 4.
Fig. 3. Case 1
Fig. 4. Case 2
We compute safety factors mi and m2 and choose the least one, which is to be the safety factor m* for this problem.
For Case 1 the admissible displacement field u = (u1, U2) describes the strain localization for a simple shear with dilatancy in a narrow linear zone of thickness h inclined at an angle y. In the Cartesian coordinates related to this zone
ui = — u0 cos (a — p) : «2 = — u0 sin (a — p).
u0
Y0 = ~T cos (a — p). h
u0
e0 = — sin (a — p). h
+
2
+
We compute the integrals in (5).
i e (U) : a0dQ = e0a0S, JQ
in
r0
The "separating" triangular domain moves as a whole, hence
here a0 = ts/k, ts is the yield point, S = hl, l = a (cos y + sin y ctg (a — y)).
/ f0 ■ uodtt = f0 ■ SA, n
here f0 = pgu0 sin a, Sa = 2Hl, H = a sin y. Therefore the safety factor mi is equal to
2ts . sin (a — y)
m1 = -min-. (6)
Kpga ueuc sin a sin y
u=0
The condition u G UC takes the form y0 ^ ve0, where v = \Jl/K2 — 4/3. Then
S0, , 1 . , N 1
— = tg (a — y) ^ — or sin (a — y) ^
Y0 ^ v ^ sjv2 + 1'
The expression (6) attains its minimum for equalities in the latter ones. v sin a — cos a
Computing sin y =--—, we obtain the formula for the safety factor
V v2 + 1
2ts 1 1
mi
Kpga sin a (v sin a — cos a)
Studying Case 2 and performing similar computations we arrive at the following formula for the safety factor
2ts . sin (ft — V) . 2ts 1 1
m2 =- min ————— and m2 —
Kpga u g u s^n ft sin V Kpga sin ft (v sin ft — cos ft)
U = 0
Analysis shows that the inequality m1 < m2 holds for ft < a. Thus, the localization zone is inclined at the angle
y = a — arcsin —; (7)
and corresponds to Fig. 3.
It should be noted that in the engineering solution obtained in [1] the localization zones were given linear due to mechanical considerations. It turns out that their directions can be uniquely determined from the variational principle (5). Note that solutions with curvilinear localization zones were earlier considered in [7].
3. Numerical modelling
The computational algorithm is based on a finite-element approximation of the model. It reduces the problem of determining a displacement field in a material with different strengths to the solution of a number of static problems of the linear elasticity theory with initial stresses.
The algorithm does not employ theorems on a limit load estimate. Using it, the limit loads can determined only approximately as the loads such that exceeding them results in intensive deformation of a material [7].
As an initial approximation we take a solution of an elastic problem. The idea behind the algorithm is to replace the defining equations
a = E : e--- n (E : e — ao),
1 + X
where X = const ^ 0, by the iterative formula (n = 1,2, 3,...)
an = E : en--Kr n (E : en-i — ao) .
1 + X v 0/
Here n is the projection operator to the cone K with respect to the norm |a|0 = Va : E-1 : a, which acts as follows (k is the volume compression modulus and n is the shear modulus), [4]:
• if t (s) ^ Kp (s), then a = s;
• if t (s) > Kp (s) and fip (s) + KkT (s) ^ 0, then a = 0;
• if t (s) > Kp (s) and fip (s) + KkT (s) > 0, then
kp (a)t . , s n /NT , s № (s) + KkT (s)
a = ^r[s + p(s) • — p(a) • ^ p (a) =--.
t (s) n + K2k
The iteration steps are performed until the norm of the difference of two approximate solutions at neighboring steps becomes less than a given accuracy. /
The internal friction parameter varies in the interval 0 < k < %/3/2, calculations are carried out for k = 0.3. The regularization parameter X in computations was taken equal to 0.0001. It turns out that its further decreasing results in a low rate of convergence of the iterative algorithm. The accuracy is equal to e = 0.0001. Numerical experiments are performed using a program written in Borland C++.
4. Results of numerical experiments
In this section we consider a graphical interpretation of results of numerical experiments for the problem of planar gravity flow of a granular medium in a convergent channel (Fig. 3).
Denote
2ts 1
fcr = m*Pg =--.--:-r. (8)
Ka sin a (v sin a — cos a)
If f < fcr then the localization of deformations is absent; the medium is too light to move through the channel. If f > fcr then the medium moves, and in this case a sliding zone must be present.
Fig. 5 (a) depicts the cross-section of a symmetric convergent channel with the following parameters: a = 3 = 63.43°, a = 2 m. The finite-element grid presented in Fig. 5 (b) consists of 859 nodes and 1596 elements, the side of a triangle is 0.18 m with area 16 m2.
By formula (8) we compute the value fcr = 0.08ts. By formula (7) the most probable angle of departure of linear zone of the strain localization is y = 45.70°.
Fig. 6 (a) shows the shear intensity field obtained in the framework of the classical elasticity theory, while Fig. 6 (b) on the base of the model of a material with different strengths with the
a b
Fig. 6. The shear intensity field
internal friction parameter k = 0.3; a black line indicates the direction of the line of deformation inclined at the angle y = 45.70°.
Fig. 7 (a) shows a cross-section of an asymmetric convergent channel with the parameters a = 76.2°, 3 = 35°, a = 2 m. The finite-element grid presented in Fig. 7 (b) includes 1035 nodes and 1922 elements, the side of a triangle is 0.18 m with area 30.92 m2.
a b
Fig. 7. Asymmetric domain
By formula (8) we have fcr = 0.06ts. By formula (7) the most probable angle of departure of linear zone of the strain localization is y = 58.49°.
Fig. 8 (a) shows the shear intensity field obtained in the framework of the classical elasticity theory, while Fig. 8 (b) on the base of the model of a material with different strengths with the internal friction parameter K = 0.3; a black line indicates the direction of the line of deformation localization inclined at the angle y = 58.49°.
Thus, numerical experiments for the problem of flow of a granular medium in a convergent channel support the presence of linear zones of the strain localization inclined at an angle y given by formula (8) and correspond to Fig. 3.
Fig. 8. The shear intensity field
The work is carried out with the support of the Russian Foundation for Basic Research (grant 18-31-00019).
References
[1] S.V.Lavrikov, A.F.Revuzhenko, Localized flows of friable material in radial channels, Soviet Mining Science, 26(1990), 1-6.
[2] A.F.Revuzhenko, S.B.Stazhevskii, E.I.Shemyakin, Asymmetry of plastic flow in converging axisymmetric channels, Dokl. Akad. Nauk SSSR, 246(1979), no. 3. 572-574 (in Russian).
[3] V.P.Myasnikov, V.M.Sadovskii, Variational principles of the theory of the limiting equilibrium of media with different strengths, Journal of Applied Mathematics and Mechanics, 68(2004), no. 3, 437-446.
[4] O.Sadovskaya, V.Sadovskii, Mathematical modeling in mechanics of granular materials, Heidelberg, Springer, Berlin-Heidelberg, 2012.
[5] A.A.Gvozdev, Analysis of bearing capacity of structures by Limit State Method. Part 1. The essence of the method and its foundations, Gosstroyisdat, Moscow, 1949 (in Russian).
[6] D.C.Drucker, W.Prager, Soil mechanics and plastic analysis or limit design, Quarterly of Applied Mathematics, 10(1952), no. 2, 157-165.
[7] O.I.Kuzovatova, V.M.Sadovskii, Modeling of a Deformation Localization in a Medium with Different Strengths, J. Sib. Fed. Univ. Math. Phys., 1(2008), no. 3. 272-283 (in Russian).
Численное моделирование медленного движения сыпучей среды
Ольга И. Кузоватова
Сибирский федеральный университет Красноярск, Российская Федерация
Аннотация. Цель данной работы — найти приближенное аналитическое решение задачи о движении сыпучей среды в сходящемся канала для разработки вычислительного алгоритма на основе метода конечных элементов и проведения численных расчетов задачи.
Ключевые слова: вариационное неравенство, разнопрочная среда, локализация деформаций.
