A comparison of discrete element simulations

合集下载

各向异性对软土力学特性影响的离散元模拟

DOI ：10.16030/ki.issn.1000-3665.202006027各向异性对软土力学特性影响的离散元模拟赵　洲1，宋　晶1,2,3，刘锐鸿1 ，杨守颖1 ，李志杰1（1. 中山大学地球科学与工程学院，广东广州　510275；2. 广东省地球动力作用与地质灾害重点实验室，广东广州　510275；3. 广东省地质过程与矿产资源探查重点实验室，广东广州　510275）摘要：软土预压工程中，初始和诱发各向异性对软土力学性质的影响十分显著，而现有研究缺乏对初始和诱发各向异性的统一研究方法。

采用离散单元法，以颗粒长宽比作为定量评价指标，构建真实形态的颗粒模型，生成5组不同沉积角的初始各向异性试样，并进行竖直和水平两方向加载的双轴模拟实验，研究了初始各向异性和诱发各向异性对软土力学特性影响；在细观层面，以颗粒为对象研究了颗粒接触形式和转动角度的变化规律，以接触为对象研究了配位数和接触法向各向异性的发展趋势，在此基础上探究抗剪强度指标与各向异性关系。

结果表明：初始和诱发各向异性共同影响试样力学性质，当加载方向和软土沉积方向垂直时，土体有最大的峰值强度。

颗粒接触形式中面面接触的比例随加载的进行逐渐增大，并影响着试样初始模量和抗剪强度，配位数和接触法向各向异性受颗粒接触形式的影响有不同的演化规律，并在加载后期趋于稳定；同时，初始各向异性试样相较各向同性试样有更大的黏聚力，诱发各向异性主要影响试样内摩擦角，进而影响试样抗剪强度。

关键词：软土；各向异性；峰值应力；抗剪强度指标；细宏观性质中图分类号：TU411.3 文献标志码：A 文章编号：1000-3665（2021）02-0070-08Discrete element simulation of the influence of anisotropy on themechanical properties of soft soilZHAO Zhou 1，SONG Jing1,2,3，LIU Ruihong 1 ，YANG Shouying 1 ，LI Zhijie1（1. School of Earth Sciences and Engineering , Sun Yat-Sen University , Guangzhou , Guangdong 510275, China ；2. Guangdong Provincial Key Lab of Geodynamics and Geohazards , Guangzhou , Guangdong 510275, China ；3. Guangdong Provincial Key Laboratory of Geological Processes and Mineral Resource Exploration ,Guangzhou , Guangdong 510275, China ）Abstract ：The initial and induced anisotropy has a significant effect on the mechanical properties of soft soil in preloading engineering. However, there is a lack of unified research methods for the initial and induced anisotropy. Discrete element method is adopted in this study, and the length-width ratio of particles is used as the quantitative evaluation index. Five types of initial anisotropy samples with different deposition angles are generated. The effects of initial anisotropy and induced anisotropy on the mechanical properties of soft soil are studied by vertical and horizontal loading. At the micro level, the contact form and rotation angle of particles are examined from the point of view of particles, and the development trend of coordination number and contact normal to anisotropy is studied from the point of view of contact. The relationship between shear strength index收稿日期：2020-06-12；修订日期：2020-08-04基金项目：国家自然科学基金项目（41877228；41402239；41877229）；广东省自然科学基金项目（2019A1515010554）；广州市科技计划项目（201904010136）第一作者：赵洲（1995-），男，硕士研究生，主要从事软土微观结构及数值模拟。

the mathematical theory of finite element method

the mathematical theory of finiteelement methodThe mathematical theory of the finite element method (FEM) is a branch of numerical analysis that provides a framework for approximating solutions to partial differential equations (PDEs) using discretization techniques. The finite element method is widely used in engineering and scientific disciplines to simulate and analyze physical phenomena.At the core of the FEM is the concept of dividing a domain into a finite number of elements, which are connected at nodes. The unknown solution within each element is approximated using a simple function, referred to as the basis function. These basis functions are usually polynomials of a certain degree, and their coefficients are determined by solving a set of linear equations.The mathematical theory of the FEM involves several key concepts and techniques. One of the fundamental principles is the variational formulation, which transforms the PDE into an equivalent variational problem. This variational problem is then discretized using the finite element approximation, resulting in a system of algebraic equations.Another important aspect is the assembly process, where the contributions from each element are combined to form the global stiffness matrix and right-hand side vector. This assembly is based on the integration of the basis functions and their derivatives over the element domains.Error estimation and convergence analysis are also essential components of the mathematical theory of the FEM. Various techniques, such as the energy method and the posteriori error estimators, are used to assess the accuracy of the finite element solution and to determine the appropriate mesh refinement for achieving convergence.Furthermore, the mathematical theory of the FEM includes the treatment ofboundary conditions, imposition of symmetries, and the development of efficient solvers for the resulting linear systems. It also addresses issues such as numerical stability,并行 computing, and adaptivity.In summary, the mathematical theory of the finite element method provides a comprehensive framework for numerically solving PDEs. It encompasses concepts such as element discretization, variational formulation, assembly, error estimation, and convergence analysis, which collectively enable the accurate and efficient simulation of a wide range of physical problems.。

Statistical mechanics of two-dimensional vortices and stellar systems

a r X iv:c o n d -ma t/212223v1[c ond-m at.stat-m ec h]1D ec22Statistical mechanics of two-dimensional vortices and stellar systems Pierre-Henri Chavanis Laboratoire de Physique Quantique,Universit´e Paul Sabatier,118,route de Narbonne 31062Toulouse,France Abstract.The formation of large-scale vortices is an intriguing phenomenon in two-dimensional turbulence.Such organization is observed in large-scale oceanic or atmo-spheric ﬂows,and can be reproduced in laboratory experiments and numerical simula-tions.A general explanation of this organization was ﬁrst proposed by Onsager (1949)by considering the statistical mechanics for a set of point vortices in two-dimensional hydrodynamics.Similarly,the structure and the organization of stellar systems (glob-ular clusters,elliptical galaxies,...)in astrophysics can be understood by developing a statistical mechanics for a system of particles in gravitational interaction as initiated by Chandrasekhar (1942).These statistical mechanics turn out to be relatively similar and present the same diﬃculties due to the unshielded long-range nature of the in-teraction.This analogy concerns not only the equilibrium states,i.e.the formation of large-scale structures,but also the relaxation towards equilibrium and the statistics of ﬂuctuations.We will discuss these analogies in detail and also point out the speciﬁcities of each system.1Introduction Two-dimensional ﬂows with high Reynolds numbers have the striking property of organizing spontaneously into coherent structures (the vortices)which dominate the dynamics [93](see Fig.1).The robustness of Jupiter’s Great Red Spot,a huge vortex persisting for more than three centuries in a turbulent shear between two zonal jets,is probably related to this general property.Some other coherent structures like dipoles (pairs of cyclone/anticyclone)and sometimes tripoles have been found in atmospheric or oceanic systems and can persist during several days or weeks responsible for atmospheric blocking.Some astrophysicists invoke the existence of organized vortices in the gaseous component of disk galaxies in relation with the emission of spiral density waves [99].It has also been proposed that planetary formation might have begun inside persistent gaseous vortices born out of the protoplanetary nebula [5,121,15,60,33](see Fig.2).As a result,hydrodynamical vortices occur in a wide of geophysical or astrophysical situations and their robustness demands a general understanding.Similarly,it is striking to observe that self-gravitating systems follow a kindof organization despite the diversity of their initial conditions and their en-vironement [9](see Fig.3).This organization is illustrated by morphological classiﬁcation schemes such as the Hubble sequence for galaxies and by simple2Pierre-Henri Chavanisrules which govern the structure of individual self-gravitating systems.For ex-ample,elliptical galaxies display a quasi-universal luminosity proﬁle described by de Vaucouleur’s R1/4law and most of globular clusters are well-ﬁtted by the Michie-King model.On the other hand,theﬂat rotation curves of spiral galaxies can be explained by the presence of a dark matter halo with a density proﬁle decreasing as r−2at large distances.The fractal nature of the interstellar medium and the large scale structures of the universe also display some form of organization.Fig.1.Self-organization of two-dimensional turbulentﬂows into large-scale vortices [93].These vortices are long-lived and dominate the dynamics.The question that naturally emerges is what determines the particular con-ﬁguration to which a self-gravitating system or a large-scale vortex settles.It is possible that their actual conﬁguration crucially depends on the conditions that prevail at their birth and on the details of their evolution.However,in view of their apparent regularity,it is tempting to investigate whether their organi-zation can be favoured by some fundamental physical principles like those of thermodynamics and statistical physics.We ask therefore if the actual states of self-gravitating systems in the universe and coherent vortices in two-dimensional2D vortices and stellar systems3Fig.2.A scenario of planet formation inside large-scale vortices presumably present in the Keplerian gaseous disk surrounding a star at its birth.Starting from a random vorticityﬁeld,a series of anticyclonic vortices appears spontaneously(upper panel). Due to the Coriolis force and to the friction with the gas,these vortices can eﬃciently trap dust particles passing nearby(lower pannel).The local increase of dust concen-tration inside the vortices can initiate the formation of planetesimals and planets by gravitational instability.This numerical simulation is taken from[15].turbulentﬂows are not simply more probable than any other possible conﬁgura-tion,i.e.if they cannot be considered as maximum entropy states.This statistical mechanics approach has been initiated by Onsager[101]for a system of point vortices and by Chandrasekhar[21]in the case of self-gravitating systems.It turns out that the statistical mechanics of two-dimensional vortices and self-gravitating systems present a deep analogy despite the very diﬀerent physical nature of these systems.This analogy was pointed out by Chavanis in[29,32,35] and further developed in[54,30,34,36,47,48].In the following,we will essentially discuss the statistical mechanics of2D vortices and refer to the review of Pad-manabhan[103](and his contribution in this book)for more details about the4Pierre-Henri Chavanisrge-scale structures in the universe as observed with the Hubble space tele-scope.The analogy with Fig.1is striking and will be discussed in detail in this paper. statistical mechanics of self-gravitating systems.We will see that the analogy be-tween two-dimensional vortices and(three-dimensional)self-gravitating systems concerns not only the prediction of the equilibrium state,i.e.the formation of large-scale structures,but also the statistics ofﬂuctuations and the relaxation towards equilibrium.This paper is organized as follows.In Sec.2,we discuss the statistical me-chanics of point vortices introduced by Onsager[101]and further developed by Joyce&Montgomery[70]and Pointin&Lundgren[107]among others(see a complete list of references in the book of Newton[98]).We discuss the existence of a thermodynamic limit in Sec.2.7and make the connexion withﬁeld theory. Statistical equilibrium states of axisymmetricﬂows are obtained analytically in Sec.2.8-2.9.The relation with equilibrium states of self-gravitating systems is shown in Sec.2.10.In Sec.3,we discuss the statistics of velocityﬂuctuations pro-duced by a random distribution of point vortices and use this stochastic approach to obtain an estimate of the diﬀusion coeﬃcient of point vortices.Application to2D decaying turbulence is considered in Sec.3.4.In Sec.4,we describe the relaxation of a point vortex in a thermal bath and analyze this relaxation in terms of a Fokker-Planck equation involving a diﬀusion and a drift.In Sec.5,we develop a more general kinetic theory of point vortices.A new kinetic equation is obtained which satisﬁes all conservation laws of the point vortex system and increases the Boltzmann entropy(H-theorem).We mention the connexion with2D vortices and stellar systems5 the kinetic theory of stars developed by Chandrasekhar[21].In Sec.6,we dis-cuss the violent relaxation of2D vortices and stellar systems.We mention the analogy between the Vlasov and the Euler equations and between the statistical approach developed by Lynden-Bell[90]for collisionless stellar systems and by Kuz’min[83],Miller[95]and Robert&Sommeria[111]for continuous vorticity ﬁelds.The concepts of“chaotic mixing”and“incomplete relaxation”are dis-cussed in the light of a relaxation theory in Sec.6.3.Application of statistical mechanics to geophysicalﬂows and Jupiter’s Great Red Spot are evocated in Sec.6.4.2Statistical mechanics of point vortices intwo-dimensional hydrodynamics2.1Two-dimensional perfectﬂowsThe equations governing the dynamics of an invisicidﬂow are the equation of continuity and the Euler equation:∂ρ∂t +(u∇)u=−16Pierre-Henri Chavaniswhere ∆=∂2xx +∂2yy is the Laplacian operator.In an unbounded domain,thisequation can be written in integral form asψ(r ,t )=−12πz × ω(r ′,t )r −r ′∂t +u ∇ω=0.(10)This corresponds to the transport of the vorticity ωby the velocity ﬁeld u .It is easy to show that the ﬂow conserves the kinetic energyE =u 22 (∇ψ)2d 2r =12 ωψd 2r ,(12)where the second equality is obtained by a part integration with the condition ψ=0on the boundary.Therefore,E can be interpreted either as the kinetic energy of the ﬂow (see Eq.(11))or as a potential energy of interaction between vortices (see Eq.(12)).2.2The point vortex gasWe shall consider the situation in which the velocity is created by a collection of N point vortices.In that case,the vorticity ﬁeld can be expressed as a sum of δ-functions in the formω(r ,t )=N i =1γi δ(r −r i (t )),(13)where r i (t )denotes the position of point vortex i at time t and γi is its circulation.According to Eqs.(9)(13),the velocity of a point vortex is equal to the sum of the velocities V (j →i )produced by the N −1other vortices,i.e.V i = j =i V (j →i )with V (j →i )=−γj |r j −r i |2.(14)2D vortices and stellar systems7 As emphasized by Kirchhoﬀ[79],the above dynamics can be cast in a Hamil-tonian formγi dx i∂y i,γidy i∂x i,(15) H=−12m).Thisis related to the particular circumstance that a point vortex is not a material particle.Indeed,an isolated vortex remains at rest contrary to a material particle which has a rectilinear motion due to its inertia.Point vortices form therefore a very peculiar Hamiltonian system.Note also that the Hamiltonian of point vortices can be either positive or negative(in the case of vortices of diﬀerent signs)whereas the kinetic energy of theﬂow is necessarily positive.This is clearlya drawback of the point vortex model.2.3The microcanonical approach of Onsager(1949)The statistical mechanics of point vortices wasﬁrst considered by Onsager[101] who showed the existence of negative temperature states at which point vortices of the same sign cluster into“supervortices”.He could therefore explain the formation of large,isolated vortices in nature.This was a remarkable anticipation since observations were very scarce at that time.Let us consider a liquid enclosed by a boundary,so that the vortices are conﬁned to an area A.Since the coordinates(x,y)of the point vortices are canonically conjugate,the phase space coincides with the conﬁguration space and isﬁnite:dx1dy1...dx N dy N= dxdy N=A N.(17)This striking property contrasts with most classical Hamiltonian systems con-sidered in statistical mechanics which have unbounded phase spaces due to the presence of a kinetic term in the Hamiltonian.As is usual in the microcanonical description of a system of N particles,we introduce the density of statesg(E)= dx1dy1...dx N dy Nδ E−H(x1,y1,...,x N,y N) ,(18)8Pierre-Henri Chavaniswhich gives the phase space volume per unit interaction energy E.The equilib-rium N-body distribution of the system,satisfying the normalization condition µ(r1,...,r N)d2r1...d2r N=1,is given byµ(r1,...,r N)=1T =dS2π i<jγiγj ln|r i−r j|.(22)2D vortices and stellar systems9 Making the change of variable x=r/R,weﬁnd that g(E,V)=V N g(E′,1) with E′=E+1βV 1+ββV 1+γ2(N−1)γ2.(25)We shall see in Sec.2.8that this negative critical inverse temperature is the minimum inverse temperature that the system can achieve.If,on the other hand,we consider a neutral system consisting of N/2vortices of circulationγand N/2vortices of circulation−γ,weﬁndP=N8π .(26)This result is well-known is plasma physics[113].The critical temperature at which the pressure vanishes is now positiveβc=8π∂t =−Ni=1γ∇δ(r−r i(t))V i.(28)Since V i=u(r i(t),t),we can rewrite the foregoing equation in the form∂ω10Pierre-Henri ChavanisSince the velocity is divergenceless,we obtain∂ω4π i=jγ2 ln|r i−r j|1=−2D vortices and stellar systems 11whereg (r 1,r 2,t )=µ(r 1,...,r N ,t )d 3r 2...d 2r N ,(36)is the two-body distribution function.In the mean-ﬁeld approximation,whichis exact in a properly deﬁned thermodynamic limit with N →+∞(see Sec.2.7),we haveg (r 1,r 2,t )=P (r 1,t )P (r 2,t ).(37)Accounting that N (N −1)≃N 2for large N ,the average energy takes the formE =−12ω ψd 2r =u 2n i !.(40)The logarithm of this number deﬁnes the Boltzmann ing Stirlingformula and considering the continuum limit in which ∆,ν→0,we get the classical formula S =−NP (r )ln P (r )d 2r ,(41)where P (r )is the density probability that a point vortex be found in the sur-face element centered on r .At equilibrium,the system is in the most probablemacroscopic state,i.e.the state that is the most represented at the microscopic level.This optimal state is obtained by maximizing the Boltzmann entropy (41)at ﬁxed energy (39)and vortex number N ,or total circulationΓ=Nγ=ω d 2r .(42)12Pierre-Henri ChavanisWriting the variational principle in the formδS−βδE−αδΓ=0,(43) whereβandαare Lagrange multipliers,it is readily found that the maximum entropy state corresponds to the Boltzmann distributionω =Ae−βγψ,(44) with inverse temperatureβ.We can account for the conservation of angular momentum L= ω r2d2r(in a circular domain)and impulse P= ω yd2r (in a channel)by introducing appropriate Lagrange multipliersΩand U foreach of these constraints.In that case,Eq.(44)remains valid provided that we replace the streamfunctionψby the relative streamfunctionψ′=ψ+Ω2π i<j ln|r i−r j|2D vortices and stellar systems13 where the potential of interaction has been normalized by R.For simplicity, we have ignored the contribution of the images but this shall not aﬀect the ﬁnal results.We now introduce the change of variables x=r/R and deﬁne the functionu(x1,...,x N)=1N2γ2.(49) In terms of these quantities,the density of states can be rewritteng(E)=2πV NNu(x1,...,x N) .(50)The proper thermodynamic limit for a system of point vortices with equal cir-culation in the microcanonical ensemble is such that N→+∞withﬁxedΛ.Wesee that the box size R does not enter in the normalized energyΛ.Therefore,the thermodynamic limit corresponds to N→+∞withγ∼N−1→0and E∼1. This is a very unusual thermodynamic limit due to the non-extensivity of thesystem.Note that the total circulationΓ=Nγremainsﬁxed in this process.For suﬃciently large N,the density of states can be writteng(E)≃ Dρe NS[ρ]δ Λ−E[ρ]) δ 1− ρ(r)d2r ,(51) withS[ρ]=− ρ(r)lnρ(r)d2r,(52) E[ρ]=−114Pierre-Henri Chavaniswhich is the normalization factor in the Gibbs measure1µ(r1,...,r N)=ln Z.(57)βUsing the notations introduced previously,we can rewrite the integral(55)in the formZ(β)=V N N0... N0N i=1d2x i eηu(x1,...,x N),(58)whereβNγ2η=2D vortices and stellar systems15 and phase transitions occur.This is the case,in particular,for the gravitational problem(see Sec.2.10).Finally,the grand canonical partition function is deﬁned byZ GC=+∞ N=0z N2(∇ξ)2−ρ(r)ξ(r)}d2r=e−12(∇ξ)2d2r+√−βγξ(r),we can easily carry out the summation on N to obtainZ GC= Dφe−12(∇φ)2−µ2eφ(r)},(67)T eff=−βγ2,µ2=−zβγ2.(68) Therefore,the grand partition function of the point vortex gas corresponds to a Liouvilleﬁeld theory with an actionS[φ]=12(∇φ)2−µ2eφ(r)}.(69)While the previous description is formally correct if we deﬁne Z and Z GC by Eqs.(55)and(64),it must be noted however that the canonical and grand canon-ical ensembles may not have a physical meaning for point vortices.In particular, it is not clear how one can impose a thermal bath at negative temperature. On the other hand,the usual procedure to derive the canonical ensemble from the microcanonical ensemble rests on a condition of additivity which is clearly lacking in the present case.2.8Axisymmetric equilibrium states in a diskLet us consider a collection of N point vortices with circulationγconﬁned within a disk of radius R.At statistical equilibrium,the streamfunctionψis solution of the Boltzmann-Poisson equation(45).If we work in a circular domain,we must in principle account for the conservation of angular momentum.This can lead to bifurcations between axisymmetric and oﬀ-axis solutions[117].We shall,16Pierre-Henri Chavanishowever,ignore this constraint for the moment in order to obtain analytical ex-pressions for the thermodynamical parameters.This is a suﬃcient approximation to illustrate the structure of the problem,which is our main concern here.If we conﬁne our attention to axisymmetric solutions,the Boltzmann-Poisson equation(45)can be written1dr r dψξddξ =λe−φ,(72)φ(0)=φ′(0)=0,(73)withλ=1ifβ<0andλ=−1ifβ>0.It turns out that this equation can be solved analytically as noticed by a number of authors.With the change of variables t=lnξandφ=2lnξ−z,Eq.(72)can be rewrittend2zdz(λe z).(74)This corresponds to the motion of aﬁcticious particle in a potential V(z)=λe z. This equation is readily integrated and,returning to original variables,weﬁnally obtaine−φ=18ξ2)2.(75)From the circulation theorem(6)applied to an axisymmetricﬂow,we have−dψ2πr.(76) whereΓ(r)= r0ω(r′)2πr′dr′is the circulation within r.Taking r=R and introducing the dimensionless variables previously deﬁned,we obtainη≡βγΓπR2(η+4)1η+4r22D vortices and stellar systems 17At positive temperatures (η>0),the vorticity is an increasing function of the distance and the vortices tend to accumulate on the boundary of the domain (Fig.4).On the contrary at negative temperatures (η<0),the vorticity is a decreasing function of the distance and the vortices tend to group themselves in the core of the domain to form a “supervortex”(Fig.5).These results are consistent with Onsager’s prediction [101].We also conﬁrm that statistical equilibrium states only exist for η>ηc =−4,as previously discussed.At this critical temperature,the central vorticity becomes inﬁnite and the solution tends to a Dirac peak:ω (r )→Γδ(r ),for η→ηc =−4.(79)00.20.40.60.81r/R1234567<ω>/ω∗η>0η=0η=20Fig.4.Statistical equilibrium states of point vortices at positive temperatures (ω∗=Γ/πR 2).The vortices are preferentially localized near the wall.The energy deﬁned by Eq.(39)can be written in the dimensionless formΛ≡2πE2η2 α0φ′(ξ)2ξdξ.(80)The integral can be carried out explicitly using Eq.(75).Eliminating αbetweenEqs.(80)and (77),we ﬁnd that the temperature is related to the energy by the equation of state Λ=1ηln418Pierre-Henri Chavanis00.20.40.60.81r/R246810<ω>/ω∗η<0η=−2η=−3η=−3.5Fig.5.Statistical equilibrium states of point vortices at negative temperatures showing a clustering.For η=−4,the vortices collapse at the center of the domain and the vorticity proﬁle is a Dirac peak.0.250.50.751Λ=2πE/Γ2−10−8−6−4−20246810η=βγΓ/2πΛ0=1/8ηc =−4Fig.6.Equilibrium phase diagram (caloric curve)for point vortices with equal circula-tion conﬁned within a disk.For simplicity,the angular momentum has not been taken into account (Ω=0).2D vortices and stellar systems1900.250.50.751Λ=2πE/Γ2−7−6−5−4−3S /NΛ0=1/8η=0η>0η<0ηc =−4Fig.7.Entropy vs energy plot for a system of point vortices with equal circulation conﬁned within a disk.The entropy (41)can also be calculated easily from the above results.Within an unimportant additive constant,it is given byS ηln 4−1+820Pierre-Henri Chavanisaccount,the density of states is given byg (E,L )=δ E −H (r 1,...,r N ) δ L −N i =1γr 2i N i =1d 2r i ,(83)and the angular velocity of the ﬂow byΩ=2T∂SL ,oneﬁnds the exact resultΩ=2N8π.(85)Therefore,the vorticity ﬁeld is determined by the Boltzmann distributionω =Ae −βγψ′,(86)where ψ′is the relative streamfunctionψ′≡ψ+Ω4βL(4+η)r 2.(87)For η=0,one hasω =γNΓLr 2.(88)For large r ,the asymptotic behavior of Eq.(86)isω ∼14L(4+η)r 2(r →+∞),(89)where we have used ψ∼−(Γ/2π)ln r at large distances.From Eq.(89),one sees that η≥−4is required for the existence of an integrable solution.Inserting the relation (86)in the Poisson equation (7),we get−∆ψ′=Ae −2πη2πLη(4+η).(90)With the change of variablesξ=γNN 2γ2−2πηdξ2+1dξ=2πηe φ−(4+η),(92)2D vortices and stellar systems21 and the vorticity(86)becomesω =N2γ22πL,if r≤(2L/Γ)1/2,(94)and ω =0otherwise.This vortex patch is the state of minimum energy at ﬁxed circulation and angular momentum.Forη→−4,one has approximatelyω =N2γ2(1−AπηγN4L(4+η)r2.(95)Theﬁrst factor is an exact solution of Eq.(92)with the second term on the right hand side neglected(see Sec.2.8).The second factor is a correction for large r, in agreement with the asymptotic result expressed by Eq.(89).The parameter A tends to inﬁnity asη→−4and is determined from the condition ω d2r=Γby the formulaπA+ln(πA)=−C−ln 1+η|r j−r i|3,(97)where F(j→i)is the force created by star j on star i.The force can be written as the gradient F=−∇Φof a gravitational potentialΦwhich is related to the stellar densityρ(r,t)=Ni=1mδ(r−r i),(98)by the Poisson equation∆Φ=4πGρ.(99)22Pierre-Henri ChavanisFurthermore,the equations of motion(Newton’s equations)can be put in the Hamiltonian formm d r i∂v i,md v i∂r i,H=1|r i−r j|.(100)In the analogy between stellar systems and two-dimensional vortices,the star densityρplays the role of the vorticityω,the force F the role of the velocity V and the gravitational potentialΦthe role of the streamfunctionψ.The crucial point to realize is that,for the two systems,the interaction is a long-range un-shielded Coulombian interaction(in D=3or D=2dimensions).This makes the connexion between point vortices and stellar systems deeper than between point vortices and electric charges for example.In particular,point vortices can organize into large scale clusters,like stars in galaxies,while the distribution of electric charges in a neutral plasma is uniform.There are,on the other hand, fundamental diﬀerences between stars and vortices.In particular,a star creates an acceleration while a vortex creates a velocity.On the other hand,the gravi-tational interaction is attractive and directed along the line joining the particles while the interaction between vortices is rotational and perpendicular to the line joining the vortices.Despite these important diﬀerences,the statistical mechanics of2D vortices and stellar systems are relatively similar.Like the point vortex gas,the self-gravitating gas is described at statistical equilibrium by the Boltzmann distri-butionρ =Ae−βΦ,(101) obtained by maximizing the Boltzmann entropy atﬁxed mass M and energy E. Its structure is therefore determined by solving the Boltzmann-Poisson equation∆Φ=4πGAe−βΦ,(102) where A andβ>0have to be related to M and E.This statistical mechanics approach has been developed principally for globular clusters relaxing towards equilibrium via two-body encounters[9].It is clear that the Boltzmann-Poisson equation(102)is similar to the Boltzmann-Poisson equation(45)for point vor-tices at negative temperatures.The density proﬁle determined by these equations is a decreasing function of the distance,which corresponds to a situation of clus-tering(see Figs.5and8).The similarity of the maximum entropy problem for stars and vortices,and the Boltzmann-Poisson equations(102)(45),is aﬁrst manifestation of the formal analogy existing between these two systems.However,due to the diﬀerent dimension of space(D=3for stars instead of D=2for vortices),the mathematical problems diﬀer in the details.First of all,the density proﬁle determined by the Boltzmann-Poisson equation(102)in D=3decreases like r−2at large distances leading to the so-called inﬁnite mass2D vortices and stellar systems 23−50510ln(ξ)−20−15−10−50ln (e −ψ)singularsphereξ−2Fig.8.Density proﬁle of the self-gravitating gas at statistical equilibrium.The dashed line corresponds to the singular solution ρ=1/2πGβr 2.problem since M = +∞0ρ4πr 2dr →+∞[19].There is no such problem forpoint vortices in two dimensions:the vorticity decreases like r −4,or even more rapidly if the conservation of angular momentum is accounted for,and the totalcirculation Γ= +∞0ω2πrdr is ﬁnite.The inﬁnite mass problem implies thatno statistical equilibrium state exists for open star clusters,even in theory.A system of particles in gravitational interaction tends to evaporate so that the ﬁnal state is just two stars in Keplerian orbit.This evaporation process has been clearly identiﬁed in the case of globular clusters which gradually lose stars to the beneﬁt of a neighboring galaxy.In fact,the evaporation is so slow that we can consider in a ﬁrst approximation that the system passes by a succession of quasiequilibrium states described by a truncated isothermal distribution function (Michie-King model)[9].This justiﬁes the statistical mechanics approach in that sense.Another way of solving the inﬁnite mass problem is to conﬁne the system within a box of radius R .However,even in that case,the notion of equilibrium poses problem regarding what now happens at the center of the conﬁguration.The equilibrium phase diagram (E,T )for bounded self-gravitating systems is represented in Fig.9.The caloric curve has a striking spiral behavior parametrized by the density contrast R =ρ(0)/ρ(R )going from 1(homogeneous system)to +∞(singular sphere)as we proceed along the spiral.There is no equilibrium state below E c =−0.335GM 2/R or T c =GMm24Pierre-Henri ChavanisΛ=−ER/GM 20.51.52.5η=βG M /R Fig.9.Equilibrium phase diagram for self-gravitating systems conﬁned within a box.For suﬃciently low energy or temperature,there is no equilibrium state and the system undergoes gravitational collapse.−1−0.500.51 1.5Λ=−ER/GM 200.511.522.533.5η=βG M /R µ=105µ=104µ=103µ=102µ=10Fig.10.Equilibrium phase diagram for self-gravitating fermions [41].The degeneracy parameter µplays the role of a small-scale cut-oﬀǫ∼1/µ.For ǫ→0,the classical spiral of Fig.9is recovered.2D vortices and stellar systems25 mass goes to zero.Therefore,the singularity contains no mass and this process alone cannot lead to a black hole.Since the T(E)curve has turning points,this implies that the microcanon-ical and canonical ensembles are not equivalent and that phase transitions will occur[103].In the microcanonical ensemble,the series of equilibria becomes un-stable after theﬁrst turning point of energy(MCE)corresponding to a density contrast of709.At that point,the solutions pass from local entropy maxima to saddle points.In the canonical ensemble,the series of equilibria becomes unsta-ble after theﬁrst turning point of temperature(CE)corresponding to a density contrast of32.1.At that point,the solutions pass from minima of free energy (F=E−T S)to saddle points.It can be noted that the region of negative speciﬁc heats between(CE)and(MCE)is stable in the microcanonical en-semble but unstable in the canonical ensemble,as expected on general physical grounds.The thermodynamical stability of isothermal spheres can be deduced from the topology of theβ−E curve by using the turning point criterion of Katz[75]who has extended Poincar´e’s theory of linear series of equilibria.The stability problem can also be reduced to the study of an eigenvalue equation associated with the second order variations of entropy or free energy as studied by Padmanabhan[102]in the microcanonical ensemble and by Chavanis[37]in the canonical ensemble.This study has been recently extended to other statisti-cal ensembles[44]:grand canonical,grand microcanonical,isobaric....The same stability limits as Katz are obtained but this method provides in addition the form of the density perturbation proﬁles that trigger the instability at the critical points.It also enables one to show a clear equivalence between thermodynamical stability in the canonical ensemble and dynamical stability with respect to the Navier-Stokes equations(Jeans problem)[37,44].These analytical methods can be extended to general relativity[38].It must be stressed,however,that the statistical equilibrium states of self-gravitating systems are at most metastable: there is no global maximum of entropy or free energy for a classical system of point masses in gravitational interaction[2].Phase transitions in self-gravitating systems can be studied in detail by in-troducing a small-scale cut-oﬀǫin order to regularize the potential.This can be achieved for example by considering a system of self-gravitating fermions (in which case an eﬀective repulsion is played by the Pauli exclusion principle) [62,8,52,41,43]or a hard spheres gas[3,103,120].Other forms of regularization are possible[59,128,45].For these systems,there can still be gravitational collapse but the core will cease to shrink when it feels the inﬂuence of the cut-oﬀ.The result is the formation of a compact object with a large mass:a“fermion ball”or a hard spheres“condensate”.The equilibrium phase diagram of self-gravitating fermions is represented in Fig.10and has been discussed at length by Chavanis [41]in the light of an analytical model.The introduction of a small-scale cut-oﬀhas the eﬀect of unwinding the classical spiral of Fig.9.For a small cut-oﬀǫ≪1,the trace of the spiral is still visible and the T(E)curve is multivalued (Fig.11).This can lead to a gravitationalﬁrst order phase transition between a gaseous phase with an almost homogeneous density proﬁle(upper branch)and a。

UDEC中文指导说明

通用离散元用户指导（U D E C 3.1）2004.9目录1 引言 (1)1.1 总论 (1)1.2 与其他方法的比较 (2)1.3 一般特性 (2)1.4 应用领域 (3)2 开始启动 (4)2.1 安装和启动程序 (4)2.1.7 内存赋值 (4)2.1.9 运行UDEC (5)2.1.10 安装测试程序 (5)2.2 简单演示－通用命令的应用 (5)2.3 概念与术语 (6)2.4 UDEC模型：初始块体的划分 (8)2.5 命令语法 (9)2.6 UDEC应用基础 (10)2.6.1 块体划分 (10)2.6.2 指定材料模型 (16)2.6.2.1 块体模型 (16)2.6.2.2 节理模型 (17)2.6.3 施加边界条件和初始条件 (19)2.6.4 迭代为初始平衡 (21)2.6.5 进行改变和分析 (24)2.6.6 保存或恢复计算状态 (25)2.6.7 简单分析的总结 (25)2.8 系统单位 (26)3 用UDEC求解问题 (27)3.1 一般性研究 (27)3.1.1 第1步：定义分析模型的对象 (28)3.1.2 第2步：产生物理系统的概念图形 (28)3.1.3 第3步：建造和运行简单的理想模型 (28)3.1.4 第4步：综合特定问题的数据 (29)3.1.5 第5步：准备一系列详细的运行模型 (29)3.1.6 第6步：进行模型计算 (29)3.1.7 第7步：提供结果和解释 (30)3.2 产生模型 (30)3.2.1 确定UDEC模型合适的计算范围 (30)3.2.2 产生节理 (32)3.2.2.1 统计节理组生成器 (32)3.2.2.2 VORONOI多边形生成器 (34)3.2.2.3 例子 (34)3.2.3 产生内部边界形状 (35)3.3 变形块体和刚体的选择 (38)3.4 边界条件 (42)3.4.1 应力边界 (42)3.4.1.1 施加应力梯度 (43)3.4.1.2 改变边界应力 (44)3.4.1.3 打印和绘图 (44)3.4.1.4 提示和建议 (45)3.4.2 位移边界 (46)3.4.3 真实边界－选择合理类型 (46)3.4.4 人工边界 (46)3.4.4.1 对称轴 (46)3.4.4.2 截取边界 (46)3.4.4.3 边界元边界 (49)3.5 初始条件 (50)3.5.1 在均匀介质中的均匀应力：无重力 (50)3.5.2 无节理介质中具有梯度变化的应力：均匀材料 (51)3.5.3 无节理介质中具有梯度变化的应力：非均匀材料 (51)3.5.4 具有非均匀单元的密实模型 (52)3.5.5 随模型变化的初始应力 (53)3.5.6 节理化介质的应力 (54)3.5.7 绘制应力等值线图 (55)3.6 加载与施工模拟 (57)3.7 选择本构模型 (62)3.7.1 变形块体材料模型 (63)3.7.2 节理材料模型 (64)3.7.3 合理模型的选择 (65)3.8 材料性质 (71)3.8.1 岩块性质 (71)3.8.1.1 质量密度 (71)3.8.1.2 基本变形性质 (71)3.8.1.3 基本强度性质 (72)3.8.1.4 峰后效应 (73)3.8.1.5 现场性质参数的外延 (77)3.8.2 节理性质 (80)3.9 提示和建议 (81)3.9.1 节理几何形状的选择 (81)3.9.2 设计模型 (81)3.9.3 检查模型运行时间 (82)3.9.4 对允许时间的影响 (82)3.9.5 单元密度的考虑 (83)3.9.6 检查模型响应 (83)3.9.7 检查块体接触 (83)3.9.8 应用体积模量和剪切模量 (83)3.9.9 选择阻尼 (84)3.9.10 给块体和节理模型指定模型和赋值 (84)3.9.11 避免圆角误差 (85)3.9.12 接触嵌入 (85)3.9.13 非联结块体 (86)3.9.14 初始化变量 (86)3.9.15 确定坍塌荷载 (86)3.9.16 确定安全系数 (86)3.10 解释 (88)3.10.1 不平衡力 (88)3.10.2 块体/网格结点的速度 (88)3.10.3 块体破坏的塑性指标 (89)3.11 模拟方法 (89)3.11.1 有限数据系统模拟 (89)3.11.2 混沌系统的模拟 (90)3.11.3 局部化、物理的不稳定性和应力路径 (91)1 引言1.1 总论通用离散元程序（UDEC，Universal Distinct Element Code）是一个处理不连续介质的二维离散元程序。

Advances in discrete element modelling of underground

REVIEW ARTICLEAdvances in discrete element modelling of underground excavationsCarlos Labra ÆJerzy Rojek ÆEugenio On˜ate ÆFrancisco ZarateReceived:5November 2007/Accepted:6May 2008/Published online:17July 2008ÓSpringer-Verlag 2008Abstract The paper presents advances in the discrete element modelling of underground excavation processes extending modelling possibilities as well as increasing computational efﬁciency.Efﬁcient numerical models have been obtained using techniques of parallel computing and coupling the discrete element method with ﬁnite element method.The discrete element algorithm has been applied to simulation of different excavation processes,using dif-ferent tools,TBMs and roadheaders.Numerical examples of tunnelling process are included in the paper,showing results in the form of rock failure,damage in the material,cutting forces and tool wear.Efﬁciency of the code for solving large scale geomechanical problems is also shown.Keywords Coupling ÁDiscrete element method ÁFinite element method ÁParallel computation ÁTunnelling1IntroductionA discrete element algorithm is a numerical technique which solves engineering problems that are modelled as a large system of distinct interacting bodies or particles that are subject to gross motion.The discrete element method (DEM)is widely recognized as a suitable tool to model geomaterials [1,2,4,8].The method presents important advantages in simulation of strong discontinuities such as rock fracturing during an underground excavation or rock failure induced by a tunnel excavation.It is difﬁcult to solve such problems using conventional continuum-based procedures such as the ﬁnite element method (FEM).The DEM makes possible the simulation of different excavation processes [5,7]allowing the determination of the damage of the rock or soil,or evaluation of cutting forces in rock excavation with roadheaders or TBMs.Different possibil-ities of DEM applications in simulation of tunnelling process are shown in the paper.Examples include new developments like evaluation of tool wear in rock cutting processes.The main problem in a wider use of this method is the high computational cost required by the simulations ﬁrst of all due to large number of discrete elements usually required.Different strategies are possible in addressing this problem.This paper will present two approaches:parall-elization and coupling the DEM and FEM.Parallelization techniques are useful for the simulation of large-scale problems,where the number of particles involved does not allow the use of a single processor,or where the single processor calculation would require an extremely long time.A shared memory parallelization of the DEM algorithm is presented in the paper.A high per-formance code for the simulation of tunnel construction problems is described and examples of the efﬁciency of thebra ÁE.On˜ate ÁF.Zarate International Center for Numerical Methods in Engineering,Technical University of Catalonia,Gran Capitan s/n,08034Barcelona,Spaine-mail:clabra@ E.On˜ate e-mail:onate@ F.Zaratee-mail:zarate@J.Rojek (&)Institute of Fundamental Technological Research,PolishAcademy of Sciences,Swietokrzyska 21,00049Warsaw,Poland e-mail:jrojek@.plActa Geotechnica (2008)3:317–322DOI 10.1007/s11440-008-0071-2code for solving large-scale geomechanical problems are shown in the paper.In many cases discontinuous material failure is localized in a portion of the domain,the rest of it can be treated as continuum.Continuous material is usually modelled more efﬁciently using the FEM.In such problems coupling of the discrete element method with the FEM can provide an optimum solution.Discrete elements are used only in a portion of the analysed domain where material fracture occurs,while outside the DEM subdomainﬁnite elements can be bining these two methods in one model of rock cutting allows us to take advantages of each method. The paper presents a coupled discrete/ﬁnite element tech-nique to model underground excavation employing the theoretical formulation initiated in[5]and further devel-oped in[6].2Discrete element method formulationThe discrete element model assumes that material can be represented by an assembly of distinct particles or bodies interacting among themselves.Generally,discrete elements can have arbitrary shape.In this work the formulation employing cylindrical(in2D)or spherical(in3D)rigid particles is used.Basic formulation of the discrete element formulation using spherical or cylindrical particles wasﬁrst proposed by Cundall and Strack[1].Similar formulation has been developed by the authors[5,7]and implemented in the explicit dynamic code Simpact.The code has a lot of original features like modelling of tool wear in rock cut-ting,thermomechanical coupling and other capabilities not present in commercial discrete element codes.Translational and rotational motion of rigid spherical or cylindrical elements is described by means of the Newton–Euler equations of rigid body dynamics:M D€r D¼F D;J D_X D¼T Dð1Þwhere r D is the position vector of the element centroid in a ﬁxed(inertial)coordinate frame,X D is the angular veloc-ity,M D is the diagonal matrix with the element mass on the diagonal,J D is the diagonal matrix with the element moment of inertia on the diagonal,F D is the vector of resultant forces,and T D is the vector of resultant moments about the element central axes.Vectors F D and T D are sums of all forces and moments applied to the element due to external load,contact interactions with neighbouring spheres and other obstacles,as well as forces resulting from damping in the system.Equations of motion(1)are inte-grated in time using the central difference scheme.The overall behaviour of the system is determined by the cohesive/frictional contact laws assumed for the inter-action between contacting rigid spheres(or discs in2D).The contact law can be seen as the formulation of the material model on the microscopic level.Modelling of rock or cohesive zones requires contact models with cohesion allowing tensile interaction force between particle.In the present work the simplest of the cohesive models,the elastic perfectly brittle model is used.This model is char-acterized by linear elastic behaviour when cohesive bonds are active:r¼k n u n;s¼k t u tð2Þwhere r and s are the normal and tangential contact force, respectively,k n and k t are the interface stiffness in the normal and tangential directions and u n and u t the normal and tangential relative displacements,respectively. Cohesive bonds are broken instantaneously when the interface strength is exceeded in the tangential direction by the tangential contact force or in the normal direction by the tensile contact force.The failure(decohesion)criterion is written as:r R n;j s j R t;ð3Þwhere R n and R t are the interface strengths in the normal and tangential directions,respectively.Breakage of cohe-sive bonds allows us to simulate fracture of material and its propagation.In the absence of cohesion the frictional contact is assumed with the Coulomb friction model.3Coupling the DEM and FEMIn the present work the so-called explicit dynamic formu-lation of the FEM is used.The explicit FEM is based on the solution of discretized equations of motion written in the current conﬁguration in the following form:M F€r F¼F ext FÀF int Fð4Þwhere M F is the mass matrix,r F is the vector of nodal displacements,F F ext and F F int are the vectors of external loads and internal forces,respectively.Similarly to the DEM algorithm,the central difference scheme is used for time integration of(4).It is assumed that the DEM and FEM can be applied in different subdomains of the same body.The DEM and FEM subdomains,however,do not need to be disjoint—they can overlap each other.The common part of the subdomains is the part where both discretization types are used with gradually varying contribution of each modelling method.This idea follows that used for molecular dynamics coupling with a continuous model in[9].The coupling of DEM and FEM subdomains is provided by additional kinematical constraints.Interface discrete elements are constrained by the displacementﬁeld of overlapping interfaceﬁnite elements.Making use of thesplit of the global vector of displacements of discrete ele-ments,r D ,into the unconstrained part,r DU ,and the constrained one,r DC ,r D ={r DU ,r DC }T ,additional kine-matic relationships can be written jointly in the matrix notation as follows:v ¼r DC ÀNr F ¼0;ð5Þwhere N is the matrix containing adequate shape functions.Additional kinematic constraints (5)can be imposed by the Lagrange multiplier or penalty method.The set of equations of motion for the coupled DEM/FEM system with the penalty coupling is as follows"M F 0000"M DU 0000"M DC 0000"J D26643775€r F €r DU €r DC _X D 8>><>>:9>>=>>;¼"F ext F À"F int F þN T k DF v "F DU"F DC Àk DF v "T D 8>><>>:9>>=>>;ð6Þwhere k DF is the diagonal matrix containing on its diagonal the values of the discrete penalty function,and globalmatrices "M F ;"M DU ;"M DC and "J D ;and global vectors"F int F ;"F ext F ;"F DU ;"F DC and "T D are obtained by aggregation of adequate elemental matrices and vectors taking into account appropriate contributions from the discrete and ﬁnite element parts.Equation (6)can be integrated in time using the standard central difference scheme.4Application of DEM to simulation of tunnelling process Fracture of rock or soil as well as interaction between a tunnelling machine and rock during an excavation process can be simulated by means of the DEM.This kind of analysis enables the comparison of the excavation process under different conditions.4.1Simulation of tunnelling with a TBMSimpliﬁed models of a tunnelling process must be used due to a high computational cost of a full-scale simulation in this case.We assume that the TBM is modelled as a cylinder with a special contact model for the tunnel face is adopted.Figure 1presents a simpliﬁed tunnelling process.The rock sample,with a diameter of 10m and a length of 7m,is discretized with randomly generated and densely com-pacted 40,988spheres.Discretization of the TBM geometry employs 1,193rigid triangular elements.Tunnelling pro-cess has been carried out with prescribed horizontal velocity 5m/h and rotational velocity of 10rev/min.Rock properties of granite are used,and the microscopic DEM parameters corresponding to the macroscopic granite properties are obtained using the methodology described in [10].A special condition is adopted to eliminate the spherical particles in the face of the tunnel.Each particle,which is in contact with the TBM and lacks cohesive contacts with other particles,is removed from the model.Thus,the advance of the TBM and the absorption of the material in the shield of the TBM is modelled.Figure 1a,c presents the displacement of the TBM and the elimination of the rock material.The area affected by the loss of cohesive contacts,resulting in material failure is shown in Fig.2.This loss of cohesion can be considered as damage ,because it produces the change of the equivalent Young modulus.4.2Simulation of linear cutting test of single disccutter Simulation of the linear cutting test was performed.A rock sample with dimensions of 13591095cm is repre-sented by an assembly of randomly generated and densely compacted 40,449spherical elements of radii ranging from 0.08to 0.60cm.The granite properties are assumed in the simulation and appropriate DEM parameters are evaluated.The disc cutter is treated as a rigid body and the parameters describing its interaction with the rock are as follows:contact stiffness modulus k n =10GPa,Coulomb friction coefﬁcient l =0.8.The velocity of the disc cutter is assumed to be 10m/s.Fig.1Simulation of TBM excavation:Evolution and elimination ofmaterialFig.2Simulation of TBM excavation:Damage over tunnel surfaceFigure 3a shows the discretization of the disc cutter.Only the area of the cutter ring in direct interaction with the rock is discretized with discrete elements due to the com-putational cost reasons.The whole model is presented in Fig.3b.The evolution of the normal cutting force during the process is depicted in Fig.4a.The values of the forces should be validated,because the boundary condition can affect the results.The evolution of the wear,using the for-mulation presented in [5],can be seen in Fig.4b.The elimination of the discrete elements,where the wear exceed the prescribed limit,permit the modiﬁcation of the disc cutter shape,which leads to a change of the interaction forces.In the present case,a low value of the wear constant is considered,in order to maintain the initial tool shape.Accumulated wear indicates the areas where the removal of the tool material is most intensive.An acceleration of the wear process using higher values of the wear constant is required in order to obtain in a short time considered in the analysis the amount of wear equivalent to real working time.5High performance simulationsOne of the main problems with the DEM simulation is the computational cost.The contact search,the force calcula-tion for each contact,and the large number of elements necessary to resolve a real life problem requires a high computational effort.High performance computation,and parallel implementation could be necessary to run simu-lations with large number of time steps.The advances of the computer capabilities during last years and the use of multiprocessors techniques enable the use of parallel computing methods for the discrete element analysis of large scale real problems.A shared memory parallel version of the code is tested.The main idea is to make a partition of the mesh of particles and use each processor for the contact calculation at different parts of the mesh.The partition process is performed using a special-ized library [3].The calculation of the cohesive contacts requires most of the computational cost.A special structure for the database,and the dynamic load balance is used in order to obtain a good performance for the simulations.Two different structures for the contact data are used in order to have a good management of the information.The ﬁrst data structure is created for the initial cohesive contacts,where a static array can be used.The other data structureisFig.3Linear cutting test simulation:a cutter ring with partial discretization;b full discretized model0.014Table 1Times for different number of processors Time (s)versus processors 124Total404.31272.93156.85Static contacts (per step)0.12790.06920.0351Dynamic contacts (per step)0.00590.00570.0055Time integration (per step)0.04260.03570.0344Speed up1.001.842.58designed for the dynamic contacts,occurring in the process of rock fragmentation,and the interaction between differ-ent bodies.The management of this kind of contact is completely dynamic,and it is not necessary to store vari-ables with the history information.Table 1presents the times of parallel simulations of a tunnelling process,which was described earlier.The main computational cost is due to the cohesive contacts evalu-ation.The results shown in the table conﬁrm that a good speed-up has been achieved.6DEM and DEM/FEM simulation of rock cutting A process of rock cutting with a single pick of a roadheader cutter-head has been simulated using discrete and hybrid discrete/ﬁnite element models.In the hybrid DEM/FEM model discrete elements have been used in the part of rock mass subjected to fracture,while the other part have been discretized with ﬁnite elements.In both models the tool is considered rigid,assuming the elasticity of the tool is irrelevant for the purpose of modelling of rock fracture.Figure 5presents results of DEM and DEM/FEM sim-ulation.Both models produce similar failures of rock during cutting.Cutting forces obtained using these two models are compared in Fig.6.Both curves show oscilla-tions typical for cutting of brittle rock.In both cases similar values of amplitudes are observed.Mean values of cutting forces agree very well.This shows that combined DEM/FEM simulation gives similar results to a DEM analysis,while being more efﬁcient numerically—computation time has been reduced by half.7Conclusions •Discrete element method using spherical or cylindrical rigid particles is a suitable tool in modelling of underground excavation processes.•Use of the model in a particular case requires calibra-tion of the discrete element model using available experimental results.•Discrete element simulations of real engineering prob-lems require large computation time and memory resources.•Efﬁciency of discrete element computation can be improved using technique of parallel computations.Parallelization makes possible the simulation of large problems.•The combination of discrete and ﬁnite elements is an effective approach for simulation of underground rock excavation.Acknowledgments The work has been sponsored by the EU project TUNCONSTRUCT (contract no.IP 011817-2)coordinated by Prof.G.Beer (TU Graz,Austria).References1.Cundall PA,Strack ODL (1979)A discrete numerical method for granular assemblies.Geotechnique29:47–65Fig.5Simulation of rock cutting:a DEM model,b DEM/FEM model2.Campbell CS(1990)Rapid granularﬂows.Annu Rev Fluid Mech2:57–923.Karypis G,Kumar V(1998)A fast and high quality multilevelscheme for partitioning irregular graphs.SIAM J Sci Comput 20:359–3924.Mustoe G(ed)(1992)Eng Comput9(2).Special issue5.On˜ate E,Rojek J(2004)Combination of discrete element andﬁnite element methods for dynamic analysis of geomechanics put Methods Appl Mech Eng193:3087–3128 6.Rojek J(2007)Modelling and simulation of complex problems ofnonlinear mechanics using theﬁnite and discrete element meth-ods(in Polish).Habilitiation Thesis,Institute of Fundamental Technological Research Polish Academy of Sciences,Warsaw7.Rojek J,On˜ate E,Zarate F,Miquel J(2001)Modelling of rock,soil and granular materials using spherical elements.In:2nd European conference on computational mechanics ECCM-2001, Cracow,26–29June8.Williams JR,O’Connor R(1999)Discrete element simulationand the contact problem.Arch Comp Meth Eng6(4):279–304 9.Xiao SP,Belytschko T(2004)A bridging domain method forcoupling continua with molecular put Methods Appl Mech Eng193:1645–166910.Zarate F,Rojek J,On˜ate E,Labra C(2007)A methodology todetermine the particle properties in2d and3d dem simulations.In:ECCOMAS thematic conference on computational methods in tunnelling EURO:TUN-2007,Vienna,Austria,27–29August。

弹性力学第4章

5Formulation and Solution StrategiesThe previous chapters have now developed the basicﬁeld equations of elasticity theory.Theseresults comprise a system of differential and algebraic relations among the stresses,strains,anddisplacements that express particular physics at all points within the body under investigation.In this chapter we now wish to complete the general formulation byﬁrst developing boundaryconditions appropriate for use with theﬁeld equations.These conditions specify the physicsthat occur on the boundary of body,and generally provide the loading inputs that physicallycreate the interior stress,strain,and displacementﬁelds.Although theﬁeld equations are thesame for all problems,boundary conditions are different for each problem.Therefore,properdevelopment of boundary conditions is essential for problem solution,and thus it is importantto acquire a good understanding of such development bining theﬁeldequations with boundary conditions then establishes the fundamental boundary value problemsof the theory.This eventually leads us into two different formulations,one in terms ofdisplacements and the other in terms of stresses.Because these boundary value problems aredifﬁcult to solve,many different strategies have been developed to aid in problem solution.Wereview in a general way several of these strategies,and later chapters incorporate many ofthese into the solution of speciﬁc problems.5.1Review of Field EquationsBefore beginning our discussion on boundary conditions we list here the basicﬁeld equationsfor linear isotropic elasticity.Appendix A includes a more comprehensive listing of allﬁeldequations in Cartesian,cylindrical,and spherical coordinate systems.Because of its ease of useand compact properties,our formulation uses index notation.Strain-displacement relations:e ij¼12(u i,jþu j,i)(5:1:1)Compatibility relations:e ij,klþe kl,ijÀe ik,jlÀe jl,ik¼0(5:1:2)83Equilibrium equations:s ij,jþF i¼0(5:1:3) Elastic constitutive law(Hooke’s law):s ij¼l e kk d ijþ2m e ije ij¼1þnEs ijÀnEs kk d ij(5:1:4)As mentioned in Section2.6,the compatibility relations ensure that the displacements arecontinuous and single-valued and are necessary only when the strains are arbitrarily speciﬁed.If,however,the displacements are included in the problem formulation,the solution normallygenerates single-valued displacements and strain compatibility is automatically satisﬁed.Thus,in discussing the general system of equations of elasticity,the compatibility relations(5.1.2)are normally set aside,to be used only with the stress formulation that we discuss shortly.Therefore,the general system of elasticityﬁeld equations refers to the15relations(5.1.1),(5.1.3),and(5.1.4).It is convenient to deﬁne this entire system using a generalized operatornotation asJ{u i,e ij,s ij;l,m,F i}¼0(5:1:5) This system involves15unknowns including3displacements u i,6strains e ij,and6stresses s ij.The terms after the semicolon indicate that the system is also dependent on two elastic materialconstants(for isotropic materials)and on the body force density,and these are to be given apriori with the problem formulation.It is reassuring that the number of equations matches thenumber of unknowns to be determined.However,this general system of equations is of suchcomplexity that solutions by using analytical methods are essentially impossible and furthersimpliﬁcation is required to solve problems of interest.Before proceeding with development ofsuch simpliﬁcations,it is usefulﬁrst to discuss typical boundary conditions connected with theelasticity model,and this leads us to the classiﬁcation of the fundamental problems.5.2Boundary Conditions and Fundamental ProblemClassiﬁcationsSimilar to otherﬁeld problems in engineering science(e.g.,ﬂuid mechanics,heat conduction,diffusion,electromagnetics),the solution of system(5.1.5)requires appropriate boundaryconditions on the body under study.The common types of boundary conditions for elasticityapplications normally include speciﬁcation of how the body is being supported or loaded.Thisconcept is mathematically formulated by specifying either the displacements or tractions atboundary points.Figure5-1illustrates this general idea for three typical cases includingtractions,displacements,and a mixed case for which tractions are speciﬁed on boundary S tand displacements are given on the remaining portion S u such that the total boundary is givenby S¼S tþS u.Another type of mixed boundary condition can also occur.Although it is generally not possible to specify completely both the displacements and tractions at the same boundarypoint,it is possible to prescribe part of the displacement and part of the traction.Typically,this 84FOUNDATIONS AND ELEMENTARY APPLICATIONStype of mixed condition involves the speciﬁcation of a traction and displacement in two different orthogonal directions.A common example of this situation is shown in Figure 5-2for a case involving a surface of problem symmetry where the condition is one of a rigid-smooth boundary with zero normal displacement and zero tangential traction.Notice that in this example the body under study was subdivided along the symmetry line,thus creating a new boundary surface and resulting in a smaller region to analyze.Because boundary conditions play a very essential role in properly formulating and solving elasticity problems,it is important to acquire a clear understanding of their speciﬁcation and use.Improper speciﬁcation results in either no solution or a solution to a different problem than what was originally sought.Boundary conditions are normally speciﬁed using the coordinate system describing the problem,and thus particular components of the displacements and tractions are set equal to prescribed values.For displacement-type conditions,such a speciﬁ-cation is straightforward,and a common example includes ﬁxed boundaries where the dis-placements are to be zero.For traction boundary conditions,the speciﬁcation can be a bit more complex.Figure 5-3illustrates particular cases in which the boundaries coincide with Cartesian or polar coordinate surfaces.By using results from Section 3.2,the traction speciﬁcation can be reduced to a stress speciﬁcation.For the Cartesian example in which y ¼constant ,Displacement Conditions Mixed ConditionsTraction ConditionsFIGURE 5-1Typical boundary conditions.T (n)y u Rigid-SmoothFIGURE 5-2Line of symmetry boundary condition.Formulation and Solution Strategies 85the normal traction becomes simply the stress component s y ,while the tangential traction reduces to t xy .For this case,s x exists only inside the region,and thus this component of stress cannot be speciﬁed on the boundary surface y ¼constant .A similar situation exists on the vertical boundary x ¼constant ,where the normal traction is now s x ,the tangential traction is t xy and the stress component s y exists inside the domain.Similar arguments can be made for polar coordinate boundary surfaces as shown.Drawing the appropriate element along the boundary as illustrated allows a clear visualization of the particular stress components that act on the surface in question.Such a sketch also allows determination of the positive directions of these boundary stresses,and this is useful to properly match with boundary loadings that might be prescribed.It is recommended that sketches similar to Figure 5-3be used to aid in the proper development of boundary conditions during problem formulation.Consider the pair of two-dimensional example problems with mixed conditions as shown in Figure 5-4.For the rectangular plate problem,all four boundaries are coordinate surfaces,andr(Cartesian Coordinate Boundaries)(Polar Coordinate Boundaries)xFIGURE 5-3Boundary stress components on coordinatesurfaces.(n)(n)(n)= σy = 0= t xy = 0(n)T x (n)T y = t xy =0,Traction Condition (Non-Coordinate Surface Boundary)(Coordinate Surface Boundaries)FIGURE 5-4Example boundary conditions.86FOUNDATIONS AND ELEMENTARY APPLICATIONSthis simpliﬁes speciﬁcation of particular boundary conditions.Theﬁxed conditions on the left edge simply require that x and y displacement components vanish on x¼0,and this speciﬁca-tion does not change even if this were not a coordinate surface.However,as per our previous discussion,the traction conditions on the other three boundaries simplify because they are coordinate surfaces.These simpliﬁcations are shown in theﬁgure for each of the traction speciﬁed surfaces.The second problem of a tapered cantilever beam has an inclined face that is not a coordinate surface.For this problem,theﬁxed end and top surface follow similar procedures as in theﬁrst example and are speciﬁed in theﬁgure.However,on the inclined face,the traction is to be zero and this does not reduce to a simple speciﬁcation of the vanishing of individual stress components.On this face each traction component is set to zero,giving the resultT(n)x¼s x n xþt xy n y¼0T(n)y¼t xy n xþs y n y¼0where n x and n y are the components of the unit normal vector to the inclined face.This is the more general type of speciﬁcation,and it should be clearly noted that none of the individual stress components in the x,y system will vanish along this surface.It should also be pointed out for this problem that the unit normal vector components are constants for all points on the inclined face.However,for curved boundaries the normal vector changes with surface position.Although these examples provide some background on typical boundary conditions,many other types will be encountered throughout the text.Several exercises at the end of this chapter provide additional examples that will prove to be useful for students new to the elasticity formulation.We are now in the position to formulate and classify the three fundamental boundary-value problems in the theory of elasticity that are related to solving the general system ofﬁeld equations(5.1.5).Our presentation is limited to the static case.Problem1:Traction problemDetermine the distribution of displacements,strains,and stresses in the interior of an elastic body in equilibrium when body forces are given and the distribution of the tractions are prescribed over the surface of the body,T(n)i(x(s)i)¼f i(x(s)i)(5:2:1) where x(s)i denotes boundary points and f i(x(s)i)are the prescribed traction values. Problem2:Displacement problemDetermine the distribution of displacements,strains,and stresses in the interior of an elastic body in equilibrium when body forces are given and the distribution of the displacements are prescribed over the surface of the body,u i(x(s)i)¼g i(x(s)i)(5:2:2) where x(s)i denotes boundary points and g i(x(s)i)are the prescribed displacement values.Formulation and Solution Strategies87Problem3:Mixed problemDetermine the distribution of displacements,strains,and stresses in the interior of an elasticbody in equilibrium when body forces are given and the distribution of the tractions areprescribed as per(5.2.1)over the surface S t and the distribution of the displacementsare prescribed as per(5.2.2)over the surface S u of the body(see Figure5-1).As mentioned previously,the solution to any of these types of problems is formidable,andfurther reduction and simpliﬁcation of(5.1.5)is required for the development of analyticalsolution methods.Based on the description of Problem1with only traction boundary condi-tions,it would appear to be desirable to express the fundamental system solely in terms ofstress,that is,J(t){s ij;l,m,F i}thereby reducing the number of unknowns in the system.Likewise for Problem2,a displacement-only formulation of the form J(u){u i;l,m,F i}wouldappear to simplify the problem.We now pursue such specialized formulations and explicitlydetermine these reducedﬁeld equation systems.5.3Stress FormulationFor theﬁrst fundamental problem in elasticity,the boundary conditions are to be given only interms of the tractions or stress components.In order to develop solution methods for this case,it is very helpful to reformulate the general system(5.1.5)by eliminating the displacementsand strains and thereby cast a new system solely in terms of the stresses.We now develop thisreformulated system.By eliminating the displacements,we must now include the compatibil-ity equations in the fundamental system ofﬁeld equations.We therefore start by using Hooke’slaw(5.1.4)2and eliminate the strains in the compatibility relations(5.1.2)to gets ij,kkþs kk,ijÀs ik,jkÀs jk,ik¼n 1þn (s mm,kk d ijþs mm,ij d kkÀs mm,jk d ikÀs mm,ik d jk)(5:3:1)where we have used the arguments of Section2.6,that the six meaningful compatibility relations are found by setting k¼l in(5.1.2).Although equations(5.3.1)represent the compatibility in terms of stress,a more useful result is found by incorporating the equilibrium equations into the system.Recall that from(5.1.3),s ij,j¼ÀF i,and also note that d kk¼3. Substituting these results into(5.3.1)givess ij,kkþ11þns kk,ij¼n1þns mm,kk d ijÀF i,jÀF j,i(5:3:2)For the case i¼j,relation(5.3.2)reduces to s ii,kk¼À1þn1ÀnF i,i.Substituting this result backinto(5.3.2)gives the desired relations ij,kkþ11þns kk,ij¼Àn1Ànd ij F k,kÀF i,jÀF j,i(5:3:3)This result is the compatibility relations in terms of the stress and is commonly called theBeltrami-Michell compatibility equations.For the case with no body forces,these relations canbe expressed as the following six scalar equations:88FOUNDATIONS AND ELEMENTARY APPLICATIONS(1þn)r2s xþ@2@x(s xþs yþs z)¼0(1þn)r2s yþ@2@y(s xþs yþs z)¼0(1þn)r2s zþ@2@z(s xþs yþs z)¼0(1þn)r2t xyþ@2@x@y(s xþs yþs z)¼0(1þn)r2t yzþ@2@y@z(s xþs yþs z)¼0(1þn)r2t zxþ@2(s xþs yþs z)¼0(5:3:4)Recall that the six developed relations(5.3.3)or(5.3.4)actually represent three independentresults as per our discussion in Section2.6.Thus,combining these results with the threeequilibrium equations(5.1.3)provides the necessary six relations to solve for the six unknownstress components for the general three-dimensional case.This system constitutes the stressformulation for elasticity theory and is appropriate for use with traction boundary conditionproblems.Once the stresses have been determined,the strains may be found from Hooke’s law(5.1.4)z,and the displacements can be then be computed through integration of(5.1.1).As perour previous discussion in Section2.2,such an integration process determines the displace-ments only up to an arbitrary rigid-body motion,and the displacements obtained are single-valued only if the region under study is simply connected.The system of equations for the stress formulation is still rather complex,and analytical solutions are commonly determined for this case by making use of stress functions.Thisconcept establishes a representation for the stresses that automatically satisﬁes the equilibriumequations.For the two-dimensional case,this concept represents the in-plane stresses in termsof a single function.The representation satisﬁes equilibrium,and the remaining compatibilityequations yield a single partial differential equation(biharmonic equation)in terms of thestress function.Having reduced the system to a single equation then allows us to employ manyanalytical methods toﬁnd solutions of interest.Further discussion on these techniques ispresented in subsequent chapters.5.4Displacement FormulationWe now wish to develop the reduced set ofﬁeld equations solely in terms of the displacements.This system is referred to as the displacement formulation and is most useful when combinedwith displacement-only boundary conditions found in the Problem2statement.This develop-ment is somewhat more straightforward than our previous discussion for the stress formulation.For this case,we wish to eliminate the strains and stresses from the fundamental system(5.1.5).This is easily accomplished by using the strain-displacement relations in Hooke’s lawto gives ij¼l u k,k d ijþm(u i,jþu j,i)(5:4:1) which can be expressed as six scalar equationsFormulation and Solution Strategies89s x¼l@u@xþ@n@yþ@w@zþ2m@u@xs y¼l@u@xþ@v@yþ@w@zþ2m@v@ys z¼l@u@xþ@v@yþ@w@zþ2m@w@zt xy¼m @u@yþ@v@x,t yz¼m @v@zþ@w@y,t zx¼m @w@xþ@u@z(5:4:2)Using these relations in the equilibrium equations gives the resultm u i,kkþ(lþm)u k,kiþF i¼0(5:4:3)which are the equilibrium equations in terms of the displacements and are referred to as Navier’s or Lame´’s equations.This system can be expressed in vector form asm r2uþ(lþm),(,Áu)þF¼0(5:4:4) or written out in terms of the three scalar equationsm r2uþ(lþm)@@x@u@xþ@v@yþ@w@zþF x¼0m r2vþ(lþm)@@y@u@xþ@v@yþ@w@zþF y¼0m r2wþ(lþm)@@z@u@xþ@v@yþ@w@zþF z¼0(5:4:5)where the Laplacian is given by r2¼(@2=@x2)þ(@2=@y2)þ(@2=@z2).Navier’s equations arethe desired formulation for the displacement problem,and the system represents three equa-tions for the three unknown displacement components.Similar to the stress formulation,thissystem is still difﬁcult to solve,and additional mathematical techniques have been developedto further simplify these equations for problem mon methods normally employthe use of displacement potential functions.It is shown in Chapter13that several suchschemes can be developed that allow the displacement vector to be expressed in terms ofparticular potentials.These schemes generally simplify the problem by yielding uncoupledgoverning equations in terms of the displacement potentials.This then allows several analyt-ical methods to be employed to solve problems of interest.Several of these techniques arediscussed in later sections of the text.To help acquire a general understanding of these results,a summaryﬂow chart of the developed stress and displacement formulations is shown in Figure5-5.Note that for thestress formulation,the resulting system J(t){s ij;l,m,F i}is actually dependent onlyon the single material constant Poisson’s ratio,and thus it could be expressed asJ(t){s ij;n,F i}.90FOUNDATIONS AND ELEMENTARY APPLICATIONS5.5Principle of SuperpositionA very useful tool for the solution to many problems in engineering science is the principle ofsuperposition.This technique applies to any problem that is governed by linear equations.Under the assumption of small deformations and linear elastic constitutive behavior,allelasticityﬁeld equations(see Figure5-5)are linear.Furthermore,the usual boundary condi-tions speciﬁed by relations(5.2.1)and(5.2.2)are also linear.Thus,under these conditions allgoverning equations are linear,and the superposition concept can be applied.It can be easilyproved(see Chou and Pagano1967)that the general statement of the principle can beexpressed as follows:Principle of Superposition:For a given problem domain,if the state{s(1)ij,e(1)ij,u(1)i}is asolution to the fundamental elasticity equations with prescribed body forces F(1)i andsurface tractions T(1)i,and the state{s(2)ij,e(2)ij,u(2)i}is a solution to the fundamentalequations with prescribed body forces F(2)i and surface tractions T(2)i,then the sta-te{s(1)ijþs(2)ij,e(1)ijþe(2)ij,u(1)iþu(2)i}will be a solution to the problem with body forcesF(1)iþF(2)i and surface tractions T(1)iþT(2)i.In order to see a more direct application of this principle,consider a simple two-dimensionalcase with no body forces as shown in Figure5-6.It can be observed that the solution to themore complicated biaxial loading case(1)þ(2)is thus found by adding the two simplerproblems.This is a common use of the superposition principle,and we make repeated use ofthis application throughout the text.Thus,once the solutions to some simple problems aregenerated,we can combine these results to generate a solution to a more complicated case withsimilar geometry.Formulation and Solution Strategies915.6Saint-Venant’s PrincipleConsider the set of three identical rectangular strips under compressive loadings as shown in Figure 5-7.As indicated,the only difference between each problem is the loading.Because the total resultant load applied to each problem is identical (statically equivalent loadings),it is expected that the resulting stress,strain,and displacement ﬁelds near the bottom of each strip would be approximately the same.This behavior can be generalized by considering an elastic solid with an arbitrary loading T (n )over a boundary portion S *,as shown in Figure 5-8.Based on experience from other ﬁeld problems in engineering science,it seems logical that the particular boundary loading would produce detailed and characteristic effects only in the vicinity of S *.In other words,we expect that at points far away from S *the stresses generally depend more on the resultant F R of the tractions rather than on the exact distribution.Thus,the characteristic signature of the generated stress,strain,and displacement ﬁelds from a given boundary loading tend to disappear as we move away from the boundary loading points.These concepts form the principle of Saint-Venant ,which can be stated as follows:Saint-Venant’s Principle:The stress,strain,and displacement ﬁelds caused by two different statically equivalent force distributions on parts of the body far away from the loading points are approximately the same.=+ ij ij i ij ij i {s (1)ij + s (2) ij , e (1) ij , + e (2) ij , u (1) i + u (2) i }FIGURE 5-6Two-dimensional superposition example.2P 2P3P 3P 3P FIGURE 5-7Statically equivalent loadings.92FOUNDATIONS AND ELEMENTARY APPLICATIONSThis statement of the principle includes qualitative terms such as far away and approxi-mately the same ,and thus does not provide quantitative estimates of the differences between the two elastic ﬁelds in question.Quantitative results have been developed by von Mises (1945),Sternberg (1954),and Toupin (1965),while Horgan (1989)has presented a recent review of related work.Some of this work is summarized in Boresi and Chong (2000).If we restrict our solution to points away from the boundary loading,Saint-Venant’s principle allows us to change the given boundary conditions to a simpler statically equivalent statement and not affect the resulting solution.Such a simpliﬁcation of the boundary conditions greatly increases our chances of ﬁnding an analytical solution to the problem.This concept therefore proves to be very useful,and we formally outline this solution scheme in the next section.5.7General Solution StrategiesHaving completed our formulation and related solution principles,we now wish to examine some general solution strategies commonly used to solve elasticity problems.At this stage we categorize particular methods and outline only typical techniques that are commonly used.As we move further along in the text,many of these methods are developed in detail and are applied in speciﬁc problem solution.We ﬁrst distinguish three general methods of solution called direct,inverse ,and semi-inverse .5.7.1Direct MethodThis method seeks to determine the solution by direct integration of the ﬁeld equations (5.1.5)or equivalently the stress and/or displacement formulations given in Figure 5-5.Boundary conditions are to be satisﬁed exactly.This method normally encounters signiﬁcant mathemat-ical difﬁculties,thus limiting its application to problems with simple geometry.EXAMPLE 5-1:Direct Integration Example:Stretching of Prismatic Bar Under Its Own WeightAs an example of a simple direct integration problem,consider the case of a uniform prismatic bar stretched by its own weight,as shown in Figure 5-9.The body forces forContinuedFIGURE 5-8Saint-Venant concept.5.7.2Inverse MethodFor this technique,particular displacements or stresses are selected that satisfy the basicﬁeld equations.A search is then conducted to identify a speciﬁc problem that would be solved by this solutionﬁeld.This amounts to determine appropriate problem geometry,boundary conditions,and body forces that would enable the solution to satisfy all conditions on the ing this scheme it is sometimes difﬁcult to construct solutions to a speciﬁc problem of practical interest.5.7.3Semi-Inverse MethodIn this scheme part of the displacement and/or stressﬁeld is speciﬁed,and the other remaining portion is determined by the fundamentalﬁeld equations(normally using direct integration) and the boundary conditions.It is often the case that constructing appropriate displacement and/or stress solutionﬁelds can be guided by approximate strength of materials theory.The usefulness of this approach is greatly enhanced by employing Saint-Venant’s principle, whereby a complicated boundary condition can be replaced by a simpler statically equivalent distribution.EXAMPLE5-3:Semi-Inverse Example:Torsion of Prismatic BarsA simple semi-inverse example may be borrowed from the torsion problem that isdiscussed in detail in Chapter9.Skipping for now the developmental details,we propose the following displacementﬁeld:ContinuedThere are numerous mathematical techniques used to solve the elasticity ﬁeld equations.Many techniques involve the development of exact analytical solutions ,while others involve the construction of approximate solution schemes .A third procedure involves the establishment of numerical solution methods .We now brieﬂy provide an overview of each of these techniques.5.7.4Analytical Solution ProceduresA variety of analytical solution methods are used to solve the elasticity ﬁeld equations.The following sections outline some of the more common methods.Power Series MethodFor many two-dimensional elasticity problems,the stress formulation leads to the use of a stress function f (x ,y ).It is shown that the entire set of ﬁeld equations reduces to a single partial differential equation (biharmonic equation)in terms of this stress function.A general mathematical scheme to solve this equation is to look for solutions in terms of a power series in the independent variables,that is,f (x ,y )¼P C mn x m y n (see Neou 1957).Use of the boundary conditions determines the coefﬁcients and number of terms to be used in the series.This method is employed to develop two-dimensional solutions in Section 8.1.Fourier MethodA general scheme to solve a large variety of elasticity problems employs the Fourier method.This procedure is normally applied to the governing partial differential equations by using。

自动化专业英语词汇大全

自动化专业英语词汇大全accelerationtransducer加速度传感器basecoordinatesystem基座坐标系acceptancetesting验收测试Bayesclassifier贝叶斯分类器accessibility可及性bearingalignment方位对准accumulatederror累积误差bellowspressuregauge波纹管压力表AC-DC-ACfrequencyconverter交-直-交变频器benefit-costanalysis收益本钱分析AC(alternatingcurrent)electricdrive交流电子传bilinearsystem双线性系统动biocybernetics生物控制论activeattitudestabilization主动姿态稳定biologicalfeedbacksystem生物反响系统actuator驱动器，执行机构blackboxtestingapproach黑箱测试法adaline线性适应元blindsearch盲目搜索adaptationlayer适应层blockdiagonalization块对角化adaptivetelemetersystem适应遥测系统Boltzmanmachine玻耳兹曼机adjointoperator伴随算子bottom-updevelopment自下而上开发admissibleerror容许误差boundaryvalueanalysis边界值分析aggregationmatrix集结矩阵brainstormingmethod头脑风暴法AHP(analytichierarchyprocess)层次分析法breadth-firstsearch广度优先搜索amplifyingelement放大环节butterflyvalve蝶阀analog-digitalconversion模数转换CAE(computeraidedengineering)计算机辅助工annunciator信号器程antennapointingcontrol天线指向控制CAM(computeraidedmanufacturing)计算机辅助anti-integralwindup抗积分饱卷制造aperiodicdecomposition非周期分解Camflexvalve偏心旋转阀aposterioriestimate后验估计canonicalstatevariable标准化状态变量approximatereasoning近似推理capacitivedisplacementtransducer电容式位移传aprioriestimate先验估计感器articulatedrobot关节型机器人capsulepressuregauge膜盒压力表assignmentproblem配置问题，分配问题CARD计算机辅助研究开发associativememorymodel联想记忆模型Cartesianrobot直角坐标型机器人associatron联想机cascadecompensation串联补偿asymptoticstability渐进稳定性catastrophetheory突变论attainedposedrift实际位姿漂移centrality集中性attitudeacquisition姿态捕获chainedaggregation链式集结AOCS(attritudeandorbitcontrolsystem)姿态轨chaos混沌道控制系统characteristiclocus特征轨迹attitudeangularvelocity姿态角速度chemicalpropulsion化学推进attitudedisturbance姿态扰动calrity清晰性attitudemaneuver姿态机动classicalinformationpattern经典信息模式attractor吸引子classifier分类器augmentability可扩大性clinicalcontrolsystem临床控制系统augmentedsystem增广系统closedlooppole闭环极点automaticmanualstation自动-手动操作器closedlooptransferfunction闭环传递函数automaton自动机clusteranalysis聚类分析autonomoussystem自治系统coarse-finecontrol粗-精控制backlashcharacteristics间隙特性cobwebmodel蛛网模型coefficientmatrix系数矩阵costatevariable共态变量cognitivescience认知科学cost-effectivenessanalysis费用效益分析cognitron认知机couplingoforbitandattitude轨道和姿态耦合coherentsystem单调关联系统criticaldamping临界阻尼combinationdecision组合决策criticalstability临界稳定性combinatorialexplosion组合爆炸cross-overfrequency穿越频率，交越频率combinedpressureandvacuumgauge压力真空currentsourceinverter电流[源]型逆变器表cut-offfrequency截止频率commandpose指令位姿cybernetics控制论companionmatrix相伴矩阵cyclicremotecontrol循环遥控compartmentalmodel房室模型cylindricalrobot圆柱坐标型机器人compatibility相容性，兼容性dampedoscillation阻尼振荡compensatingnetwork补偿网络damper阻尼器compensation补偿，矫正dampingratio阻尼比compliance柔顺，顺应dataacquisition数据采集compositecontrol组合控制dataencryption数据加密computablegeneralequilibriummodel可计算一datapreprocessing数据预处理般均衡模型dataprocessor数据处理器conditionallyinstability条件不稳定性DCgenerator-motorsetdrive直流发电机-电动机configuration组态组传动connectionism连接机制Dcontroller微分控制器connectivity连接性decentrality分散性conservativesystem守恒系统decentralizedstochasticcontrol分散随机控制consistency一致性decisionspace决策空间constraintcondition约束条件decisionsupportsystem决策支持系统consumptionfunction消费函数decomposition-aggregationapproach分解集结法context-freegrammar上下文无关语法decouplingparameter解耦参数continuousdiscreteeventhybridsystemdeductive-inductivehybridmodelingmethod演simulation连续离散事件混合系统仿真绎与归纳混合建模法continuousduty连续工作制delayedtelemetry延时遥测controlaccuracy控制精度derivationtree导出树controlcabinet控制柜derivativefeedback微分反响controllabilityindex可控指数describingfunction描述函数controllablecanonicalform可控标准型desiredvalue希望值[control]plant控制对象,被控对象despinner消旋体controllinginstrument控制仪表destination目的站controlmomentgyro控制力矩陀螺detector检出器controlpanel控制屏，控制盘deterministicautomaton确定性自动机controlsynchro控制[式]自整角机deviation偏差controlsystemsynthesis控制系统综合deviationalarm偏差报警器controltimehorizon控制时程DFD数据流图cooperativegame合作对策diagnosticmodel诊断模型coordinabilitycondition可协调条件diagonallydominantmatrix对角主导矩阵coordinationstrategy协调策略diaphragmpressuregauge膜片压力表coordinator协调器differenceequationmodel差分方程模型cornerfrequency转折频率differentialdynamicalsystem微分动力学系统----differentialgame微分对策economicindicator经济指标differentialpressurelevelmeter差压液位计eddycurrentthicknessmeter电涡流厚度计differentialpressuretransmitter差压变送器effectiveness有效性differentialtransformerdisplacementtransducereffectivenesstheory效益理论差动变压器式位移传感器elasticityofdemand需求弹性differentiationelement微分环节electricactuator电动执行机构digitalfiler数字滤波器electricconductancelevelmeter电导液位计digitalsignalprocessing数字信号处理electricdrivecontrolgear电动传动控制设备digitization数字化electrichydraulicconverter电-液转换器digitizer数字化仪electricpneumaticconverter电-气转换器dimensiontransducer尺度传感器electrohydraulicservovale电液伺服阀directcoordination直接协调electromagneticflowtransducer电磁流量传感器disaggregation解裂electronicbatchingscale电子配料秤discoordination失协调electronicbeltconveyorscale电子皮带秤discreteeventdynamicsystem离散事件动态系统electronichopperscale电子料斗秤discretesystemsimulationlanguage离散系统仿elevation仰角真语言emergencystop异常停顿discriminantfunction判别函数empiricaldistribution经历分布displacementvibrationamplitudetransducer位endogenousvariable内生变量移振幅传感器equilibriumgrowth均衡增长dissipativestructure耗散构造equilibriumpoint平衡点distributedparametercontrolsystem分布参数控equivalencepartitioning等价类划分制系统ergonomics工效学distrubance扰动error误差disturbancecompensation扰动补偿error-correctionparsing纠错剖析diversity多样性estimate估计量divisibility可分性estimationtheory估计理论domainknowledge领域知识evaluationtechnique评价技术dominantpole主导极点eventchain事件链dose-responsemodel剂量反响模型evolutionarysystem进化系统dualmodulationtelemeteringsystem双重调制遥exogenousvariable外生变量测系统expectedcharacteristics希望特性dualprinciple对偶原理externaldisturbance外扰dualspinstabilization双自旋稳定factbase事实dutyratio负载比failurediagnosis故障诊断dynamicbraking能耗制动fastmode快变模态dynamiccharacteristics动态特性feasibilitystudy可行性研究dynamicdeviation动态偏差feasiblecoordination可行协调dynamicerrorcoefficient动态误差系数feasibleregion可行域dynamicexactness动它吻合性featuredetection特征检测dynamicinput-outputmodel动态投入产出模型featureextraction特征抽取econometricmodel计量经济模型feedbackcompensation反响补偿economiccybernetics经济控制论feedforwardpath前馈通路economiceffectiveness经济效益fieldbus现场总线economicevaluation经济评价finiteautomaton有限自动机economicindex经济指数FIP(factoryinformationprotocol)工厂信息协议firstorderpredicatelogic一阶谓词逻辑harmoniousstrategy和谐策略fixedsequencemanipulator固定顺序机械手heuristicinference启发式推理fixedsetpointcontrol定值控制hiddenoscillation隐蔽振荡FMS(flexiblemanufacturingsystem)柔性制造系hierarchicalchart层次构造图统hierarchicalplanning递阶规划flowsensor/transducer流量传感器hierarchicalcontrol递阶控制flowtransmitter流量变送器homeostasis内稳态fluctuation涨落homomorphicmodel同态系统forcedoscillation强迫振荡horizontaldecomposition横向分解formallanguagetheory形式语言理论hormonalcontrol内分泌控制formalneuron形式神经元hydraulicstepmotor液压步进马达forwardpath正向通路hypercycletheory超循环理论forwardreasoning正向推理Icontroller积分控制器fractal分形体，分维体identifiability可辨识性frequencyconverter变频器IDSS(intelligentdecisionsupportsystem)智能frequencydomainmodelreductionmethod频域决策支持系统模型降阶法imagerecognition图像识别frequencyresponse频域响应impulse冲量fullorderobserver全阶观测器impulsefunction冲击函数，脉冲函数functionaldecomposition功能分解inching点动FES(functionalelectricalstimulation)功能电刺激incompatibilityprinciple不相容原理functionalsimularity功能相似incrementalmotioncontrol增量运动控制fuzzylogic模糊逻辑indexofmerit品质因数gametree对策树inductiveforcetransducer电感式位移传感器gatevalve闸阀inductivemodelingmethod归纳建模法generalequilibriumtheory一般均衡理论industrialautomation工业自动化generalizedleastsquaresestimation广义最小二inertialattitudesensor惯性姿态敏感器乘估计inertialcoordinatesystem惯性坐标系generationfunction生成函数inertialwheel惯性轮geomagnetictorque地磁力矩inferenceengine推理机geometricsimilarity几何相似infinitedimensionalsystem无穷维系统gimbaledwheel框架轮informationacquisition信息采集globalasymptoticstability全局渐进稳定性infraredgasanalyzer红外线气体分析器globaloptimum全局最优inherentnonlinearity固有非线性globevalve球形阀inherentregulation固有调节goalcoordinationmethod目标协调法initialdeviation初始偏差grammaticalinference文法推断initiator发起站graphicsearch图搜索injectionattitude入轨姿势gravitygradienttorque重力梯度力矩input-outputmodel投入产出模型grouptechnology成组技术instability不稳定性guidancesystem制导系统instructionlevellanguage指令级语言gyrodriftrate陀螺漂移率integralofabsolutevalueoferrorcriterion绝对gyrostat陀螺体误差积分准那么Halldisplacementtransducer霍尔式位移传感器integralofsquarederrorcriterion平方误差积分准hardware-in-the-loopsimulation半实物仿真那么harmoniousdeviation和谐偏差integralperformancecriterion积分性能准那么integrationinstrument积算仪器localasymptoticstability局部渐近稳定性integrity整体性localoptimum局部最优intelligentterminal智能终端logmagnitude-phasediagram对数幅相图interactedsystem互联系统，关联系统longtermmemory长期记忆interactivepredictionapproach互联预估法，关联lumpedparametermodel集总参数模型预估法Lyapunovtheoremofasymptoticstability李雅普interconnection互联诺夫渐近稳定性定理intermittentduty断续工作制macro-economicsystem宏观经济系统internaldisturbance内扰magneticdumping磁卸载ISM(interpretivestructuremodeling)解释构造建magnetoelasticweighingcell磁致弹性称重传感器模法magnitude-frequencycharacteristic幅频特性invariantembeddingprinciple不变嵌入原理magnitudemargin幅值裕度inventorytheory库伦论magnitudescalefactor幅值比例尺inverseNyquistdiagram逆奈奎斯特图manipulator机械手inverter逆变器man-machinecoordination人机协调investmentdecision投资决策manualstation手动操作器isomorphicmodel同构模型MAP(manufacturingautomationprotocol)制造iterativecoordination迭代协调自动化协议jetpropulsion喷气推进marginaleffectiveness边际效益job-lotcontrol分批控制Mason'sgainformula梅森增益公式joint关节masterstation主站Kalman-Bucyfiler卡尔曼-布西滤波器matchingcriterion匹配准那么knowledgeaccomodation知识顺应maximumlikelihoodestimation最大似然估计knowledgeacquisition知识获取maximumovershoot最大超调量knowledgeassimilation知识同化maximumprinciple极大值原理KBMS(knowledgebasemanagementsystem)知mean-squareerrorcriterion均方误差准那么识库管理系统mechanismmodel机理模型knowledgerepresentation知识表达meta-knowledge元知识ladderdiagram梯形图metallurgicalautomation冶金自动化lag-leadcompensation滞后超前补偿minimalrealization最小实现Lagrangeduality拉格朗日对偶性minimumphasesystem最小相位系统Laplacetransform拉普拉斯变换minimumvarianceestimation最小方差估计largescalesystem大系统minorloop副回路lateralinhibitionnetwork侧抑制网络missile-targetrelativemovementsimulator弹体leastcostinput最小本钱投入-目标相对运动仿真器leastsquarescriterion最小二乘准那么modalaggregation模态集结levelswitch物位开关modaltransformation模态变换librationdamping天平动阻尼MB(modelbase)模型库limitcycle极限环modelconfidence模型置信度linearizationtechnique线性化方法modelfidelity模型逼真度linearmotionelectricdrive直线运动电气传动modelreferenceadaptivecontrolsystem模型参linearmotionvalve直行程阀考适应控制系统linearprogramming线性规划modelverification模型验证LQR(linearquadraticregulatorproblem)线性二modularization模块化次调节器问题MEC(mosteconomiccontrol)最经济控制loadcell称重传感器motionspace可动空间MTBF(meantimebetweenfailures)平均故障间隔orderparameter序参数时间orientationcontrol定向控制MTTF(meantimetofailures)平均无故障时间originator始发站multi-attributiveutilityfunction多属性效用函数oscillatingperiod振荡周期multicriteria多重判据outputpredictionmethod输出预估法multilevelhierarchicalstructure多级递阶构造ovalwheelflowmeter椭圆齿轮流量计multiloopcontrol多回路控制overalldesign总体设计multi-objectivedecision多目标决策overdamping过阻尼multistatelogic多态逻辑overlappingdecomposition交叠分解multistratumhierarchicalcontrol多段递阶控制Padeapproximation帕德近似multivariablecontrolsystem多变量控制系统Paretooptimality帕雷托最优性myoelectriccontrol肌电控制passiveattitudestabilization被动姿态稳定Nashoptimality纳什最优性pathrepeatability路径可重复性naturallanguagegeneration自然语言生成patternprimitive模式基元nearest-neighbor最近邻PR(patternrecognition)模式识别necessitymeasure必然性侧度Pcontrol比例控制器negativefeedback负反响peaktime峰值时间neuralassembly神经集合penaltyfunctionmethod罚函数法neuralnetworkcomputer神经网络计算机perceptron感知器Nicholschart尼科尔斯图periodicduty周期工作制noeticscience思维科学perturbationtheory摄动理论noncoherentsystem非单调关联系统pessimisticvalue悲观值noncooperativegame非合作博弈phaselocus相轨迹nonequilibriumstate非平衡态phasetrajectory相轨迹nonlinearelement非线性环节phaselead相位超前nonmonotoniclogic非单调逻辑photoelectrictachometrictransducer光电式转速nonparametrictraining非参数训练传感器nonreversibleelectricdrive不可逆电气传动phrase-structuregrammar短句构造文法nonsingularperturbation非奇异摄动physicalsymbolsystem物理符号系统non-stationaryrandomprocess非平稳随机过程piezoelectricforcetransducer压电式力传感器nuclearradiationlevelmeter核辐射物位计playbackrobot示教再现式机器人nutationsensor章动敏感器PLC(programmablelogiccontroller)可编程序逻Nyquiststabilitycriterion奈奎斯特稳定判据辑控制器objectivefunction目标函数plugbraking反接制动observabilityindex可观测指数plugvalve旋塞阀observablecanonicalform可观测标准型pneumaticactuator气动执行机构on-lineassistance在线帮助point-to-pointcontrol点位控制on-offcontrol通断控制polarrobot极坐标型机器人openlooppole开环极点poleassignment极点配置operationalresearchmodel运筹学模型pole-zerocancellation零极点相消opticfibertachometer光纤式转速表polynomialinput多项式输入optimaltrajectory最优轨迹portfoliotheory投资搭配理论optimizationtechnique最优化技术poseovershoot位姿过调量orbitalrendezvous轨道交会positionmeasuringinstrument位置测量仪orbitgyrocompass轨道陀螺罗盘posentiometricdisplacementtransducer电位器orbitperturbation轨道摄动式位移传感器positivefeedback正反响realizability可实现性,能实现性powersystemautomation电力系统自动化realtimetelemetry实时遥测predicatelogic谓词逻辑receptivefield感受野pressuregaugewithelectriccontact电接点压力表rectangularrobot直角坐标型机器人pressuretransmitter压力变送器rectifier整流器pricecoordination价格协调recursiveestimation递推估计primalcoordination主协调reducedorderobserver降阶观测器primaryfrequencyzone主频区redundantinformation冗余信息PCA(principalcomponentanalysis)主成分分析法reentrycontrol再入控制principleofturnpike大道原理regenerativebraking回馈制动，再生制动priority优先级regionalplanningmodel区域规划模型process-orientedsimulation面向过程的仿真regulatingdevice调节装载productionbudget生产预算regulation调节productionrule产生式规那么relationalalgebra关系代数profitforecast利润预测relaycharacteristic继电器特性remotemanipulator遥控操作器PERT(programevaluationandreviewtechnique)方案评审技术remoteregulating遥调programsetstation程序设定操作器remotesetpointadjuster远程设定点调整器proportionalcontrol比例控制rendezvousanddocking交会和对接proportionalplusderivativecontroller比例微分控reproducibility再现性制器resistancethermometersensor热电阻protocolengineering协议工程resolutionprinciple归结原理prototype原型resourceallocation资源分配pseudorandomsequence伪随机序列responsecurve响应曲线pseudo-rate-incrementcontrol伪速率增量控制returndifferencematrix回差矩阵pulseduration脉冲持续时间returnratiomatrix回比矩阵pulsefrequencymodulationcontrolsystem脉冲reverberation回响调频控制系统reversibleelectricdrive可逆电气传动pulsewidthmodulationcontrolsystem脉冲调宽revoluterobot关节型机器人控制系统revolutionspeedtransducer转速传感器PWMinverter脉宽调制逆变器rewritingrule重写规那么pushdownautomaton下推自动机rigidspacecraftdynamics刚性航天动力学QC(qualitycontrol)质量管理riskdecision风险分析quadraticperformanceindex二次型性能指标robotics机器人学qualitativephysicalmodel定性物理模型robotprogramminglanguage机器人编程语言quantizednoise量化噪声robustcontrol鲁棒控制quasilinearcharacteristics准线性特性robustness鲁棒性queuingtheory排队论rollgapmeasuringinstrument辊缝测量仪radiofrequencysensor射频敏感器rootlocus根轨迹rampfunction斜坡函数rootsflowmeter腰轮流量计randomdisturbance随机扰动rotameter浮子流量计，转子流量计randomprocess随机过程rotaryeccentricplugvalve偏心旋转阀rateintegratinggyro速率积分陀螺rotarymotionvalve角行程阀ratiostation比值操作器rotatingtransformer旋转变压器reachability可达性Routhapproximationmethod劳思近似判据reactionwheelcontrol反作用轮控制routingproblem路径问题sampled-datacontrolsystem采样控制系统socioeconomicsystem社会经济系统samplingcontrolsystem采样控制系统softwarepsychology软件心理学saturationcharacteristics饱和特性solararraypointingcontrol太阳帆板指向控制scalarLyapunovfunction标量李雅普诺夫函数solenoidvalve电磁阀SCARA(selectivecomplianceassemblyrobotsource源点arm)平面关节型机器人specificimpulse比冲scenarioanalysismethod情景分析法speedcontrolsystem调速系统sceneanalysis物景分析spinaxis自旋轴s-domains域spinner自旋体self-operatedcontroller自力式控制器stabilitycriterion稳定性判据self-organizingsystem自组织系统stabilitylimit稳定极限self-reproducingsystem自繁殖系统stabilization镇定，稳定self-tuningcontrol自校正控制Stackelbergdecisiontheory施塔克尔贝格决策理论semanticnetwork语义网络stateequationmodel状态方程模型semi-physicalsimulation半实物仿真statespacedescription状态空间描述sensingelement敏感元件staticcharacteristicscurve静态特性曲线sensitivityanalysis灵敏度分析stationaccuracy定点精度sensorycontrol感觉控制stationaryrandomprocess平稳随机过程sequentialdecomposition顺序分解statisticalanalysis统计分析sequentialleastsquaresestimation序贯最小二乘statisticpatternrecognition统计模式识别估计steadystatedeviation稳态偏差servocontrol伺服控制，随动控制steadystateerrorcoefficient稳态误差系数servomotor伺服马达step-by-stepcontrol步进控制settlingtime过渡时间stepfunction阶跃函数sextant六分仪stepwiserefinement逐步精化shorttermplanning短期方案stochasticfiniteautomaton随机有限自动机shorttimehorizoncoordination短时程协调straingaugeloadcell应变式称重传感器signaldetectionandestimation信号检测和估计strategicfunction策略函数signalreconstruction信号重构stronglycoupledsystem强耦合系统similarity相似性subjectiveprobability主观频率simulatedinterrupt仿真中断suboptimality次优性simulationblockdiagram仿真框图supervisedtraining监视学习simulationexperiment仿真实验supervisorycomputercontrolsystem计算机监控simulationvelocity仿真速度系统simulator仿真器sustainedoscillation自持振荡singleaxletable单轴转台swirlmeter旋进流量计singledegreeoffreedomgyro单自由度陀螺switchingpoint切换点singlelevelprocess单级过程symbolicprocessing符号处理singlevaluenonlinearity单值非线性synapticplasticity突触可塑性singularattractor奇异吸引子synergetics协同学singularperturbation奇异摄动syntacticanalysis句法分析sink汇点systemassessment系统评价slavedsystem受役系统systematology系统学slower-than-real-timesimulation欠实时仿真systemhomomorphism系统同态slowsubsystem慢变子系统systemisomorphism系统同构socio-cybernetics社会控制论systemengineering系统工程----tachometer转速表turbineflowmeter涡轮流量计targetflowtransmitter靶式流量变送器Turingmachine图灵机taskcycle作业周期two-timescalesystem双时标系统teachingprogramming示教编程ultrasoniclevelmeter超声物位计telemechanics远动学unadjustablespeedelectricdrive非调速电气传动telemeteringsystemoffrequencydivisiontypeunbiasedestimation无偏估计频分遥测系统underdamping欠阻尼telemetry遥测uniformlyasymptoticstability一致渐近稳定性teleologicalsystem目的系统uninterruptedduty不连续工作制，长期工作制teleology目的论unitcircle单位圆temperaturetransducer温度传感器unittesting单元测试templatebase模版库unsupervisedlearing非监视学习tensiometer X力计upperlevelproblem上级问题texture纹理urbanplanning城市规划theoremproving定理证明utilityfunction效用函数therapymodel治疗模型valueengineering价值工程thermocouple热电偶variablegain可变增益，可变放大系数thermometer温度计variablestructurecontrolsystem变构造控制thicknessmeter厚度计vectorLyapunovfunction向量李雅普诺夫函数three-axisattitudestabilization三轴姿态稳定velocityerrorcoefficient速度误差系数threestatecontroller三位控制器velocitytransducer速度传感器thrustvectorcontrolsystem推力矢量控制系统verticaldecomposition纵向分解thruster推力器vibratingwireforcetransducer振弦式力传感器timeconstant时间常数vibrometer振动计time-invariantsystem定常系统，非时变系统viscousdamping粘性阻尼timeschedulecontroller时序控制器voltagesourceinverter电压源型逆变器time-sharingcontrol分时控制vortexprecessionflowmeter旋进流量计time-varyingparameter时变参数vortexsheddingflowmeter涡街流量计top-downtesting自上而下测试WB(waybase)方法库topologicalstructure拓扑构造weighingcell称重传感器TQC(totalqualitycontrol)全面质量管理weightingfactor权因子trackingerror跟踪误差weightingmethod加权法trade-offanalysis权衡分析Whittaker-Shannonsamplingtheorem惠特克-香transferfunctionmatrix传递函数矩阵农采样定理transformationgrammar转换文法Wienerfiltering维纳滤波transientdeviation瞬态偏差workstationforcomputeraideddesign计算机辅助设计工作站transientprocess过渡过程transitiondiagram转移图w-planew平面transmissiblepressuregauge电远传压力表zero-basedbudget零基预算transmitter变送器zero-inputresponse零输入响应trendanalysis趋势分析zero-stateresponse零状态响应triplemodulationtelemeteringsystem三重调制zerosumgamemodel零和对策模型遥测系统z-transformz变换。

On the stochastic pendulum with Ornstein-Uhlenbeck noise

a r X i v :c o n d -m a t /0407198v 1 [c o n d -m a t .s t a t -m e c h ] 8 J u l 2004On the stochastic pendulum with Ornstein-Uhlenbeck noiseKirone MallickService de Physique Th´e orique,Centre d’´Etudes de Saclay,91191Gif-sur-Yvette Cedex,France ∗Philippe MarcqInstitut de Recherche sur les Ph´e nom`e nes Hors ´Equilibre,Universit´e de Provence,49rue Joliot-Curie,BP 146,13384Marseille Cedex 13,France †(Dated:March 16,2004)We study a frictionless pendulum subject to multiplicative random noise.Because of destructive interference between the angular displacement of the system and the noise term,the energy ﬂuctua-tions are reduced when the noise has a non-zero correlation time.We derive the long time behavior of the pendulum in the case of Ornstein-Uhlenbeck noise by a recursive adiabatic elimination pro-cedure.An analytical expression for the asymptotic probability distribution function of the energy is obtained and the results agree with numerical stly,we compare our method to other approximation schemes.PACS numbers:05.10.Gg,05.40.-a,05.45.-aI.INTRODUCTIONThe behavior of a nonlinear dynamical system is strongly modiﬁed when randomness is taken into account:noise can shift bifurcation thresholds,create new phases (noise-induced transitions),or even generate spatial patterns [1,2,3,4,5].The interplay of noise with nonlinearity gives rise to a variety of phenomena that constantly motivate new research of theoretical and practical signiﬁcance.Stochastic resonance [6]and biomolecular Brownian motors [7,8]are celebrated examples of nonlinear random systems of current interest.In particular,stochastic ratchets have generated a renewed interest in the study of simple mechanical systems subject to random interactions,the common ancestor of such models being Langevin’s description of Brownian motion.Many unexpected phenomena appear when one generalizes Langevin’s equations to include,e.g.,inertial terms,nonlinearities,external (multiplicative)noise or noise with ﬁnite correlation time;each of these new features opens a ﬁeld of investigations that calls for speciﬁc techniques or approximation schemes [9].Nonlinear oscillators with parametric noise are often used as paradigms for the study of these various eﬀects and related mathematical methods [10].The advantage of such models is that they have an appealing physical interpretation and appear as building blocks in many diﬀerent ﬁelds;they can be simulated on a computer or constructed as real electronic or mechanical systems [11,12].Moreover,the mathematical apparatus needed to analyze them remains relatively elementary (as compared to the perturbative ﬁeld-theoretical methods required for spatio-temporal systems [13])and can be expected to yield exact and rigorous results.In the present work,we study the motion of a frictionless pendulum with parametric noise,which can be physically interpreted as a randomly vibrating suspension axis.We show that the long time behavior of a stochastic pendulum driven by a colored noise with ﬁnite correlation time is drastically diﬀerent from that of a pendulum subject to white noise.Whereas the average energy of the white-noise pendulum is a linear function of time,that of the colored-noise pendulum grows only as the the square-root of time.Our analysis is based on a generalization of the averaging technique that we have used previouly for nonlinear oscillators subject to white noise:in [14,15,16]an eﬀective dynamics for the action variable of the system is derived after integrating out the fast angular variable,and is then exactly solved.However,this averaging technique as such ceases to apply to systems with colored noise when the time scale of the fast variable becomes smaller than the correlation time of the noise.Correlations between the fast variable and the noise modify the long time scaling behavior of the system and therefore must be taken into account.We shall develop here a method that systematically retains these correlation terms before the fast variable is averaged out.This will allow us to derive analytical expressions for the asymptotic probability distribution function (P.D.F.)of the energy of the stochastic pendulum,and to deduce the long time behavior of the system.Our analytical results are veriﬁed by numerical simulations.Finally,we shall compare our method and results with some known approximation schemes used for multivariate systems with colored noise.We shall show in particular thatsmallcorrelation time expansions cannot explain the anomalous diﬀusion exponent when truncated at anyﬁnite order.The partial summation of Fokker-Planck type terms,used to derive a‘best eﬀective Fokker-Planck equation’for colored noise[17,18]leads to results that agree with ours.We emphasize that the noise considered in this work has aﬁnite correlation time and its auto-correlation function does not have long time tails.This article is organized as follows.In section II,we analyse the case where the parametricﬂuctuations of the pendulum are modeled by Gaussian white noise,and explain heuristically why colored noise leads to anomalous scaling of the energy.In section III,we study the case of Ornstein-Uhlenbeck noise and explain how a recursive adiabatic elimination of the fast variable can be performed.This allows us to derive analytical results for the P.D.F. of the energy,which we validate with direct numerical simulations.In section IV,we compare our results with eﬀective Fokker-Planck approaches.Concluding remarks are presented in section V.II.THE PENDULUM WITH PARAMETRIC NOISEThe dynamics of a non-dissipative classical pendulum with parametric noise can be described by the following system of stochastic diﬀerential equations˙Ω=−(ω2+ξ(t))sinθ,(1)˙θ=Ω,(2) whereθrepresents the angular displacement andΩthe angular velocity.The energy E of the system isE=Ω2∂t =−∂∂Ω ω2sinθP t +D∂Ω2.(5)From Eq.(2),we observe that the angular variableθvaries rapidly as compared toΩ.Thus,following[14],we assume that,in the long time limit,the angleθis uniformly distributed over[0,2π].This allows us to average Eq.(5)over the angular variable and to derive an eﬀective Fokker-Planck equation for the marginal distribution˜P t(Ω):∂˜P t4∂2˜P twhere we have replaced sin2θby its mean value1/2.From this eﬀective Fokker-Planck equation,we readily deduce thatΩis a Gaussian variable with P.D.F.˜P t (Ω)=1πD texp −Ω222exp −2E4t,(9)and also the skewness andﬂatness factorsS(E)= E3 2)Γ(1Γ(5√E2 2=Γ(92)2)2=35whereαis an unknown exponent to be determined.Weﬁnd from Eq.(2)thatθ∼tα+1,and Eq.(1)can then be written as˙Ω≃ξ(t)sin tα+1.(13) (We could have retained the deterministic term−ω2sinθbut it would not aﬀect this scaling analysis.)We now take ξ(t)to be a discrete dichotomous noise with correlation timeτand with values±1.The previous equation,then, becomesΩ(t)∼t/τk=1ǫk kτ(k−1)τsin xα+1d x,withǫk=±1.(14)We estimate the last integral by an integration by parts:(α+1) t+τt xαsin xα+1xα t+τt+O(t−α−1),(15) and obtainΩ2 ∼t/τk=1 kτ(k−1)τsin xα+1d x 2∼t/τ k=112τe−|t−t′|/τ,(17)whereτis the correlation time of the noise.This noiseξcan be generated from white noise via the Ornstein-Uhlenbeck equation˙ξ=−1τη(t),(18) whereη(t)is a white noise of auto-correlation function Dδ(t−t′).In the stationary limit,t,t′≫τ,the solution of Eq.(18)satisﬁes Eq.(17).The pendulum with Ornstein-Uhlenbeck noise is thus written as a three-dimensional stochastic dynamical system coupled to a white noiseη(t):˙Ω=−ω2sinθ−ξsinθ,(19)˙θ=Ω,(20)˙ξ=−1τη(t).(21) The Fokker-Planck equation for the three-dimensional P.D.F.P t(θ,Ω,ξ)is given by∂P t∂θ(ΩP t)+∂τ∂2τ2∂2P tA.Zeroth-order averagingWe show here that the averaging procedureused in section II A forwhite noiseleads toerroneous results for colored noise.Averaging the Fokker-Planck equation (22)over the fast angular variable θ,we ﬁnd that the marginal P.D.F.˜Pt (Ω,ξ)obeys the Ornstein-Uhlenbeck diﬀusion equation ∂˜Pt τ∂2τ2∂2˜Pt τξ+12to yield Eq.(6).In the next subsections,we develop an averaging scheme that allows us to eliminate adiabatically the fast variable while retaining the correlation terms.The idea is to deﬁne recursively a new set of dynamical variables that embodies the correlations order by order.This will enable us to derive sound asymptotic results for the pendulum with colored noise.In this scheme,Eq.(24)appears as a zeroth order approximation,and its correct interpretation is not that Ωis conserved but that its variations are slower than that of normal diﬀusion.B.First-order averagingMultiplying both sides of Eq.(19)by Ω,and using Eq.(20),we obtainΩ˙Ω=−Ω ω2sin θ+ξsin θ=−ω2˙θsin θ−ξ(˙θsin θ).(25)Introducing the energy E of the system deﬁned in Eq.(3),this equation becomesdE dtΩ2dt=ddt(E −ξcos θ)=cos θτη.(28)This leads us to deﬁne a new dynamical variable E 1E 1=E −ξcos θ=Ω2τξ−cos θ2 E 1+(ω2+ξ)cos θ ,(31)˙ξ=−1τη(t ).(32)FIG.2:Stochastic pendulum for ω=1.0and (D ,τ)=(1,0.1),(1,1),(1,10),(10over 104realizations.We plot the ratio E 2 /(D t/τ2)vs.time t .The dotted line in the ﬁgure corresponds to the asymptotic value 0.25predicted by second-order averaging (Eq.(55)).Note that convergence is slower for smaller values of the correlation time τ.The advantage of this system as compared to the previous one (19,20,21)is that the white noise η(t )now appears in the equation for the dynamical variable E 1and this white noise contribution will survive the averaging process.The Fokker-Planck equation for the P.D.F.P t (E 1,θ,ξ)associated with Eqs.(30,31and 32)reads:∂P t ∂E 1cos θ∂θ(Ω(E 1,θ,ξ)P t )+1∂ξ(ξP t )+D ∂E 21−2cos θ∂2P t ∂ξ2 .(33)Averaging this equation with respect to θleads to the following evolution equation for the marginal distribution˜Pt (E 1,ξ):∂˜Pt τ∂4τ2∂2˜Pt 2τ2∂2˜Pt √τξ+1√D twith E ≥0.(37)From this expression we calculate the ﬁrst two moments of the energyE =1π D t 2≃0.564D t2,(38) E 2 =1τ2.(39)These expressions provide scaling relations between the averages,the time t and the dimensional parameters of the problem,D and τ.We have veriﬁed numerically that these scalings are correct (see Fig.2).However,the prefactorsthat appear in Eqs.(38)and(39)are pure numbers and do not agree with the results of our numerical simulations. We conclude that Eq.(38)is exact at leading order as it gives the correct asymptotic scaling for the energy,E∝t1/2 but fails to provide the prefactors.The reason is that some correlations betweenθandξhave still been neglected in the averaging procedure.Carrying out the calculations to the next higher order will enable us to derive the correct expressions for the prefactors.C.Second-order averagingMore precise results can indeed be derived by applying recursively the procedure described above.In the Langevin equation(30)for E1,we perform one more‘integration by parts’and obtain˙E 1=ξcosθΩ−cosθdt ξsinθτΩ˙ξ+˙Ωτξ−cosθτΩ=Ω2τΩ.(41)Using Eqs.(19)and(21),Eq.(40)becomes˙E 2=−ξsinθτΩ2(ω2+ξ)− cosθ+sinθτ.(42)In this equation we must express the variableΩin terms of E2,θandξ.Inverting the relation(41),we deduce thatΩ=(2E2)1(2E2)1τ(2E2)+O 122 ,(43) where we have retained terms up to the order1/E2.From Eq.(42),we deduce the Langevin equation for E2˙E 2=J E(E2,θ,ξ)+D E(E2,θ,ξ)η(t)E3τ2(2E2)1τ(2E2)(ω2+ξ),(45)D E(E2,θ,ξ)=− cosθ+sinθ2 .(46)This equation,combined with Eqs.(20and21),deﬁnes a three-dimensional stochastic system for the variables (E2,θ,ξ).The Fokker-Planck equation for the P.D.F.P t(E2,θ,ξ)is∂P t∂E2(J E P t)−∂τ∂2τ2 ∂∂E2(D E P t)+∂2∂E2D E∂P t∂ξ2 .(47)We now integrate out the fast angular variableθfrom Eq.(47),retaining only the leading term in the average of the expression∂∂E2(D E P t),(recalling that∂/∂E2scales as E−12,the contribution of the subdominant terms is ofthe order of E−5/22and is negligible in the long time limit).We thus obtain the following evolution equation for themarginal distribution˜P t(E2,ξ):∂˜P t∂E2 ω2ξ+ξ2τ∂2τ21∂E22+∂2˜P tAlthough the cross-derivative terms between E2andξvanish,the two variables are coupled through the drift term. From Eq.(48)we derive an eﬀective,two-dimensional,stochastic system in E2andξ:˙E 2=−ω2ξ+ξ2√τξ+18τ2E2+12τη2(t)−18τ2E2+12τη2(t).(52)The variables E2andξare decoupled at large times:the eﬀective problem is thus one-dimensional.In the long time limit,the variable E2is identical to the energy E up toﬁnite terms.We thus obtain the asymptotic P.D.F.of the pendulum’s energy by explicitely solving the Fokker-Planck associated with Eq.(52)˜P t (E)=2√Γ 12exp−τ2E22π4)2 D tτ21/2,(54)E2 =1τ2.(55) Besides the skewness and theﬂatness factors areS(E)=Γ(74)4)3/2=6√Γ(1 4)Γ(1Γ(5FIG.3:τ=1.(a)the average E and the ratio E /(D t/τ2)1/2(inset),(b)the skewness andﬂatness factors of E vs.time t.The asymptotic behavior of all observables presented agrees with Eqs.(54–57)(dotted lines in theﬁgures),irrespective of the value ofω.the predictions of Eqs.(54-57).In particular,we notice that the asymptotic P.D.F.˜P t(E),and therefore the moments of the energy,do not depend on the value of the mean frequencyω(see Fig.3).Our averaging technique thus provides sound asymptotic results for the energy of the stochastic pendulum:this technique not only yields the correct scalings but also leads to analytical formulae for the large time behavior of the energy.This averaging method could be carried over to the third order to calculate the subdominant corrections to the P.D.F.of the energy.However,we shall not pursue this course any further:the calculations become very unwieldy and the agreement between the analytical results and the numerical computations is already very satisfactory at the second order.We emphasize that the white noise and the colored noise cases fall in two distinct universality classes because the long time scaling exponents are diﬀerent.The colored noise scaling will always be observed after a suﬃciently long time provided that the correlation timeτis non-zero.However,the eﬀect of this correlation time appears only when the period T of the pendulum is less thanτ.This period,which is proportional toΩ−1,decreases with time.At short times,the period T is much greater thanτand white noise scalings are observed.At large times,T≪τand colored noise scalings are satisﬁed.The crossover time t c is reached when T∼Ω−1∼τ;using Eq.(9),which is valid for t<t c,we obtaint c∼(Dτ2)−1.(58) Hence,when the correlation timeτbecomes vanishingly small,the crossover time diverges to inﬁnity and the colored noise regime is not reached(simulation times needed to observe the colored noise scalings become increasingly long). In Fig.4,we plot the behavior of the mean energy vs time for D=1.0andτ=0.01,0.1and1.Forτ=1,the colored noise scaling regime is obtained from the very beginning and the curve for E has a slope1/2in log-log scale.For τ=0.1,at short times E ∝t whereas at long times E ∝t1/2.The crossover is observed around t c∼100=τ−2. Forτ=0.01,the curve for E has a slope1and the colored noise regime is not reached in this simulation(t c∼104). The averaging method provides analytical results for t≫t c and t≪t c.The intermediate time regime is described by a crossover function of the scaled variable t/t c[15],which cannot be analyzed by our technique.In the next section,we compare our method with two well-known approximations based on eﬀective colored noise Fokker-Planck equations.One of these approaches will enable us to derive approximate formulae for t c and for the crossover function.PARISON WITH OTHER APPROXIMATION SCHEMESThe main diﬃculty for the study of a Langevin equation with colored noise stems from its non-Markovian character [1]:there exists no closed Fokker-Planck equation that describes the evolution of the P.D.F.of the dynamical variables. The process can be embedded into a Markov system if the colored noise itself is treated as a variable.However,this mathematical trick increases the dimension of the problem by one and the noise has to be integrated out in the end: exact calculations can be carried out with this method only in the case of linear problems.FIG.4:Stochastic pendulum forω=0.0, D=1.0andτ=0.01,0.1and energy vs. time.The crossover between white noise behavior( E ∼t)and colored noise behavior( E ∼t1/2)occurs later for smaller correlation timeτ.A colored noise master equation can be rigorously derived(e.g.,with the help of functional methods)but it involves correlations between the dynamical variables and the noise[23,24].The equation of motion for these correlations involves higher order correlations,and so on.Since this hierarchy must be stopped at some stage,this question is a genuine closure problem.Many diﬀerent approaches have been devoted to derive eﬀective Fokker-Planck equations, such as short correlation time expansions[17,18,25,26],uniﬁed colored noise approximation[27],projection methods [28]and self-consistent decoupling Ans¨a tze[29](for a general review see[30]).In section III,to derive the long time behavior of the stochastic pendulum with colored noise,we did not use any closure approximation but started from the exact Fokker-Planck equation for the system,treating the noise as an auxiliary variable.Thus,we did not make any hypothesis onτand our analytical formulae are valid for any value of the correlation time(for t larger than the crossover time given in Eq.(58)).In this section,we compare our results with those that can be derived from some well-known approximation schemes.A.Small correlation time expansionAn approximate evolution equation for the P.D.F.of a Langevin equation with colored noise can be derived in the case of short correlation times(i.e.,in the white noise limit)by expanding the colored noise master equation around the Markovian point[26,30].This procedure leads to a Fokker-Planck type equation,with eﬀective drift and diﬀusion coeﬃcients.Applying to our system(19,20,21)the smallτexpansion derived in[26]for arbitrary stochastic equations with colored noise,we obtain,atﬁrst order inτ,the eﬀective Fokker-Planck equation for P t(θ,Ω)∂P t∂θ(ΩP t)−3Dτ∂Ω+D(1−τΩ)∂Ω2+Dτ∂θ∂Ω sin2θP t .(59)To simplify the discussion,we have taken the mean frequencyωequal to zero.After integrating out the rapid variations ofθan averaged Fokker-Planck equation is obtained which is similar to Eq.(6),but with an eﬀective diﬀusion given by D eﬀ=D(1−τΩ)(because D eﬀbecomes negative for large values ofΩ,Eq.(59)is valid only over the restricted region of positive diﬀusion).From dimensional analysis,we observe that such an eﬀective Fokker-Planck equation leads to a normal diﬀusive behavior,Ω∼t1/2and E∼t,and therefore cannot account for the results we have obtained.B.Best Eﬀective Fokker-Planck equationWe now turn to another smallτapproximation,which intends to improve theﬁrst-order eﬀective Fokker-Planck equation(59)by summing contributions of the type Dτn(where n is an integer larger than1).The resulting equation has been christened‘Best Fokker-Planck Equation’(B.F.P.E.)by its proponents[17,18,31].Although this approach is not free from drawbacks and is known to lead in some cases to unphysical results[32],we show that for the system studied in this work,the B.F.P.E.leads to results that agree with ours.An approximate evolution equation for the P.D.F.P t(Ω,θ)is given by the second-order cumulant expansion[1]of the(stochastic)Liouville equation associated with Eqs.(19and20):∂P t∂θ(ΩP t),(61)L1(t)P t=−∂∂t =−∂2∂21+(τΩ)2P t +Dτ∂Ωsinθ∂(1+(τΩ)2)2P t .(63)Integrating out the fast variableθ,we obtain an averaged B.F.P.E for˜P t(Ω),the probability distribution of the slow variable,∂˜P t4∂1+(τΩ)2∂˜P t4 Ω4 +14t.(65)When t is small,the quadratic term dominates over the quartic term,and we recover Ω2 =2 E ≃D4 Ω4 ≃1τ2t.(66)This result is identical to our Eq.(55),which was validated by numerical results(see Fig.2).The identity(65)can also be used to derive an approximate scaling function for the mean energy.Let us deﬁne theﬂatnessφof˜P t(Ω)asΩ4 =φ Ω2 2.(67) Rigorously speakingφis a function of time,but it remains a number of order1.For simplicity,let us assume thatφis constant.Substituting Eq.(67)in Eq.(65),and solving for Ω2 ,we obtainE =1Dτ2φt+1−1V.CONCLUSIONA nonlinear pendulum subject to parametric noise undergoes a noise-induced diﬀusion in phase space.The charac-teristics of this motion depend on the nature of the randomness:when the noise is white the energy of the pendulum grows linearly with time,whereas it varies only as the square root of time when the noise is colored.This change of behavior is due to destructive interference between the displacement of the pendulum and the noise term:the eﬀect of the colored noise is partially averaged out by fast angular variations as soon as the period of the pendulum becomes smaller than the correlation time τof the noise.We have carried out an analytical study of this model by deﬁning recursively new coordinates in the phase space and averaging out the fast angular variable.At zeroth order,the only information obtained is that the motion is subdiﬀusive;at ﬁrst order,this procedure provides the correct scalings;at second order,a quantitative agreement with numerical simulations is reached.We emphasize that our method is diﬀerent from the usual approximations that involve eﬀective Fokker-Planck equations in which the colored noise does not appear as an auxiliary variable.Our averaging procedure integrates out the fast dynamical variable and leads to an eﬀective stochastic dynamics for the slow variables with colored noise.Whereas usual eﬀective Fokker-Planck equations are valid only for small noise and for short correlation times,we do not make any hypothesis on the amplitude of the noise or on the correlation time τ.However,the asymptotic subdiﬀusive regime is reached earlier for larger values of τ.Our results agree with those derived after averaging the ‘Best Fokker-Planck Equation’(B.F.P.E.)approximation (which is obtained from a summation of a truncated cumulant expansion of the Liouville equation).This approximation is also useful to draw a qualitative physical picture of the system and allows to calculate approximate crossover functions between short time (white noise)and long time (colored noise)regimes.The recursive averaging scheme that we have used here for the stochastic pendulum can be extended to other systems coupled with colored noise.In particular,we believe that thanks to this method a precise mathematical analysis of the long time behavior of a nonlinear oscillator subject to multiplicative or additive colored noise can be carried out [33].AcknowledgmentsWe thank S.Mallick for his helpful comments on the manuscript and F.Moulay for useful discussions.APPENDIX A:DERIV ATION OF EQ.(63)In this appendix,we derive the B.F.P.E.(Eq.(63))following the procedure of [17,18].We evaluate the right hand side of Eq.(60)by applying the following operator formulaexp(A )B exp(−A )=B +[A,B ]+13![A,[A,[A,B ]]]+...,(A1)with A =L 0and B =L 1(t −x ).We thus have to calculate some commutators of the two operators L 0and L 1(t −x )deﬁned in Eqs.(61)and (62).By induction,we derive the following expression for the n -th commutatorT n =[L 0,[...,[L 0,[L 0,L 1(t −x )]]...]=ξ(t −x ) ∂∂θH (n )2(Ω,θ) ,(A2)where the functions H (n )1and H (n )2satisfy the recursion relationsH (n )1=−Ω∂H (n −1)1∂θ.(A4)The ﬁrst few terms can be calculated explicitely and we obtainH (0)1=−sin θand H (0)2=0,(A5)H (1)1=Ωcos θand H (1)2=−sin θ,(A6)H (2)1=Ω2sin θandH (2)2=2Ωcos θ.(A7)The general solution for the recursion(A4)is readily found:H(n)1=(−1)n−1Ωn sin(θ+nπ2).(A9)From Eqs.(A1)and(A2),we deduce the following identityL1(t)exp(L0x)L1(t−x)exp(−L0x) = ∞n=0x n∂Ωsinθ ∞n=0x n∂ΩH(n)1(Ω,θ)+∂2τ∂∂ΩH1(Ω,θ,t)+∂n!H(n)1=∞n=0 t/τ0d x x n e−x2),(A12)H2(Ω,θ,t)=∞n=0 t0d xx n e−x/τn!(−1)n−1τn+1nΩn−1cos(θ+nπ1+(τΩ)2,(A14) H2(Ω,θ,∞)=−τ2(1−(τΩ)2)sinθ−2τΩcosθ[18]Lindenberg K and West B J1984Physica A12825[19]Bourret R C,Frisch U and Pouquet A1973Physica65303[20]Van Kampen N G1976Phys.Rep.24171[21]Lindenberg K,Seshadri V and West B J1981Physica A105445[22]Lindenberg K,Seshadri V and West B J1980Phys.Rev.A222171[23]H¨a nggi P1985in Stochastic Processes Applied to Physics edited by Pesquera L and Rodriguez M A(Singapore:WorldScientiﬁc)[24]H¨a nggi P1989in Noise in Dynamical Systems Vol.1edited by Moss F and Mc Clintock P V E(Cambridge:CambridgeUniversity Press)[25]San Miguel M and Sancho J M1980Phys.Lett.A7697[26]Ramirez-Piscina L and Sancho J M1988Phys.Rev.A374469[27]H’walisz L,Jung P,H¨a nggi P,Talkner P and Schimansky-Geier L1989Z.Phys.B77471[28]Fox R F1986Phys.Rev.A33467[29]Fronzoni L,Grigolini P,H¨a nggi P,Moss F,Mannella R and Mc Clintock P V E1986Phys.Rev.A333320[30]H¨a nggi P and Jung P1995Adv.Chem.Phys.89239[31]Peacock-Lopez E,de la Rubia F J,Lindenberg K and West B J1989Phys.Lett.A13696[32]H¨a nggi P,Marchesoni F and Grigolini P1984Z.Phys.B56333[33]Mallick K and Marcq P in preparation。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

A comparison of discrete element simulations and experiments forF sandpiles _composed of spherical particlesYanjie Li a ,Yong Xu a,*,Colin Thornton ba Department of Applied Mechanics,China Agricultural University,Beijing 100083,China bSchool of Engineering,University of Birmingham,Edgbaston,Birmingham B152TT,UKReceived 11April 2005;received in revised form 1August 2005;accepted 6September 2005Available online 18October 2005AbstractDiscrete element simulations,with the particle–particle interaction model based on classical contact mechanics theory between two non-adhesive spheres,were carried out and compared with F sandpile _experiments using spherical particles in order to assess the validation of the simulation data.The contact interaction model is a combination of Hertzian theory for the normal interaction and Mindlin–Deresiewicz theory for the tangential interaction.To ensure the consistency of the simulations with the experiments,the measurement of sliding friction,a key parameter in DEM simulations,was highlighted.A simple experimental method for establishing the value of the friction coefficient was proposed and used in measuring the friction of the rough glass beads and steel balls to be modelled in the simulations.The simulations were carried out for two cases according to particle arrangements:the first is quasi-two-dimensional (Q2D),with a smaller flat cuboidal box containing the spherical particles inside another box for discharge,and the second case is axisymmetric (3D).For both cases,simulations and experiments were carried out for assemblies of polydisperse rough glass beads under the same parisons were made that showed that the profiles and hence the measured angles of repose,in each case,were in good agreement,thus supporting the validity of the discrete element model used.Further numerical–experimental comparisons were carried out for 3D conical piles using smooth monodisperse steel balls and the same conclusions were obtained.D 2005Elsevier B.V .All rights reserved.Keywords:Discrete element method;Particles;Sandpiles;Angle of repose;Base pressure distribution1.IntroductionThe Discrete Element Method has become a powerful numerical method for analysing discontinuous media since the important work by Cundall and Strack [1].Apart from the advances of applications in civil engineering with block element modelling,the particulate discrete element method for granular materials has developed rapidly in analysing both macroscopic and microscopic behaviours in many applica-tions in chemical engineering,agricultural engineering,etc.Progress has involved both model improvement and intensive simulations.Many investigators have proposed models modified from the original for specific applications.Thornton and Yin [2],Thornton [3,4]developed new particle interac-tion models for spherical particles based on contact mechanics for both nonadhesive and adhesive particles as well as for elastoplastic interactions.Oda et al.[5]modified the conventional linear spring–dashpot model by adding rolling resistance with which he analysed shear band phenomena.Han et al.[6]combined the discrete element method with the finite element method to simulate peen-forming processes.Williams et al.[7]developed a contact detection algorithm for particles with arbitrary geometries.Cleary [8]reported DEM simulations illustrating applications to various industrial processes.However,many researchers have done DEM simulations with their own models without experimental validation.Only a small number of studies have presented their predicted results in comparison with experimental data with enough specified parameters.The ambiguous parameter specification leads to theoretical uncertainty and thus restricts its applicability.For instance,the model based on classical contact mechanics between0032-5910/$-see front matter D 2005Elsevier B.V .All rights reserved.doi:10.1016/j.powtec.2005.09.002*Corresponding author.Tel.:+861062736514;fax:+861062736514.E-mail address:xuyong@ (Y .Xu).Powder Technology 160(2005)219–228/locate/powtectwo spheres is theoretically rational,however,up to date there is little experimental evidence to provide complete validation.A typical problem related to many natural phenomena and various engineering applications are F sandpiles_.The superfi-cial shape of a heap is actually a macroscopic result formed from various effects of the internal mechanisms involving many physical factors.Therefore this approach is helpful for understanding the microscopic mechanics and the interactive mechanics.On the other hand,because laboratory work on sandpiles is relatively easy,it is a useful benchmark problem for evaluating the applicability of a specific model used in DEM simulations.Zhou et al.[9,10]studied sandpiles of glass beads using both experiments and simulations.They used a modified DEM model in which the effect of rolling resistance was considered.In their consequent study both the rolling and sliding friction effects were examined,together with the influence of density and size of the particles. Comparisons between the numerical and simulated results were quite good.Matuttis et al.[11]used the Discrete Element Method to simulate dense packing and heaps and analysed the stress distribution with both spherical and non-spherical particles.In the simulations reported below,sliding friction is considered but rolling friction is ignored.To ensure that the simulations and corresponding experiments have the same condition,an easy-to-use experimental methodology for measuring sliding friction is adopted to establish the value of the friction coefficient to be specified for the simulations.2.Discrete element model for dry spheres without adhesionThe discrete model used in this approach is for dry spherical particles without adhesion,the interactions include normal and tangential inter-particle forces.The interactions between spheres1and2,with radius R1and R2and center coordinates (x1,y1,z1)and(x2,y2,z2)respectively,are shown in Fig.1.The normal force is determined by Hertzian theory asN¼4E4R*ðÞ1=2a3=2ð1Þwhere a,R*and E*are defined as a¼R1þR2Àﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃx2Àx1ðÞ2þy2Ày1ðÞ2þz2Àz1ðÞ2q;a>0ðÞð2Þ1 R4¼1R1þ1R2;1E4¼1Àm21E1þ1Àm22E2ð3ÞIn an incremental scheme with a time step,D t,anincremental relative approach of the two spheres is D a,then the incremental normal contact force is calculated byD N¼2E4a D að4Þwhere a¼ﬃﬃﬃﬃﬃﬃﬃﬃa R*pis the radius of the contact area.The tangential force is modelled by the theory of Mindlin and Deresiewicz[12].When two contacting surfaces are subjected to an increasing tangential displacement,d,then relative slip is initiated at the perimeter and progresses inward over an annular area of the contact surface.The incremental tangential force D T due to the incremental tangential displace-ment D d depends not only on the loading history but also on the variation of the normal force.Therefore D T may be obtained from the following equation[20]D T¼8G4a h k D dþÀ1ðÞk l D N1Àh kðÞð5Þin which h k depends on the loading status.That is,if |D T|<l D N,then there is no slip so h k=1,otherwise the slip effects need to be considered ash k¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ1ÀTþl D Nl N3rk¼0loadingðÞﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ1ÀÀ1ðÞk TÀT kðÞþ2l D N2l3sk¼1;2unloading and reloadingðÞ8>>><>>>:ð6Þwhere l is the sliding friction coefficient,and G*is defined as,1G4¼2Àv1G1þ2Àv2G2:ð7ÞIn Eq.(6),T k represents the historical tangential force from which unloading or reloading commenced and it needs to be updated as T k=T kÀ(À1)k l D N to allow for the effect of the variation of the normal force.For dynamic processes it is necessary to consider the elastic wave propagation across particles,the time for load transfer from one particle to adjacent contacting particles and the need not to transmit energy across a system faster than in nature.For linear contact stiffnesses the critical time step is related to the ratio of the contact spring stiffness to particle density.For non-linear springs(e.g.Hertz)the critical time-step cannot be calculated a priori.However,it was shown by Miller andFig.1.The forces between two spheres without adhesion.Y.Li et al./Powder Technology160(2005)219–228 220Pursey (see Ref.[13])that the Rayleigh waves account for 67%of the radiated energy in comparison with the dilational (7%)and distortional (26%)waves.In the simulations,it is assumed that all of the energy is transferred by the Rayleigh waves.This is a good approximation since the difference between the Rayleigh wave speed and the distortional wave speed is very small and the energy transferred by the dilational wave is negligible.In addition,the average time of arrival of the Rayleigh wave at any contact is the same irrespective the location of the contact point.For simplicity,the critical time-step is based on the average particle size and a fraction of this is used in the simulations.Therefore,the critical time-step is given by the following equation.D t c ¼k R¯b ﬃﬃﬃﬃq G r ð8Þwhere R¯is the average particle radius,q is the particle density,G is the particle shear modulus and b can be obtained [13]from2Àb 2ÀÁ4¼161Àb 2ÀÁ1Àb 21À2m 21Àm ðÞ !ð9Þwhich can be approximated [20]by b ¼0:8766þ0:163mð10Þwhere m is the Poisson ratio of the particle.For the simulations presented below the critical time-steps were 3.98E À4s (glass)and 2.14E À4s (steel)and the actual time-steps used are listed in the relevant tables.3.Experimental method for measuring sliding friction In order to compare the simulated results with the corresponding experiments,all parameters involved in the simulation have to be the same as those in the experiments,such as particle number with size,friction between particle–particle and particle–wall and other physical properties.Among them friction is a key factor which needs to be established by means of experiments.According to modern theories of friction,sliding friction depends not only on the surface roughness,but also on the sliding velocity and contacting area and is related to molecular interactions.It is rather complex and there is no common or convenient formula available.Therefore the classical friction theory is applied for simplicity in this approach,which states that l ¼F =Nð11Þin which F is the frictional force and N is the normal force,the effect of contact area is ignored,therefore it is suitable for point–point contact.According to the classical friction theory,the frictional effect between two contacting particles with the same material and surface condition can be considered as point–point contact for the real contact area is very small.This is also the case between a particle and a plane-wall.Therefore the sliding friction between two particles can be measured for a particle and a wall (with the same material and surface condition)instead of a pair of spheres.Based on this concept a methodology was proposed and a simple device for measuring sliding friction of two particle surfaces was set up.Fig.2(a)and (b)show a sketch and the actual device,which was used for measuring the friction of rough glass beads as reported below.The ‘‘particle-board’’(the upper one)was designed as a group of particles of the same size that were glued firmly to a small plate through a foam-rubber pad so that all the lowest points of the spheres make contact with the base plate (the lower one)on a common plane to form multi-point contacts.The base plate is made of the same material as the particles and the upper surface has the same roughness.When theparticle-board,(b)Fig.2.Experimental setup for measuring sliding friction.(a)Sketch;(b)actual device.50100150200250F r i c t i o n a l f o r c e F (g )Pressure force N (g)Fig.3.Relationship of friction force N and normal force F .Table 1Experimental values of sliding frictionGlass ball –glass plate Glass ball –perspex plate Steel ball –steel plate 0.15450.13330.2140Y.Li et al./Powder Technology 160(2005)219–228221with the additional vertical load,is driven by a pull force via a dragline to slide at a low constant velocity along the base plate,the frictional force F(equal to the pull force)and the vertical normal force N can be measured and the sliding friction can be calculated using Eq.(11).This design allows larger contact forces and hence the measurement can be more accurate and stable.For a DEM simulation,the friction between a particle and the wall of a container and between two particles with different material properties need to be specified.The same methodology as above can be applied using a particle-board and another base plate with specific properties as the user ing the device shown in Fig.2b,three groups of experiments were carried out: (1)Rough glass ball–rough glass plate for inter-particlefriction;(2)Smooth steel ball–smooth steel plate for inter-particlefriction;(3)Rough glass ball–Perspex plate for friction between theglass ball and the container.Typical results of the measured frictional forces F and the normal forces N for the three cases are shown in Fig.3,from which it is seen that the frictional force is approximately proportional to the normal force for all three cases,suggesting that the effect of the contact area on sliding friction is not significant and therefore the classical friction theory is applicable.The measured values of sliding friction to be used in the simulations were established by mathematically averag-ing the results and are given in Table1.4.Quasi-two-dimensional(Q2D)sandpiles of glass balls 4.1.Experimental setupThe device for a Q2D experiment consists of two flat boxes of Perspex glass with different dimensions.The smaller is the inner box with two optional outlets for discharge via the base. One outlet is at the side and the other is central.The larger outer box allows the inner box to move vertically and collects the discharged particles as they form a heap.Fig.4(a)and(b) provide the dimensions of the two boxes used.The depth of the inner box is approximately5.5times the average particle size. This is reasonable to reduce the frictional influence of the front and the back walls and to avoid excessive2D kinematic constraints on the individual particle movements.In the Q2D experiments,3000rough glass beads were used with a mean diameter of11.175mm,the particle size distribution is given in Table2.The physical parameters used in the simulations are listed in Table3.The values of the physical parameters are either found in relevant standard handbooks or through experiment.The particles of different sizes were mixed as a homoge-neous poly-dispersed system.All the mixed particles were poured into the inner box,which had been positioned inside the outer box.Then the inner box was vertically raised gently to allow the glass beads to flow slowly into the outer box and settle down to form a wedge.The progressive profiles for Q2D side piling experiment are shown in Fig.6(a),(c)and(e),and the whole process of experiments were taken until both the particles in the outer box and those remaining in the inner box were stationary.The repose angles were then measured from the inclined planes of the wedges for both the side discharge and the central discharge experiments.Each experiment was repeated several times to eliminate errors from random influences.The final experimental value of the angle of repose was obtained by averaging the measured results obtained from several experiments.(a)(b)Fig.4.Box dimensions used for Q2D experiments.(a)Outer box;(b)inner box.Table2Size distribution of glass beadsGroup123Diameter(m)0.0130.0120.011Number7501500750Table3Parameters used in the simulationsVariable ValueParticle density(Mg/m3) 2.456Young’s modulus of particle(GPa)55Young’s modulus of container(GPa)200Particle–wall sliding friction0.1333Particle–particle sliding friction0.1545Poisson’s ratio of particle0.25Poisson’s ratio of container0.3Time step used(s) 3.98EÀ5Y.Li et al./Powder Technology160(2005)219–2282224.2.Simulation methodThe original state was prepared by random generation of all the different sized particles,from the largest to the smallest,in the inner box.A gravitational field was introduced so that the particles rained down,contacted each other and settled down to form a bed.Then the actual sandpile simulation was started.The inner box was raised at a very low speed (0.003¨0.005m/s)so that the particles could flow downwards through the outlet and discharge into the outer box to form a wedge.Fig.5(a)and (b)show the particle configuration and the particle velocity field respectively atan(a)(b)Fig.5.Central discharge simulation.(a)Particle configuration;(b)particle velocity field.(a) 0.0 s (b) 0.0s(c) 5.0 s (d) 5.0s(e) 15.0 s (f) 15.0 sFig.6.Progressive comparisons between experimental and simulated profiles for side discharge.(a),(c),(e)Are the instant profiles at 0.0,5.0and 15.0s during experiment;(b),(d),(f)are the corresponding simulated ones.Y.Li et al./Powder Technology 160(2005)219–228223instant time for central discharge.Fig.6(b),(d)and (f)show the simulated progressive profiles for Q2D side discharge in comparison with those from experiment shown in Fig.6(a),(c)and (e)respectively.It is evident that there is no significant difference between them.4.3.Results and discussionThe simulated profiles for central discharge and the experimental result are shown in Fig.7(a)and (b)respectively,while those for side discharge are shown in Fig.8(a)and (b)respectively.Comparing the simulated profiles for central discharge shown in Fig.7(a)with those obtained from experiments in Fig.7(b),it can be seen that the simulated results are in good agreement with the experimental ones despite the small difference in projective geometry (The graphic code is designed by point projection).Comparing both profiles for side discharge,Fig.8(a)and (b),the same conclusion can also be made,suggesting that the simulated results together with the DEM model are correct.Note that the comparison is for the profiles in the outer box because the profiles in the inner box are affected by the friction between the base wall and the particles whereas the profiles in the outer box only depend on the inherent interactions,the influence from the base wall is eliminated.However,when we examine the profiles for side discharge it is observed that neither the simulated nor the experimental profiles are paring these with those obtained for central discharge it is noted that many more particles remain inthe inner box than for central discharge.As a result,insufficient particles were discharged into the outer box in order to eliminate wall effects.Therefore,the measurement was taken over the middle segment to provide an average value.Based on the pile profiles,values of repose angle from both simulations and experiments are given in Table 4.From Table 4,it is seen that the two groups of results are quite close,thereby justifying the use of the DEM technique,with the contact interaction laws described above,for simulating this type of problem.5.Conical piles5.1.Experimental techniqueThe experimental device used to produce conical piles is shown in Fig.9(a),which consists of a round funnel for central discharge and different sized circular baffles,shown in Fig.9(b)and (c)respectively.A suitable round baffle was used for preventing particle rolling and thus ensuring that the particles build up to a typical conical pile.Preliminary experiments indicated that,without such a baffle the slope of the pile can not keep constant,it varies from a higher value gradually to a lower value.This is because it is significantly dependent on the frictional situation of the base plate,especially for spherical particles.Only by using the baffles can the friction effect on the repose angles be ignored.Two different types of spherical particles were used in these experiments.One type is the rough glass beads used in the Q2D experiments with the same physical parameters and size distribution.The other type are mono-sized steel beads usedin(a)(b)parison of simulation and experiment,central discharge (glass beads).(a)Simulated;(b)experimental.(a)(b)parison of simulation and experiment,side discharge (glass beads).(a)Simulated;(b)experimental.Y.Li et al./Powder Technology 160(2005)219–228224bicycle bearings,whose particle properties are listed in Table 5.Details of the baffle are given in Table 6.The experimental procedure is simple.In order to create an axisymmetric upright conical pile of particles,the funnel is fixed with a sleeve so as to move vertically along a metal pole fixed with a tripod,as shown in Fig.9(a).During discharge the funnel,with particles,is lifted up along the vertical symmet-rical axis formed by the intersection point of the tripod and the centre of the baffle.The diameter of the funnel outlet is chosen as more than 5–6times the maximum particle size so as to allow the particles to flow out avoiding arching.The funnel,with particles,is lifted slowly up from the lowest position until the discharged particles form a full pile within the round baffle.The experiments with the poly-dispersed rough glass beads and the mono-dispersed steel beads were repeated several times and typical results are shown in Figs.10(b)and 11(b)parison between simulation and experiment A DEM simulation of conical piling is not different from the Q2D one because the latter is also three-dimensional.The only difference is that specifications of the funnel and the round baffle are approximated by specifying a significant large number of plane walls.The comparison on the progressive profiles for conical piling between experiment and simulation is similar toQ2D therefore it is omitted in this paper for simplicity.The final profile simulated for rough glass beads is shown in Fig.10(a)and those for smooth steel beads is shown in Fig.11(a)with the corresponding experimental results shown in Figs.10(b)and 11(b)respectively.The final values of the repose angles for glass beads and steel beads are listed in Table 7.Comparing the simulated profile,Fig.10(a),and experi-mental one,Fig.10(b),for the rough glass beads,it is observed that the simulated profile is in good agreement with the experimental one.The values of repose angle are almost identical (23.2-and 22.7-).Similarly,the simulated profile for the smooth steel beads,Fig.11(a),and the corresponding experimental one,Fig.11(b),are also in good agreement (24.5-and 24.2-).This is further validation of the DEM code used,at least for the sandpile problem.When comparing the conical piling of the rough glass beads with the Q2D case,either simulation or experimental,it is clear that the values of the angle of repose are always smaller than those of Q2D piling which are 3–3.5-larger than those for conical piling.This can be explained in the following manner.The frictional effect due to the front and back walls can affectTable 4Angle of repose for Q2D piles of glass ballsExperimental resultsSimulation results Central discharge 26.3-26.9-Side discharge26.4-27.2-235mm55mm(a)(b)(c)Fig.9.Experimental device for conical piling.(a)Device set up;(b)funnel;(c)round baffle.Table 5Parameters used for steel balls VariableValueNumber of particles7000Average diameter of particles (mm) 6.33Particle density (Mg/m 3)7.8Young’s modulus of particles (GPa)200Poisson’s ratio of particles0.25Particle –particle sliding friction coefficient 0.2140Time step (s)2.14E À5Y.Li et al./Powder Technology 160(2005)219–228225the results obtained from the Q2D simulations and experi-ments,as demonstrated by Zhou et al.[10].However,the more significant difference between the 2D wedges and the 3D conical piles is the difference in kinematic constraint.It is well established that,although not the same,there is a strong correlation between the angle of repose and the angle of internal shearing resistance (angle of bulk friction)observed during quasi-static shear deformation of granular materials.Cornforth [14]was the first to demonstrate,for sand,that the maximum angle of internal shearing resistance is higher in plane strain than in axisymmetric compression.Recently,from DEM simulations of a polydisperse system of elastic spheres (E =70GPa,m =0.3,l =0.5)in a periodic cell Gong [15]has shown that the maximum angle of internal shearing resistance is 4–4.5-higher in plane strain than in axisymmetric compression.Furthermore,at large strains when the systems are deforming at constant volume,the plane strain angle of internal shearing resistance is ca.2-higher than that obtained in axisymmetric compression,irrespective of the initial packing density.It is,therefore,concluded that the difference in the repose angle obtained for the Q2D wedges and the conical piles is due to the different kinematic constraint/freedom existing in plane strain and axisymmetric flow conditions.5.3.Base pressure distributionFollowing the publication by Wittmer et al.[16]there followed a somewhat heated debate between physicists and engineers [17]concerning whether or not there should be adip in the distribution of normal pressure below the centre of a F sandpile _.However,most of the work reported in the literature,both experimental and numerical studies,have been restricted to two-dimensional situations.Results of 3D experiments (conical F sandpiles _)were reported by Huntley [18]but the data were inconclusive.A possible dip was identified for heaps composed of sand or glass beads but not for lead shot.It appears that the general consensus is that it depends on the depositional technique.Ignoring any effect due to deflection/settlement of the underlying base material then gravitational pluvial deposition from a point source,as in ore storage tipping,may result in a pressure dip but this does not occur when the F sandpile _is constructed in layers,as in highway embankments and dam construction.From the results of the simulations of conical F sandpiles _reported above,both the radial distribution of the base wall contact normal forces and the base normal pressure distribu-tions were examined and are shown in Figs.12and 13for steel and glass particles respectively.From Figs.12(a)and 13(a),it can be seen that the maximum force decreases with radial distance due to the conical shape of the pile and,for any radial distance,the magnitudes of the normal forces are randomly distributed between the maximum value and zero.This is due to the nature of force transmission in granular media,Thornton [19].For any annulus there will be contacts with the base that are the end-points of the strong force chains (i.e.in the strong force sub-network)and others that belong to the weak force sub-network.Figs.12(b)and 13(b)show the base pressure distributions for steel balls and glass balls respectively.In order to obtain the pressure,the circular base plate was divided into eight or tenTable 6Baffle details VariableMaterial type or value Container materialPerspex/plasticInner diameter u Âheight h (mm)395Â10(glass balls)/250Â10(steel balls)Young’s modulus (GPa) 2.0Poisson’s ratio0.3(a)(b)Fig.10.Profiles of conical piling of rough glass beads.(a)Simulated;(b)experimental.(a) (b)Fig.11.Profiles of conical piling of smooth steel beads.(a)Simulated;(b)experimental.Table 7Angle of repose for conical piling of glass balls and steel beadsSimulation Experiment Rough glass beads 22.7-23.2-Smooth steel beads24.5-24.2-Y.Li et al./Powder Technology 160(2005)219–228226。