CN113408124A - Particle system servo control method without changing boundary shape - Google Patents

Particle system servo control method without changing boundary shape Download PDF

Info

Publication number
CN113408124A
CN113408124A CN202110668184.9A CN202110668184A CN113408124A CN 113408124 A CN113408124 A CN 113408124A CN 202110668184 A CN202110668184 A CN 202110668184A CN 113408124 A CN113408124 A CN 113408124A
Authority
CN
China
Prior art keywords
particle
stress
particles
circle
particle system
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110668184.9A
Other languages
Chinese (zh)
Other versions
CN113408124B (en
Inventor
石崇
郝李坤
陈晓
张一平
司宪志
董家豪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202110668184.9A priority Critical patent/CN113408124B/en
Publication of CN113408124A publication Critical patent/CN113408124A/en
Application granted granted Critical
Publication of CN113408124B publication Critical patent/CN113408124B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a particle system servo control method without changing boundary shapes, and aims to solve the problem that the particle system servo in the prior art can change the external outline of a model. The method comprises the following steps: generating a plurality of particles in a geometric area according to a preset porosity to obtain an initial particle system; setting measurement circles in an initial particle system, and calculating the average stress inside each measurement circle and the average stress of the particle system; controlling the free movement of particles in the particle system by using particle flow software, and updating the size of the particles in each measuring circle according to the average stress in the measuring circle in the free movement process of the particles to obtain an updated particle system; and carrying out iterative calculation according to the error rate between the updated average stress of the particle system and the target servo stress to obtain an optimal particle system model. The method can ensure the accuracy of the state control of the particle flow numerical model under the condition of not changing the boundary of the model.

Description

Particle system servo control method without changing boundary shape
Technical Field
The invention relates to a particle system servo control method without changing boundary shape, belonging to the technical field of mechanical analysis.
Background
Numerical simulation is an important means for engineering mechanical analysis, wherein a particle flow method is an important means for researching a medium deformation failure mechanism. The particle flow is generally composed of boundary walls (two-dimensional lines and three-dimensional surfaces) and internal particles (two-dimensional discs and three-dimensional spheres), a certain contact model is arranged between the particles and the particle-wall, a time step is set to be small enough to allow the particle system to move freely, the contact force and the motion state of the particles are continuously updated at each time step, and when all the particles are in close contact and uniformly distributed at each position, the corresponding contact bonding model is given. The particle system obtained by bonding can simulate the mechanical properties of rock, concrete and other rock and soil materials.
Therefore, whether a particle flow numerical model can reasonably reflect the properties of the geotechnical medium depends on the state of the particle system during bonding, and the conventional method for controlling the state of the particle system is to move the wall of the particle system inwards or outwards according to a certain rule in a normal direction, so that the average stress of the particle system meets a certain design value, the process is called servo, the average stress is also called servo stress, and the particle system is considered to meet the requirement when the average stress on the boundary wall of the particle system is equal to the servo stress and the unbalance force of the system is close to 0. However, the existing servo control operation may change the outer contour of the model, resulting in distortion of the model range. If no servo is performed, the state of the particle system cannot be controlled, so that the assigned microscopic parameters cannot be matched with the macroscopic parameters. Therefore, when it is inconvenient to control the model profile strictly or to use boundary wall servos, the state of the particle system cannot be controlled precisely (each part is just close to the servostress), and the given parameters often cannot represent macroscopic mechanical properties.
Disclosure of Invention
Aiming at the problem that the servo of a particle system can change the external outline of a model in the prior art, the invention provides a particle system servo control method without changing the shape of a boundary, which not only can ensure that all parts of the particle system are close to servo stress, but also does not need to change the position of a boundary wall of the model, ensures the accuracy of state control of a particle flow numerical model, and can accurately reflect the mechanical properties of rock and soil materials.
In order to solve the technical problems, the invention adopts the following technical means:
the invention provides a particle system servo control method without changing boundary shape, which is characterized by comprising the following steps:
a, generating a plurality of particles in a preset geometric area according to a preset porosity to obtain an initial particle system;
b, setting regular measuring circles in an initial particle system according to the particle size of the generated particles, and calculating the average stress inside each measuring circle and the average stress of the particle system;
c, controlling free movement of particles in a particle system by using particle flow software, and updating the size of the particles in each measuring circle according to the average stress in the measuring circle in the free movement process of the particles;
d, judging and deleting edge particles according to the updated particle size to obtain an updated particle system;
and E, calculating the error rate between the updated average stress of the particle system and the target servo stress, and repeating the steps C to E according to the error rate to carry out iterative calculation to obtain an optimal particle system model.
Further, the step a comprises the following steps:
carrying out geometric boundary limitation according to the model building requirement to obtain a closed geometric region;
generating a plurality of particles according to a preset particle radius range and initial porosity, wherein the centers of the particles of all the particles are located in a closed geometric area;
judging the position of each particle by using a ray method, leading out a ray from the center of the particle of each particle to the right side, and acquiring the number of intersection points of the ray and the boundary of a geometric region, wherein the direction vector of the ray is (1, 0);
when the number of the intersection points corresponding to one particle is 0 or even, deleting the particle, otherwise, keeping the particle;
the retained particles are used to form the initial particle system.
Further, in step B, the method for setting a regular measurement circle in the initial particle system according to the particle size of the generated particles is as follows:
obtaining a measurement circle radius according to the maximum particle radius in the initial particle system;
obtaining a rectangular range corresponding to the geometric area according to the coordinates of the geometric area, and setting a plurality of tangent measuring circles in the rectangular range according to the radius of the measuring circles;
extending a ray from the center of each measuring circle to a vector (1,0) in the right direction, and acquiring the number of intersection points of the ray and the boundary of the geometric area;
and when the number of the intersection points corresponding to one measuring circle is 0 or even, deleting the measuring circle, otherwise, keeping the measuring circle.
Furthermore, the value range of the radius of the measuring circle is [5 r%max,10*rmax], wherein ,rmaxThe maximum particle radius.
Further, in step B, the calculation method of the average stress inside each measurement circle and the average stress of the particle system includes the following steps:
judging the particle-particle contact and the particle-wall contact in the particle system according to the distance between the centers of the particles, the distance between the centers of the particles and the boundary of the geometric area and the particle radius, and acquiring the total number of the contact in the particle system;
traversing all contacts in the particle system, calculating the distance d from the m-th contact to the k-th measurement circlemkWherein m is equal to [1, number ∈ [ ]]，k∈[1,Ncircle]Ncircle is the number of measurement circles in the particle system;
when d ismk≤RkThen, the mth contact is judged to be located inside the kth measurement circle, where RkTo representThe radius of the kth measurement circle;
calculating the stress tensor of the inner part of the measuring circle according to all the contacts of the inner part of the measuring circle, wherein the calculation formula is as follows:
Figure BDA0003117773050000041
wherein ,
Figure BDA0003117773050000042
representing the stress tensor of the k-th measurement circle in the direction of ij, AkFor the area of the kth measurement circle, NckFor the kth measuring number of contacts inside the circle, F(C)As contact force vector, L(c)A vector connecting the centroids of two particles in contact with each other,
Figure BDA0003117773050000043
is a tensor product, i, j being 1,2 in two dimensions and 1,2,3 in three dimensions;
and calculating the average stress inside the measuring circle according to the stress tensor inside the measuring circle, wherein the calculation formula is as follows:
Figure BDA0003117773050000051
wherein ,
Figure BDA0003117773050000052
represents the mean stress inside the kth measurement circle;
and averaging the average stress in all the measurement circles to obtain the average stress of the particle system, wherein the calculation formula is as follows:
Figure BDA0003117773050000053
wherein ,σsMeans the average stress of the particle system.
Further, step C includes the steps of:
performing iterative calculation on the particle system according to a preset time step by using particle flow software PFC2D/3D, updating contact force of each contact, position and particle speed of each particle in the particle system, and realizing free movement of the particles;
after each free movement of the particles, calculating the difference value between the average stress inside the measuring circle and the target servo stress one by one:
Figure BDA0003117773050000054
wherein ,
Figure BDA0003117773050000055
represents the difference between the average stress inside the kth measurement circle and the target servo stress, stress represents the target servo stress,
Figure BDA0003117773050000056
denotes the mean stress inside the kth measurement circle, k ∈ [1, Ncirle-]Ncircle is the number of measurement circles in the particle system;
according to the difference
Figure BDA0003117773050000057
Calculating the radius change amount of the particles in the kth measurement circle:
Figure BDA0003117773050000058
wherein ,drkRepresenting the radius change of the particles in the kth measurement circle, alpha is an approximation coefficient, and alpha belongs to [0,1 ]]，KnRepresents the normal contact stiffness between the particles;
according to the radius change drkUpdating the size of each particle in the measurement circle, and the formula is as follows:
Rn＝R0n-drk (6)
wherein ,R0nIndicating before updateRadius of particle, RnIndicating the updated particle radius.
Further, in the two-dimensional case, the method for determining whether the particle belongs to the measurement circle is as follows:
traversing all the particles in the particle system, and calculating the distance from the center of each particle to the center of the measuring circle:
Figure BDA0003117773050000061
wherein ,dnRepresents the distance from the center of the nth particle to the center of the kth measuring circle, (x)k,yk) As the center coordinates of the kth measurement circle, (x)0n,y0n) The coordinate of the center of the nth particle is shown, wherein n is 1,2, …, and MP are the number of particles in the particle system;
when d isn＜RkThen, the nth particle is judged to be located inside the kth measuring circle, wherein RkIndicating the radius of the kth measurement circle.
Further, in the two-dimensional case, the step D includes the steps of:
traversing the updated contacts corresponding to all the particles, obtaining particle-wall contacts, and calculating the distance between the center of each particle and the corresponding particle-wall contact:
Figure BDA0003117773050000062
wherein ,d2Denotes the distance between the center of the particle and the particle-wall contact, (x)c,yc) Is the coordinate of the particle-wall contact, (x)1,y1) The updated coordinates of the center of the circle of the particles are obtained;
when d is2＜λ*RnWhen the particle is updated, the particle is considered to be embedded in the boundary of the geometric region, the particle is deleted, otherwise, the particle is kept updated, wherein lambda is an empirical parameter, lambda belongs to (0,1), and R belongs tonRepresents the updated particle radius;
and forming a renewed particle system by using the reserved renewed particles.
Further, step E comprises the steps of:
calculating an error rate based on the updated average stress of the particle system and the target servo stress:
Error＝|σs-stress|/stress (9)
where Error represents the Error rate between the average stress of the updated particle system and the target servo stress, σsRepresenting the average stress of the updated particle system, stress representing the target servo stress;
and when Error is greater than E, repeating the steps C-E, carrying out iterative updating on the particle system, and when Error is less than or equal to E, ending the iteration to obtain an optimal particle system model, wherein E is a preset Error tolerance value.
Further, the method comprises the following steps:
and obtaining an average stress change curve of the particle system in the process of iteratively updating the particle system, and adjusting an approximation coefficient alpha according to the average stress change curve.
The following advantages can be obtained by adopting the technical means:
the invention provides a particle system servo control method without changing boundary shape, which comprises the steps of controlling random generation of particles with certain porosity through the outer contour of a geometric area to obtain an initial particle system, arranging measurement circles in the particle system and calculating average stress, then utilizing a certain amount of time steps to calculate particles to deform freely, continuously adjusting the particle radius in each measurement circle through the difference between the average stress and target servo stress, finally enabling each part of the particle system to approach servo stress to obtain an optimal particle system model, and endowing bonding model parameters on the basis of the optimal particle system model to calculate so as to ensure accurate control of the model state, thereby ensuring the correctness of a numerical calculation result.
Compared with the prior art, the method realizes the amplification or reduction of the average stress in the model by changing the radius of the particles in the measurement circle, and finally iterates to the servo stress without changing the boundary position and the shape of the model all the time; in the servo control process, the particles clamped on the boundary of the geometric area due to particle amplification are directly deleted through contact judgment, and the integral stress control cannot be influenced. The method can quickly and effectively obtain the particle system meeting the requirements, is simple to realize, and is particularly suitable for generating the particle system in a very irregular model range.
Drawings
FIG. 1 is a flow chart illustrating the steps of a servo control method for a grain system without changing the shape of the boundary according to the present invention;
FIG. 2 is a schematic diagram of a geometric region of a two-dimensional multi-segment model according to an embodiment of the present invention;
FIG. 3 is a schematic illustration of particle contact in an embodiment of the present invention;
FIG. 4 is a schematic diagram of a measurement circle in a two-dimensional multi-segment model according to an embodiment of the present invention;
FIG. 5 is a graph showing the average stress variation in example 1 of the present invention;
FIG. 6 is a schematic diagram of a geometric region surrounded by three-dimensional closed triangles according to an embodiment of the present invention;
FIG. 7 is a schematic illustration of a particle system in a three-dimensional model according to an embodiment of the present invention;
FIG. 8 is a schematic view of contacts in a three-dimensional model according to an embodiment of the present invention;
FIG. 9 is a schematic view of a measurement circle in a three-dimensional model according to an embodiment of the present invention;
FIG. 10 is a graph showing the average stress variation in example 2 of the present invention.
Detailed Description
The technical scheme of the invention is further explained by combining the accompanying drawings as follows:
the invention relates to the construction aspect of a disc (sphere) particle system during the numerical simulation pretreatment of a particle flow in the technical field of civil engineering and the like, and mainly provides a particle system servo control method without changing the boundary shape, as shown in figure 1, which specifically comprises the following steps:
step A, generating a plurality of particles in a preset geometric area according to a preset porosity to obtain an initial particle system, wherein the specific operation is as follows: step A01,And carrying out geometric boundary limitation according to the model building requirement to obtain a closed geometric region. Taking a two-dimensional case as an example, a geometric area is shown in fig. 2, and the boundary of the geometric area is defined by U straight line segments (as shown in (a) in fig. 2), each straight line segment can be referred to as a wall, and the intersection point of adjacent walls is referred to as a vertex (as shown in (c) in fig. 2). Traversing coordinates of all vertexes of the geometric area to obtain the minimum value x of the x axial coordinateminMaximum value x of x axial coordinatemaxY minimum value of the axial coordinate yminMaximum value y of y-axis coordinatemax(ii) a Using xmin、xmax、ymin、ymaxA rectangular area can be constructed.
Step A02, presetting a certain initial porosity and a certain minimum radius rminMaximum radius rmaxAnd randomly generating a plurality of particles in a rectangular range according to the radius range and the initial porosity of the particles, wherein the centers of the particles of all the particles are required to be positioned in a closed geometric area. Under the two-dimensional condition, the initial porosity value ranges from 0.15 to 0.20, and under the three-dimensional condition, the initial porosity value ranges from 0.3 to 0.4.
And A03, performing position judgment on each particle by using a ray method, specifically, traversing all the particles generated in the step A02, leading out a ray from the center of the particle to the right side of each particle, wherein the ray direction vector is (1,0), and acquiring the number of intersection points of the ray and the boundary (U line segments) of the geometric region.
A04, when the number of the intersection points corresponding to a particle is 0 or even, deleting the particle; when the number of the intersection points corresponding to one particle is an odd number, the particle is reserved; the remaining particles are shown as (b) in fig. 2, and all lie within the boundary of the geometric region, and form an initial particle system, and the number of particles in the initial particle system is represented as MP.
And B, setting regular measuring circles (two-dimensional are circular discs and three-dimensional are round balls) in the initial particle system according to the particle size of the generated particles, and calculating the average stress in each measuring circle and the average stress of the particle system.
Step B01 obtaining a measurement based on the maximum particle radius in the initial particle systemThe radius of the circle is a circular area for testing a stress field in a small range, and in the method, the value range of the radius R of the circle is 5-10 times of the radius of the maximum particle, namely [5 x Rmax,10*rmax], wherein ,rmaxThe maximum particle radius.
And step B02, obtaining a rectangular range corresponding to the geometric area according to the coordinates of the geometric area, arranging a plurality of tangent measuring circles in the rectangular range according to the radius of the measuring circle, and simultaneously ensuring that more than 20 contacts exist in each measuring circle.
In the two-dimensional case, the rectangular range can be represented as [ x ]min,xmax]，[ymin，ymax]Assuming that the measuring circles are tangent circular areas which are regularly arranged horizontally and vertically, the number of the measuring circles along the x-axis direction is nx＝(xmax-xmin) V (2.0R), the number of circles measured along the y-axis being ny＝(ymax-ymin) V (2.0 × R) measurement circles. N is required to be placed in the rectangular rangex*nyThe coordinates (x) of the center of the circle of the measurement circle of the v-th row and the w-th columnvw,yvw) Can be expressed as:
xvw＝xmin+R+(w-1)*2.0*R (10)
yvw＝ymin+R+(v-1)*2.0*R (11)
and step B03, traversing all the measuring circles, extending a ray from the center of each measuring circle to the vector (1,0) in the right direction, and acquiring the number of intersection points of the ray and the boundary of the geometric area.
And step B04, when the number of the intersection points corresponding to one measuring circle is 0 or even, deleting the measuring circle, otherwise, keeping the measuring circle, and as shown in FIG. 4, keeping the rest measuring circles within the boundary limit of the geometric area.
And step B05, judging the particle-particle contact and the particle-wall contact in the particle system according to the distance between the centers of the particles, the distance between the centers of the particles and the boundary of the geometric area and the particle radius, and obtaining the total number of the contact in the particle system. When the distance between the centers of any two particles is less than the sum of the radii of the two particles, it means that there is a contact between the two particles, i.e. particle-particle contact; for any particle, if the distance between the center of the particle and any boundary of the geometric area is less than or equal to the radius of the particle, it indicates that there is a contact between the particle and the boundary of the geometric area, i.e. particle-wall contact.
The interparticle (particle-particle, particle-wall) effects were set following a linear contact model, which included two parameters, normal stiffness, tangential stiffness, with a particle density of 1000kg/m3 by default. As shown in fig. 3, the normal stiffness describes the relationship between the total normal force and the total normal displacement, and the tangential stiffness describes the relationship between the shear stress increment and the shear displacement increment, and the specific expression is as follows:
Figure BDA0003117773050000111
wherein ,
Figure BDA0003117773050000112
denotes the total normal force, k, of the m-th contactnIndicating normal stiffness, UnDenotes the normal displacement of the contact, nmThe normal vector representing the m-th touch,
Figure BDA0003117773050000113
denotes the shear stress increment, k, of the m-th contactsThe tangential stiffness is expressed in terms of the stiffness,
Figure BDA0003117773050000114
denotes the shear displacement increment, m ∈ [1, number]。
Step B06, traversing all contacts in the particle system, and calculating the distance d from the m contact to the k measuring circlemk：
Figure BDA0003117773050000115
wherein ,(xm,ym) Is the m < th > oneCoordinates of contact, (x)k,yk) As the center coordinates of the kth measurement circle, k ∈ [1, Ncirle ]]Ncircle is the number of measurement circles in the particle system.
Step B07, when dmk≤RkThen, the mth contact is judged to be located inside the kth measurement circle, where RkIndicating the radius of the kth measurement circle. The contact number N in the kth measuring circle can be obtained by traversing all the contactsckIn the present invention, the radii of all measurement circles are equal.
Step B08, calculating the stress tensor of the measuring circle according to all the contacts in the measuring circle, wherein the specific calculation formula is as follows:
Figure BDA0003117773050000121
wherein ,
Figure BDA0003117773050000122
representing the stress tensor of the k-th measurement circle in the direction of ij, AkFor the area of the kth measurement circle,
Figure BDA0003117773050000123
F(C)as contact force vector, L(c)A vector connecting the centroids of two particles in contact with each other,
Figure BDA0003117773050000124
is a tensor product and specifies a negative value for the compressive stress. In two dimensions, i, j is 1,2, which respectively represent x and y coordinate directions in a two-dimensional coordinate system; in three dimensions, i, j is 1,2, and 3, and represents x, y, and z coordinate directions in a three-dimensional coordinate system, respectively. When the value of i is equal to j,
Figure BDA0003117773050000125
for a positive stress, when i is not equal to j,
Figure BDA0003117773050000126
is a shear stress.
Step B09, calculating the average stress inside the measurement circle according to the stress tensor inside the measurement circle, and only considering the average of the normal stress when calculating the average stress inside the measurement circle by using the measurement circle, therefore, the specific calculation formula is as follows:
Figure BDA0003117773050000127
wherein ,
Figure BDA0003117773050000128
the mean stress inside the kth measurement circle is indicated.
Step B10, averaging the average stress in all the measuring circles to obtain the average stress of the particle system, wherein the specific calculation formula is as follows:
Figure BDA0003117773050000129
wherein ,σsMeans the average stress of the particle system.
And C, controlling the free movement of the particles in the particle system by using particle flow software, and updating the size of the particles in each measuring circle according to the average stress in the measuring circle in the free movement process of the particles.
And step C01, enabling the particle system to freely deform and move by utilizing a mature particle flow program PFC, setting a certain time step, such as 10000 steps, utilizing particle flow software PFC2D/3D, carrying out iterative calculation on the particle system according to the preset time step, and updating the contact force of each contact (particle-particle, particle-wall) in the particle system, the position of each particle and the particle speed, so that the stress field of the particle system is changed, the stress gradually tends to be stable, and the free movement of the particles is realized.
Step C02, after each particle free motion, calculating the average stress inside the measuring circle one by one according to the method in the step B, and calculating the difference value between the average stress inside each measuring circle and the target servo stress, wherein the specific calculation formula is as follows:
Figure BDA0003117773050000131
wherein ,
Figure BDA0003117773050000132
the difference between the average stress inside the kth measurement circle and the target servo stress is shown, and stress represents the target servo stress.
Step C03, according to the difference
Figure BDA0003117773050000133
Calculating the radius change amount of the particles in the kth measurement circle:
Figure BDA0003117773050000134
wherein ,drkExpressing the radius change of particles in the kth measurement circle, wherein alpha is an approximation coefficient, the value range of alpha is 0-1, and KnIndicating the normal contact stiffness between the particles.
Step C04, the radius change quantity of all the particle radiuses in the same measuring circle is assumed to be completely the same, and the change quantity d is determined according to the radius change quantityrkAnd updating the size of each particle in the measurement circle, wherein the specific formula is as follows:
Rn＝R0n-drk (19)
wherein ,R0nDenotes the particle radius before renewal, RnIndicating the updated particle radius.
Taking the two-dimensional case as an example, the method for determining whether the particle belongs to the measurement circle in step C03 is as follows:
(1) traversing all the particles in the particle system, and calculating the distance from the center of each particle to the center of the measuring circle:
Figure BDA0003117773050000141
wherein ,dnRepresents the distance from the center of the nth particle to the center of the kth measuring circle, (x)k,yk) As the center coordinates of the kth measurement circle, (x)0n,y0n) The coordinate of the center of the nth particle is n, 1,2, … and MP.
(2) When d isn＜RkThen, the nth particle is judged to be located inside the kth measurement circle.
And D, judging and deleting edge particles according to the updated particle size to obtain an updated particle system. In the two-dimensional case, the specific operation of step D is as follows:
step D01, traversing the updated contacts corresponding to all the particles, judging the contact type, obtaining the particle-wall contact, and extracting the coordinates (x) of the particle-wall contactc,yc) And calculating the distance between the center of each particle and its corresponding particle-wall contact:
Figure BDA0003117773050000142
wherein ,d2Denotes the distance between the center of the particle and the particle-wall contact, (x)1,y1) The updated coordinates of the center of the particle circle.
Step D02, when D2＜λ*RnWhen the particle is considered to be embedded in the boundary of the geometric region after updating, the particle is deleted, otherwise the particle after updating is reserved, wherein, lambda is an empirical parameter, and lambda epsilon (0,1) is preferably 0.9.
And D03, forming a renewed particle system by using the reserved renewed particles.
And E, calculating the error rate between the updated average stress of the particle system and the target servo stress, and repeating the steps C to E according to the error rate to carry out iterative calculation to obtain an optimal particle system model. The specific operation of step E is as follows:
and E01, respectively updating the particle radius of the Ncircle measurement circle positions and deleting unreasonable boundary particles, then calculating a certain time step (such as 10000 steps) by utilizing mature particle flow software PFC, enabling the particle system to be basically balanced, and extracting the unbalanced force of the updated particle system. And B, calculating the average stress of all the measuring circles by using the updated particles according to the method in the step B, and obtaining the average stress of the updated particle system.
Calculating an error rate based on the updated average stress of the particle system and the target servo stress:
Error＝|σs-stress|/stress (22)
where Error represents the Error rate between the average stress of the updated particle system and the target servo stress.
And E02, when Error is greater than E, repeating the steps C-E, iteratively updating the particle system until the Error is less than or equal to E and the unbalance force of the system is small enough (default is 1E-5), and stopping iteration to obtain an optimal particle system model, wherein E is a preset Error tolerance value, and the value range of E is 0.001-0.05.
And step F, acquiring an average stress change curve of the particle system in the process of iteratively updating the particle system, wherein if the curve gradually approaches to the target servo stress, the process is correct, and the value of the approximation coefficient can be controlled according to the change amplitude of the average stress. In the embodiment of the present invention, the default approximation coefficient α is 0.1, when the average stress obtained by updating the particle system in each iteration is farther from the target servo stress, the approximation coefficient α may be appropriately increased, and when the average stress obtained by updating the particle system in each iteration is closer to the target servo stress, the approximation coefficient α may be appropriately decreased.
And when Error is less than or equal to e, obtaining the radius of each particle (a two-dimensional disc and a three-dimensional sphere), wherein the range (the minimum radius size-the maximum radius size) of all particle radii is the optimal particle radius range of the model.
To verify the effectiveness of the method of the invention, 2 specific examples are given below:
example 1:
when a complex two-dimensional model with uneven profile is built, as shown in FIG. 2, a compact particle system is built to simulate a strongly-bonded rock mass, and the model consists of 10 sections of two-dimensional wallsIs enclosed to form an area of 1336m2. The inside is desirably filled with particles and contacted with each other, and the servo mean stress is designed to be-1 MPa (the minus sign indicates compressive stress). The steps of establishing the model by applying the method of the invention are as follows:
(1) traversing the x and y coordinates of all the vertexes to obtain the minimum value x of the x-axis coordinate of the modelmin7.58m, x axial maximum xmax76.91m, y axial minimum value ymin-6.76m, y axial maximum ymax46.51 m; setting the minimum radius r of the random particlesmin0.2m, maximum radius rmax0.3m, 0.17 for initial porosity and 5636 for particle number MP.
(2) Setting a stiffness parameter kn＝1e8，ksAnd (3) 7, wherein the time step of the equilibrium iteration is 10000 steps, so that the particle system reaches the initial equilibrium, and the total contact number in the particle system after the equilibrium is 10696. The target servo stress value stress is set to-1.00 MPa, and the error rate e is calculated to 0.001.
(3) Obtaining the radius R of the measured circle from the radius of the particles as 8.0Rmax. The number of the measurement circles is set along the x-axis and the y-axis, and the measurement circles outside the model are removed by the ray method, so that the distribution of the measurement circles is shown in fig. 4. And circulating each measuring circle, respectively calculating the average stress, the difference value with the target servo stress and the radius change quantity, and updating the particle radius in the measuring circle.
(4) And judging and deleting edge grains according to the updated grain size, calculating an error rate between the error rate and the target servo stress, and performing iterative updating according to the error rate and the error rate e until the requirement is met. The actual run in example 1 was exited after 20 iterations.
(5) The average stress variation curve is recorded, as shown in fig. 5, and it can be seen from the graph that the curve gradually approaches the target value, which satisfies the requirement. If α is 0.3, the iteration step needs only 10 times to converge, and if α is 0.8, the iteration step needs only 6 times to converge. The radius range of the optimal particle meeting the servo stress is 0.1823-0.3069 m, the unbalanced force is smaller than 1e-5, and the average stress (servo stress) of the optimal particle system model is-0.999 MPa.
Example 2:
the cubic capacity of a certain slip mass is 5000 thousands of cubes, a geometric space enclosed by triangles is shown in figure 6, the design is filled by three-dimensional spheres with the radius of 2.8m-4.0m, and the thickness (vertical direction) of the boundary position gradually transits from 0m to nearly 200m inside the slip mass. To establish a compact particle system simulation bonding rock-soil body, the design adopts-0.05 MPa confining pressure (negative sign represents pressure stress) to promote the particle system to tend to be uniform. The steps of establishing the model by applying the method of the invention are as follows:
(1) traversing the x, y and z coordinates of all the vertexes to obtain the minimum value x of the x axial coordinate of the modelminAxial maximum x of-254 m, xmax1511m, y axial minimum ymin448m, y axial maximum value ymax846m, z axial minimum zmin2146m, z axial maximum zmax3225 m; setting the minimum radius r of the random particlesmin2.8m, maximum radius rmaxThe initial particle porosity was taken to be 0.38 and the number of particles MP was 208263. All the particles are traversed, and the particles in the model are judged by using a ray method to obtain a particle system, as shown in fig. 7.
(2) Setting the Normal stiffness kn＝2e8，ksThe time step of the equilibrium iteration is 10000 steps, the particle system reaches the initial equilibrium, and the total number of contacts in the equilibrium system is 782174, as shown in fig. 8. The target servo stress value stress is set to-0.05 MPa, and the error rate e is calculated to 0.01.
(3) The measured circle radius R is 5.0 rmax calculated from the particle radius. The number of the measurement circles is set along the directions of the x axis, the y axis and the z axis, and the measurement circles positioned outside the model are judged and removed by the ray method, so that the distribution of the measurement circles is shown in fig. 9, and the number of the measurement circles is 255. And circulating each measuring circle, respectively calculating the average stress, the difference value with the target servo stress and the radius change quantity, and updating the particle radius in the measuring circle.
(4) Traversing all updated particles, judging all contacts surrounding the particles, two or more particle-wall contacts in three dimensions, according to d2＜0.9R0nJudging and deleting edge particles, and deleting boundary ballsAbout 1000.
(5) Error rates between the average stress target servo stresses are calculated from the updated particle system, and iterative updates are made according to the error rates and e until the requirements are met. The actual run in example 2 was exited after 9 iterations.
(6) The average stress variation curve of the three-dimensional model is recorded, as shown in fig. 10, it can be seen that the curve gradually approaches the target value, the requirement is met, and the model servo is completed. The optimal particle radius range satisfying the servo stress is finally obtained to be 2.8146M-4.0146M, the unbalanced force is smaller than 1e-5 at the moment, and the average stress (servo stress) of an optimal particle system model is-0.499M Pa.
Compared with the prior art, the method realizes the amplification or reduction of the average stress in the model by changing the radius of the particles in the measurement circle, and finally iterates to the servo stress without changing the boundary position and the shape of the model all the time; in the servo control process, the particles clamped on the boundary of the geometric area due to particle amplification are directly deleted through contact judgment, and the integral stress control cannot be influenced. The method can quickly and effectively obtain the particle system meeting the requirements, and the adhesion model parameters are given to calculate on the basis of the optimal particle system model to ensure the accurate control of the model state, thereby ensuring the correctness of the numerical calculation result, and being particularly suitable for the generation of the particle system in a very irregular model range.
The above description is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, several modifications and variations can be made without departing from the technical principle of the present invention, and these modifications and variations should also be regarded as the protection scope of the present invention.

Claims (10)

1. A particle system servo control method without changing boundary shape is characterized by comprising the following steps:
a, generating a plurality of particles in a preset geometric area according to a preset porosity to obtain an initial particle system;
b, setting regular measuring circles in an initial particle system according to the particle size of the generated particles, and calculating the average stress inside each measuring circle and the average stress of the particle system;
c, controlling free movement of particles in a particle system by using particle flow software, and updating the size of the particles in each measuring circle according to the average stress in the measuring circle in the free movement process of the particles;
d, judging and deleting edge particles according to the updated particle size to obtain an updated particle system;
and E, calculating the error rate between the updated average stress of the particle system and the target servo stress, and repeating the steps C to E according to the error rate to carry out iterative calculation to obtain an optimal particle system model.
2. The servo control method of claim 1, wherein the step A comprises the steps of:
carrying out geometric boundary limitation according to the model building requirement to obtain a closed geometric region;
generating a plurality of particles according to a preset particle radius range and initial porosity, wherein the centers of the particles of all the particles are located in a closed geometric area;
judging the position of each particle by using a ray method, leading out a ray from the center of the particle of each particle to the right side, and acquiring the number of intersection points of the ray and the boundary of a geometric region, wherein the direction vector of the ray is (1, 0);
when the number of the intersection points corresponding to one particle is 0 or even, deleting the particle, otherwise, keeping the particle;
the retained particles are used to form the initial particle system.
3. The servo control method for a particle system without changing the boundary shape according to claim 1, wherein in the step B, the method of setting the regular measurement circle in the initial particle system according to the particle size of the generated particle comprises:
obtaining a measurement circle radius according to the maximum particle radius in the initial particle system;
obtaining a rectangular range corresponding to the geometric area according to the coordinates of the geometric area, and setting a plurality of tangent measuring circles in the rectangular range according to the radius of the measuring circles;
extending a ray from the center of each measuring circle to a vector (1,0) in the right direction, and acquiring the number of intersection points of the ray and the boundary of the geometric area;
and when the number of the intersection points corresponding to one measuring circle is 0 or even, deleting the measuring circle, otherwise, keeping the measuring circle.
4. A method as claimed in claim 3, wherein the radius of the measurement circle is in the range of [5 x r [ ]max,10*rmax], wherein ,rmaxThe maximum particle radius.
5. The servo control method for a particle system without changing the boundary shape according to claim 1, wherein in the step B, the calculation method of the average stress inside each measurement circle and the average stress of the particle system comprises the following steps:
judging the particle-particle contact and the particle-wall contact in the particle system according to the distance between the centers of the particles, the distance between the centers of the particles and the boundary of the geometric area and the particle radius, and acquiring the total number of the contact in the particle system;
traversing all contacts in the particle system, calculating the distance d from the m-th contact to the k-th measurement circlemkWherein m is equal to [1, number ∈ [ ]]，k∈[1,Ncircle]Ncircle is the number of measurement circles in the particle system;
when d ismk≤RkThen, the mth contact is judged to be located inside the kth measurement circle, where RkRepresents the radius of the kth measurement circle;
calculating the stress tensor of the inner part of the measuring circle according to all the contacts of the inner part of the measuring circle, wherein the calculation formula is as follows:
Figure FDA0003117773040000031
wherein ,
Figure FDA0003117773040000032
representing the stress tensor of the k-th measurement circle in the direction of ij, AkFor the area of the kth measurement circle, NckFor the kth measuring number of contacts inside the circle, F(C)As contact force vector, L(c)A vector connecting the centroids of two particles in contact with each other,
Figure FDA0003117773040000033
is a tensor product, i, j being 1,2 in two dimensions and 1,2,3 in three dimensions;
and calculating the average stress inside the measuring circle according to the stress tensor inside the measuring circle, wherein the calculation formula is as follows:
Figure FDA0003117773040000034
wherein ,
Figure FDA0003117773040000035
represents the mean stress inside the kth measurement circle;
and averaging the average stress in all the measurement circles to obtain the average stress of the particle system, wherein the calculation formula is as follows:
Figure FDA0003117773040000041
wherein ,σsMeans the average stress of the particle system.
6. The servo control method of claim 1, wherein step C comprises the steps of:
performing iterative calculation on the particle system according to a preset time step by using particle flow software PFC2D/3D, updating contact force of each contact, position and particle speed of each particle in the particle system, and realizing free movement of the particles;
after each free movement of the particles, calculating the difference value between the average stress inside the measuring circle and the target servo stress one by one:
Figure FDA0003117773040000042
wherein ,
Figure FDA0003117773040000043
represents the difference between the average stress inside the kth measurement circle and the target servo stress, stress represents the target servo stress,
Figure FDA0003117773040000044
denotes the mean stress inside the kth measurement circle, k ∈ [1, Ncirle-]Ncircle is the number of measurement circles in the particle system;
according to the difference
Figure FDA0003117773040000045
Calculating the radius change amount of the particles in the kth measurement circle:
Figure FDA0003117773040000046
wherein ,drkRepresenting the radius change of the particles in the kth measurement circle, alpha is an approximation coefficient, and alpha belongs to [0,1 ]]，KnRepresents the normal contact stiffness between the particles;
according to the radius change drkUpdating the size of each particle in the measurement circle, and the formula is as follows:
Rn＝R0n-drk
wherein ,R0nDenotes the particle radius before renewal, RnIndicating the updated particle radius.
7. The servo control method of claim 6, wherein the method for determining whether the particle belongs to the circle under two dimensions comprises:
traversing all the particles in the particle system, and calculating the distance from the center of each particle to the center of the measuring circle:
Figure FDA0003117773040000051
wherein ,dnRepresents the distance from the center of the nth particle to the center of the kth measuring circle, (x)k,yk) As the center coordinates of the kth measurement circle, (x)0n,y0n) The coordinate of the center of the nth particle is shown, wherein n is 1,2, …, and MP are the number of particles in the particle system;
when d isn＜RkThen, the nth particle is judged to be located inside the kth measuring circle, wherein RkIndicating the radius of the kth measurement circle.
8. The servo control method of claim 5, wherein step D comprises the following steps in two dimensions:
traversing the updated contacts corresponding to all the particles, obtaining particle-wall contacts, and calculating the distance between the center of each particle and the corresponding particle-wall contact:
Figure FDA0003117773040000052
wherein ,d2Denotes the distance between the center of the particle and the particle-wall contact, (x)c,yc) Is the coordinate of the particle-wall contact, (x)1,y1) To be updatedThe coordinates of the center of the particle;
when d is2＜λ*RnWhen the particle is updated, the particle is considered to be embedded in the boundary of the geometric region, the particle is deleted, otherwise, the particle is kept updated, wherein lambda is an empirical parameter, lambda belongs to (0,1), and R belongs tonRepresents the updated particle radius;
and forming a renewed particle system by using the reserved renewed particles.
9. The servo control method of claim 1, wherein step E comprises the steps of:
calculating an error rate based on the updated average stress of the particle system and the target servo stress:
Error＝|σs-stress|/stress
where Error represents the Error rate between the average stress of the updated particle system and the target servo stress, σsRepresenting the average stress of the updated particle system, stress representing the target servo stress;
and when Error is greater than E, repeating the steps C-E, carrying out iterative updating on the particle system, and when Error is less than or equal to E, ending the iteration to obtain an optimal particle system model, wherein E is a preset Error tolerance value.
10. The method of claim 6, further comprising the steps of:
and obtaining an average stress change curve of the particle system in the process of iteratively updating the particle system, and adjusting an approximation coefficient alpha according to the average stress change curve.
CN202110668184.9A 2021-06-16 2021-06-16 Granular system servo control method without changing boundary shape Active CN113408124B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110668184.9A CN113408124B (en) 2021-06-16 2021-06-16 Granular system servo control method without changing boundary shape

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110668184.9A CN113408124B (en) 2021-06-16 2021-06-16 Granular system servo control method without changing boundary shape

Publications (2)

Publication Number Publication Date
CN113408124A true CN113408124A (en) 2021-09-17
CN113408124B CN113408124B (en) 2023-08-22

Family

ID=77684425

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110668184.9A Active CN113408124B (en) 2021-06-16 2021-06-16 Granular system servo control method without changing boundary shape

Country Status (1)

Country Link
CN (1) CN113408124B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114757119A (en) * 2022-04-01 2022-07-15 河海大学 Two-dimensional hydraulic calculation method for constructing pipe network by using outer-wrapped polygon
CN115795779A (en) * 2022-09-13 2023-03-14 信阳水投引九济石工程管理有限公司 Method for tracking crack formation process through measuring point displacement

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011110755A2 (en) * 2010-03-08 2011-09-15 Ronald Tai Method for deterministic modeling of powders and systems of particles
CN110688748A (en) * 2019-09-19 2020-01-14 湘潭大学 Single-particle discrete element numerical sample modeling method with random shape
CN111610091A (en) * 2020-05-11 2020-09-01 太原理工大学 Automatic calibration method for discrete element Hertz contact parameter during simulation of geotechnical material
CN112231963A (en) * 2020-10-28 2021-01-15 温州大学 Discrete element model modeling method considering temperature change and application thereof
CN112417715A (en) * 2020-10-26 2021-02-26 山东大学 Rock mass fracture simulation method and system under true triaxial servo loading state

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2011110755A2 (en) * 2010-03-08 2011-09-15 Ronald Tai Method for deterministic modeling of powders and systems of particles
CN110688748A (en) * 2019-09-19 2020-01-14 湘潭大学 Single-particle discrete element numerical sample modeling method with random shape
CN111610091A (en) * 2020-05-11 2020-09-01 太原理工大学 Automatic calibration method for discrete element Hertz contact parameter during simulation of geotechnical material
CN112417715A (en) * 2020-10-26 2021-02-26 山东大学 Rock mass fracture simulation method and system under true triaxial servo loading state
CN112231963A (en) * 2020-10-28 2021-01-15 温州大学 Discrete element model modeling method considering temperature change and application thereof

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114757119A (en) * 2022-04-01 2022-07-15 河海大学 Two-dimensional hydraulic calculation method for constructing pipe network by using outer-wrapped polygon
CN114757119B (en) * 2022-04-01 2023-08-08 河海大学 Two-dimensional hydraulic calculation method for constructing pipe network by utilizing outsourcing polygons
CN115795779A (en) * 2022-09-13 2023-03-14 信阳水投引九济石工程管理有限公司 Method for tracking crack formation process through measuring point displacement
CN115795779B (en) * 2022-09-13 2024-02-23 信阳水投引九济石工程管理有限公司 Method for tracking crack formation process through measuring point deflection

Also Published As

Publication number Publication date
CN113408124B (en) 2023-08-22

Similar Documents

Publication Publication Date Title
CN113408124B (en) Granular system servo control method without changing boundary shape
CN106960070B (en) Seepage simulation method for reconstructing coal body based on finite element-discrete element CT
Falco et al. Generation of 3D polycrystalline microstructures with a conditioned Laguerre-Voronoi tessellation technique
WO2019010859A1 (en) Modeling method for high-compactness discrete particle heterogeneous system
Ruiz-Gironés et al. Generation of curved high-order meshes with optimal quality and geometric accuracy
CN110457785B (en) Material information mapping method for material point method of structural large deformation response
CN106650018B (en) High volume fraction RVE model generation method for VCFEM analysis
WO2021186723A1 (en) Lamination modeling route generation device, lamination modeling route generation method, and machine learning device
CN116227373A (en) Concrete mechanical arm operation simulation optimization method and system based on DEM-CFD
CN109522589B (en) Non-analytic method for simulating two-phase flow of pipeline particles and electronic equipment
CN114091225A (en) Two-dimensional arbitrary-shape aggregate particle discrete element numerical sample modeling method
CN112100890B (en) Calculation method of fracture horizontal seam initiation expansion full-three-dimensional mathematical model
Alauzet Efficient moving mesh technique using generalized swapping
US20070150231A1 (en) Method and apparatus for implementing multi-grid computation for multi-cell computer models with embedded cells
CN115758766A (en) Continuous-discontinuous numerical model connecting method
CN115130299A (en) Method for designing profile curve of single-pivot flexible-wall spray pipe
CN114757119A (en) Two-dimensional hydraulic calculation method for constructing pipe network by using outer-wrapped polygon
CN104847414B (en) Structured dynamic mesh modeling method for vortex type fluid machine
Hattangady Automatic remeshing in 3‐D analysis of forming processes
CN113695449A (en) Gas cylinder hot spinning machining method based on finite element modeling
CN116663336B (en) Electromagnetic modeling method based on cylindrical side surface conformal grid generation
CN111310277A (en) Modeling method for pipeline transfer characteristics of atmospheric data sensing system, aircraft and storage medium
CN113886964B (en) Aircraft pneumatic analysis method based on boundary substitution hybrid model
CN113297781B (en) Particle cluster modeling method based on internal acting force
CN115906714B (en) Plate PECVD reaction chamber airflow simulation method and simulation equipment thereof

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant