Research method, program and equipment for realizing solid deformable crushing based on LBM-DEM coupling algorithm, and storage medium

Abstract

According to the method, a fluid computational domain is divided into grids at equal intervals by adopting an LBM method, each grid is endowed with nine directions, and all the grids in the computational domain are subjected to collision and migration operation, so that the flow of a flow field is simulated; the method comprises the following steps: dividing a solid structure into rectangular units with the same size by adopting a DEM method, performing collision or bonding state identification on each rectangular unit, and performing interaction force solution calculation so as to simulate deformation and damage of the solid structure; an immersion moving boundary method is used for processing the boundary of fluid and solid, transmitting the acting force between the fluid and the solid structure, identifying all grid division node types in a computational domain, calculating the volume fraction proportion of nodes of the boundary type, and then calculating and solving the acting force of the fluid on the solid structure. The mechanical property of the breakable solid material in the fluid can be correctly described, and the method can be used for solving related engineering problems in the fluid-solid coupling field and the like when the solid is deformed or damaged.

Description

Technical Field

[0001] The present invention belongs to the technical field of fluid mechanics numerical calculation, and specifically relates to a research method, program, device and storage medium for realizing solid deformable crushing based on LBM-DEM coupling algorithm. Background Art

[0002] Computer numerical simulation is a comprehensive application technology, which has great application value in teaching, scientific research, design, production, management, decision-making and other departments. The cost of penetration, explosion and other tests is extremely expensive and there are certain risks. However, the use of numerical simulation not only has great economic benefits, but also can accelerate the progress of theoretical and experimental research.

[0003] In the fluid-solid coupling problem, fluid-solid two-phase flow is a mixed fluid flow composed of solid and gas or liquid phases, which is widely used in thermal power, chemical industry, metallurgy, water conservancy and marine engineering. In order to reduce the test cost and predict the interaction between the two phases, we can provide more reliable theoretical guidance for the scientific formulation of the test plan, the optimal location of the measuring point during the test, the determination of the instrument range, etc.

[0004] In the current field of fluid-solid coupling, there are a series of difficulties in traditional numerical simulation, such as moving boundaries and the realization of solid phase freedom based on the assumption of continuous media. We use the mesoscopic physical background of the lattice Boltzmann (LBM) method, particle characteristics and the characteristics of discrete element (DEM) nonlinear continuous media to deal with fluid-solid coupling problems, but the existing LBM-DEM coupling algorithm technology can only simulate the movement and collision of simple single or multiple rigid solids in the fluid, but the solid cannot be deformed or destroyed, and it is impossible to characterize the properties and movement trends of breakable materials in the fluid. Summary of the invention

[0005] The purpose of the present invention is to provide a research method, program, equipment and storage medium for realizing solid deformable crushing based on LBM-DEM coupling algorithm.

[0006] The research method for realizing solid deformable crushing based on LBM-DEM coupling algorithm includes the following steps:

[0007] Step 1: For the problem of movement, collision and deformable breakage of a long strip solid in a liquid-filled tube, determine the total prediction time, obtain the parameter information of the fluid and the long strip solid, and construct a two-dimensional fluid calculation domain. The upper and lower boundaries of the fluid calculation domain are the upper and lower walls of the liquid-filled tube, and have a left boundary and a right boundary; determine the initial layout of the long strip solid in the fluid calculation domain;

[0008] Step 2: Divide the long strip solid into multiple segments evenly in the length direction; divide the fluid calculation domain evenly into grids, each grid is the same square grid, and determine the calculation direction of each grid;

[0009] Step 3: Calculate the time step of the LBM main cycle and the duration of the DEM secondary cycle based on the parameter information of the fluid and the long strip solid;

[0010] Step 4: Initialize the segments of the long strip solid to be in a bonding state, and the speed, rotation angular velocity, displacement and azimuth of each long strip solid segment are all zero; initialize the distribution vector of each grid in each calculation direction in the fluid calculation domain;

[0011] Step 5: Execute the LBM main loop and use the immersed moving boundary method to process the boundaries of the fluid and the solid. Except for the boundary grids of the fluid calculation domain, perform collision operations on other grids. For each grid that performs collision operations, calculate the collision term based on the distribution vectors of each calculation direction of the grid. The collision term of the boundary grid of the fluid calculation domain is a zero vector. Calculate the total force of the fluid on the long strip solid and the torque on each long strip solid segment based on the collision term. Finally, perform migration operations on each grid in the fluid calculation domain and update the distribution vectors of each calculation direction of each grid.

[0012] Step 6: Execute the DEM sub-cycle until the duration of the DEM sub-cycle is reached;

[0013] In each iteration of the DEM cycle, each long strip solid segment is traversed to obtain the contact conditions between each long strip solid segment and form a contact list; for each combination of two long strip solid segments in the contact list, it is determined whether they are in a bonding state; if two long strip solid segments are determined to be out of the bonding state in a certain calculation, then in subsequent calculations, the two will not be re-bonded with any long strip solid segment;

[0014] If it is judged that the two long strip solid segments are in a bonding state, the limit normal stress at the contact surface of the two segments is calculated, and it is judged whether the limit normal stress is less than the limit tensile strength and the limit compressive strength of the long strip solid material; if it is less than, it is judged that no fracture occurs between the two long strip solid segments, and the force vector and torque received by each long strip solid segment are calculated; otherwise, it is judged that a fracture occurs between the two long strip solid segments, and there is no force vector and torque between the two;

[0015] If there is an overlapping area between the two long strip solid segments, it is determined that a collision occurs between the two, and the force vector and torque received by each long strip solid segment are calculated;

[0016] For each long strip solid segment, according to the total force of the fluid on the long strip solid, the torque of the fluid on the long strip solid segment, and the force vector and torque of the long strip solid segment from other long strip solid segments, update the speed, angular velocity of rotation, displacement in the x-axis direction, displacement in the y-axis direction and azimuth of the long strip solid segment;

[0017] Step 7: Repeat steps 5 to 6 according to the time step of the LBM main loop until the total prediction time is reached, output the calculation results, and complete the prediction of the movement, collision and deformable breakage of the long solid in the liquid-filled tube.

[0018] Furthermore, in step 1, the total prediction time T is determined, and parameter information of the fluid and the long strip solid is obtained, including the flow rate u of the fluid. 0 , density ρ 0 , viscosity visco and relaxation factor τ, the length of the long solid L 1 、Width L 2 , density ρ s , the elastic modulus E, tangential modulus G, Poisson's ratio poisson, and normal stiffness coefficient k of the long strip solid material n , tangential stiffness coefficient k t , normal viscosity coefficient k nv , tangential viscosity coefficient k tv , damping coefficient η, friction coefficient μ, ultimate tensile strength, ultimate compressive strength.

[0019] Furthermore, in step 2, the long strip solid is evenly divided into N s The mass of each long solid segment is m Ns The moment of inertia is I Ns ; Take the intersection of the left boundary and the lower boundary of the fluid calculation domain as the origin, take the lower boundary as the x-axis, the positive direction of the x-axis is the flow direction of the fluid, take the y-axis of the left boundary, and establish a plane rectangular coordinate system; divide the fluid calculation domain into grids, each grid is a square grid with a side length of dx, and the index of each grid is (a x ,a y ), represents the ath x Row a y A grid of columns, a x =1,2,...,N x , a y =1,2,...,N y , N x With N y are the number of rows and columns of the grid respectively, and the total area of ​​the fluid calculation domain is N x ·dx·N y ·dx;

[0020] For each grid, take its center point as the 0th direction, the vertical line from the center point to the right side of the grid as the 1st direction, the vertical line from the center point to the left side of the grid as the 2nd direction, the vertical line from the center point to the top of the grid as the 3rd direction, the vertical line from the center point to the bottom of the grid as the 4th direction, the vertical line from the center point to the upper right vertex of the grid as the 5th direction, the vertical line from the center point to the upper left vertex of the grid as the 6th direction, the vertical line from the center point to the lower left vertex of the grid as the 7th direction, and the vertical line from the center point to the lower right vertex of the grid as the 8th direction, and construct the vector e i for:

[0021]

[0022] Furthermore, the method for calculating the time step dt of the LBM main cycle and the time step dtime of the DEM secondary cycle in step 3 is specifically:

[0023]

[0024] in,

[0025] dtime 0 Correction is made so that the corrected dtime satisfies dt=N 2 dtime, N 2 Is a positive integer, and dtime≤dtime 0 .

[0026] Furthermore, in step 4, t=0 is initialized, and the segments of the long strip solid are initialized to be in a bonding state, and the speed, rotational angular velocity, displacement and azimuth of each long strip solid segment are all 0, that is, ω dtime (j,0)=0, θ 0 (j,0) = 0; Initialize the distribution vector f in 9 directions of each grid in the fluid calculation domain i (a x ,a y ,0) and the equilibrium distribution vector

[0027] in,

[0028] Furthermore, the LBM main loop is executed in step 5, which specifically includes the following steps:

[0029] Step 5.1: Use the immersed moving boundary method to deal with the boundary between the fluid and the solid. For the grid at the interface between the fluid and the solid, calculate the proportion of the solid area in the grid ε(a x ,a y,t) and speed U s (a x ,a y ,t);

[0030] If grid(a x ,a y ) only intersects with the jth segment of the long strip solid, then its velocity U s (a x ,a y ,t) is:

[0031]

[0032] Among them, l P (a x ,a y ,j,t) is the grid at the current time t (a x ,a y ) points to the center of mass of the j-th long solid segment;

[0033] If grid(a x ,a y ) intersects with multiple long strip solid segments, then the velocity U at the intersection with each long strip solid segment is calculated separately. s (a x ,a y ,t), and then take the average value;

[0034] Step 5.2: Except for the boundary grid of the fluid calculation domain, perform collision operations on other grids. For each grid that performs collision operations, calculate the collision terms Ω in the nine directions. i (a x ,a y ,t); the collision term of the boundary grid of the fluid calculation domain is Ω i (a x ,a y ,t)=(0,0);

[0035]

[0036] in:

[0037]

[0038]

[0039] Among them, -i is the opposite direction of the i-th direction,

[0040] Step 5.3: Calculate the total force F exerted by the fluid on the long solid bar f (t) and the torque T on each long solid segment f(j,t);

[0041]

[0042] Step 5.4: Perform migration operations on each grid in the fluid computational domain. For each grid, calculate the updated distribution vector f in the nine directions. i (a x ,a y ,t+dt);

[0043] f i (a x ,a y ,t+dt)=f i (a x ,a y ,t)+Ω i (a x ,a y ,t).

[0044] Furthermore, the DEM cycle is performed in step 6, specifically including the following steps:

[0045] Step 6.1: Initialize t 2 = 0, set the iteration step size Δt of the DEM sub-cycle; ω 0 (j,t+dt)=ω dtime (j,t), θ 0 (j,t+dt)=θ dtime (j,t);

[0046] Step 6.2: Traverse each long strip solid segment, obtain the contact status between each long strip solid segment, and form a contact list; for each combination of two long strip solid segments in the contact list, determine whether they are in a bonding state; if two long strip solid segments are determined to be out of the bonding state in a certain calculation, then in subsequent calculations, the two will not be re-bonded with any long strip solid segment;

[0047] If we judge the jth 1 Segment long solid segment and the jth 2 If the segments of the long strip solid are in a bonding state, calculate the ultimate normal stress at the contact surface of the two segments and determine whether the ultimate normal stress is less than the ultimate tensile strength and ultimate compressive strength of the long strip solid material; if so, determine the jth 1 Segment long solid segment and the jth 2 If there is no break between the segments of the long strip solid, execute step 6.3; otherwise, determine the jth 1 Segment long solid segment and the jth 2There is a break between the long strip solid segments, and there is no force vector and torque between the two long strip solid segments, that is, the force vector and torque are both zero;

[0048] No. 1 Segment long solid segment and the jth 2 The deformation at the bonding surface of the long strip solid segment is expressed as The corresponding normal deformation component is The tangential deformation component is The limiting normal stresses are and At the upper and lower ends of the bonding surface, the normal stress The calculation method is:

[0049]

[0050] If the j 1 Segment long solid segment and the jth 2 There is overlapping area between the segments of a long solid Then determine the jth 1 Segment long solid segment and the jth 2 If there is a collision between the long strip solid segments, execute step 6.4;

[0051] Step 6.3: For the jth 1 Segment long solid segment and the jth 2 The force vectors on the two long strip solid segments are equal in magnitude and opposite in direction; 1 The force vector of a long solid segment The calculation method is:

[0052]

[0053] in, The jth 1 The normal elastic force, normal viscous force, tangential elastic force and tangential viscous force on the long strip solid segment;

[0054]

[0055] No. 1 The torque on a long solid segment The calculation method is:

[0056]

[0057] in, For the jth 1 Segment long solid segment and the jth2 The centroid of the bonding surface of the long strip solid segment points to the jth 1 Segment: The vector of the centroid of a segment of a long solid;

[0058] Step 6.4: For the jth collision 1 Segment long solid segment and the jth 2 The force vectors on the two long strip solid segments are equal in magnitude and opposite in direction; 1 The force vector of a long solid segment The calculation method is:

[0059]

[0060] in, The overlapping area The angle between the line connecting the center of mass of and the origin of the plane rectangular coordinate system of the fluid calculation domain and the x-axis; and The jth 1 Normal contact force and tangential contact force on the long strip solid segment;

[0061]

[0062] in, For the jth 1 Segment long solid segment and the jth 2 Normal vector of the contact surface of the segmented solid strip; For the jth 1 Segment long solid segment and the jth 2 The tangent vector of the contact surface of the segmented solid strip; For the jth 1 Segment long strip solid segment relative to the jth 2 The velocity vector of a segment of a long solid;

[0063]

[0064] in, The overlapping area The center of mass points to the jth 1 Segment: The vector of the centroid of a segment of a long solid; The overlapping area The center of mass points to the jth 2 Segment: The vector of the centroid of a segment of a long solid;

[0065]

[0066] No. 1 The torque on a long solid segment The calculation method is:

[0067]

[0068] Step 6.5: Update the speed of each solid strip segment Angular velocity Displacement in the x-axis direction Displacement in the y-axis direction and azimuth

[0069]

[0070] Step 6.6: If t 2 <dtime, then let t 2 =t 2 +Δt, return to step 6.2; otherwise, output t 2 = the speed of each long solid segment at dtime The average angular velocity of rotation is ω dtime (j, t+dt), displacement in the x-axis direction Displacement in the y-axis direction and azimuth angle θ dtime (j,t+dt).

[0071] A computer device / equipment / system comprises a memory, a processor and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above-mentioned research method for realizing solid deformable crushing based on LBM-DEM coupling algorithm.

[0072] A computer-readable storage medium stores a computer program / instruction, which, when executed by a processor, implements the steps of the research method for realizing solid deformable crushing based on the LBM-DEM coupling algorithm.

[0073] A computer program product includes a computer program / instruction, which, when executed by a processor, implements the steps of the research method for realizing solid deformable crushing based on the LBM-DEM coupling algorithm.

[0074] The beneficial effects of the present invention are:

[0075] The present invention uses the LBM method to divide the fluid calculation domain into equally spaced grids, each grid is given 9 directions and all grids in the calculation domain are subjected to collision and migration operations to simulate the flow of the flow field; the DEM method is used to divide the solid structure into equal-sized rectangular units, and the collision or bonding state of each rectangular unit is identified, and the interaction force is calculated to simulate the deformation and destruction of the solid structure; the immersed moving boundary method is used to process the boundary between the fluid and the solid, transfer the force between the fluid and the solid structure, identify all grid division node types in the calculation domain, and calculate the volume fraction ratio of the boundary type nodes, and then calculate and solve the force of the fluid on the solid structure. The present invention can correctly describe the mechanical properties of fragile solid materials in fluids, and can be used to solve related engineering problems such as the field of fluid-solid coupling where solids are deformed or destroyed. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 Schematic diagram of the calculation direction of the grid (using the D2Q9 discrete model).

[0077] Figure 2 Schematic diagram of the bounding box search algorithm.

[0078] Figure 3 Schematic diagram of the IMB method model and node density distribution probability calculation.

[0079] Figure 4 Identify schematic diagrams for node search algorithms.

[0080] Figure 5 Schematic diagram of LBM collision step.

[0081] Figure 6 Schematic diagram of LBM migration steps.

[0082] Figure 7 Schematic diagram of a simple Hooke's law contact model.

[0083] Figure 8 Schematic diagram of the parallel bonding contact model.

[0084] Fig. 9 It is the overall flow chart of the present invention.

[0085] Fig.10 Schematic diagram of the geometric structure of the calculation model in an embodiment of the present invention.

[0086] Fig.11 It is a diagram of calculation results in an embodiment of the present invention.

[0087] Fig.12 FIG. 4 is a graph showing the displacement of the top of the thin plate over time in an embodiment of the present invention. DETAILED DESCRIPTION

[0088] The present invention is further described below in conjunction with the accompanying drawings.

[0089] In the current LBM-DEM method field, the existing technology can only simulate the movement and collision of a simple single or multiple rigid solids in a fluid, but the solid cannot be deformed or destroyed, and it is impossible to characterize the properties and movement trends of the breakable material in the fluid. The technical solution of the present invention solves the current problem, realizes that the breakable solid in the fluid can be deformed and destroyed into multiple individuals, and correctly describes the mechanical properties of the breakable solid material in the fluid.

[0090] The present invention is mainly used to solve the fluid-solid coupling problem of deformable and broken solids in computational fluid dynamics, and can be used in many fluid-solid coupling fields such as marine engineering, environmental engineering, and transportation engineering. For example, in the polar ocean, ice is a solid material that is easily broken, and there are many large and small icebergs-ice ridges in the water. When a ship or underwater vehicle is sailing, the ice ridges may be split due to the influence of temperature, flow field disturbances, etc., and ice blocks may fall off, posing a hidden danger to the safe navigation of ships or underwater vehicles. The present invention can predict whether an ice ridge may break during navigation and the movement trajectory of the ice blocks after breaking, so that we can avoid it in advance to avoid collisions and safety accidents. The invention provides a certain reference value.

[0091] like Fig. 9 As shown in FIG. 1 , a research method for realizing solid deformable crushing based on the LBM-DEM coupling algorithm uses the mesoscopic physical background of the lattice Boltzmann (LBM) method, the particle characteristics and the characteristics of the discrete element (DEM) nonlinear continuous medium to deal with the fluid-solid coupling problem, including the following steps:

[0092] Step 1: For the problem of long solids moving, colliding and deformable breakage in a liquid-filled tube, determine the total predicted time T and obtain the flow velocity u of the fluid 0 , density ρ 0 , viscosity visco and relaxation factor τ; construct a two-dimensional fluid calculation domain, determine the initial arrangement of the long strip solid in the fluid calculation domain, the upper and lower boundaries of the fluid calculation domain are the upper and lower walls of the liquid-filled tube, and the fluid calculation domain has a left boundary and a right boundary;

[0093] Take the intersection of the left boundary and the lower boundary of the fluid calculation domain as the origin, take the lower boundary as the x-axis, the positive direction of the x-axis is the flow direction of the fluid, take the y-axis of the left boundary, and establish a plane rectangular coordinate system; divide the fluid calculation domain into grids, each grid is a square grid with a side length of dx, and the index of each grid is (a x ,a y ), represents the ath x Row a y A grid of columns, a x =1,2,...,Nx , a y =1,2,...,N y , the total area of ​​the fluid calculation domain is N x ·dx·N y ·dx;

[0094] Step 2: If Figure 1 As shown, for each grid, take its center point as the 0th direction, the vertical line from the center point to the right side of the grid as the 1st direction, the vertical line from the center point to the left side of the grid as the 2nd direction, the vertical line from the center point to the top of the grid as the 3rd direction, the vertical line from the center point to the bottom of the grid as the 4th direction, the vertical line from the center point to the upper right vertex of the grid as the 5th direction, the vertical line from the center point to the upper left vertex of the grid as the 6th direction, the vertical line from the center point to the lower left vertex of the grid as the 7th direction, and the vertical line from the center point to the lower right vertex of the grid as the 8th direction, construct the vector e i for:

[0095]

[0096] Step 3: Get the length L of the long solid 1 、Width L 2 , density ρ s , elastic modulus E, tangential modulus G, Poisson's ratio poisson, normal stiffness coefficient k of the material n , tangential stiffness coefficient k t , normal viscosity coefficient k nv , tangential viscosity coefficient k tv , damping coefficient η, friction coefficient μ, ultimate tensile strength, ultimate compressive strength; divide the long strip solid evenly into N s The mass of each segment of the long solid strip is The moment of inertia is Calculate the time step dt of the LBM main cycle and the time step dtime of the DEM secondary cycle 0 ;

[0097]

[0098]

[0099] in, dtime 0 Correction is made so that the corrected dtime satisfies dt=N 2 dtime, N 2 Is a positive integer, and dtime≤dtime 0 ;

[0100] Step 4: Execute the LBM main loop, initialize t = 0, initialize the long strip solid segments to be in a bonded state, and the speed, rotation angular velocity, displacement and azimuth of each long strip solid segment are all 0, that is, ω dtime (j,0)=0, θ 0 (j,0) = 0; Initialize the distribution vector f in 9 directions of each grid in the fluid calculation domain i (a x ,a y ,0) and the equilibrium distribution vector

[0101] in,

[0102] Step 5: Use the immersed moving boundary method to process the boundary between the fluid and the solid. For the grid at the interface between the fluid and the solid, calculate the proportion of the solid area in the grid ε(a x ,a y ,t) and speed U s (a x ,a y ,t);

[0103] like Figure 3 As shown, the present invention uses the immersed moving boundary method (IMB) to process the boundary between fluid and solid. The IMB method characterizes the no-slip condition of the interface between fluid and solid by the density distribution probability of the nodes covered by solid particles.

[0104] The main idea is to divide the entire computational domain into four types of nodes. First, nodes that are close to the fluid-solid boundary but belong to the solid are determined as solid boundary nodes (yellow nodes in the figure); nodes that belong to the solid except these yellow nodes are determined as solid internal nodes (red nodes in the figure); nodes that are close to the fluid-solid boundary and belong to the fluid on the outside are determined as fluid boundary nodes (cyan nodes in the figure); all other fluid nodes are determined as normal fluid nodes (blue nodes in the figure). Then the area proportion of the grid solid is calculated, and then the weighted score is calculated.

[0105] like Figure 4 As shown in the figure, for any polygon, the specific method of dividing the four types of nodes is:

[0106] 1. First, obtain the minimum circumscribed rectangular boundary of the target polygon;

[0107] 2. After numbering the nodes, perform interval division and node identification;

[0108] ① Calculate the coordinates of the intersections of the row and the polygon, and sort these intersections in ascending order according to the x-axis coordinate values, for example Figure 4 In the intersection of a row and a rectangular unit, the increasing order is b→e;

[0109] ② Combine the intersection points into intervals (be intervals) according to the arrangement order. The maximum and minimum grid points in the interval are identified as solid boundary nodes (grid points c and d); the remaining nodes in the interval are identified as solid internal nodes; the grid points outside the interval adjacent to the solid boundary nodes are identified as fluid boundary nodes (grid points a and f);

[0110] ③Put the identified four types of nodes into four sets respectively;

[0111] 3. Calculate the coordinates of the intersections of all numbered columns and polygons, repeat step 2, and obtain a set of 4 types of nodes;

[0112] 4. Perform deduplication processing on each set of the four types of nodes, and the result is the final result.

[0113] If grid(a x ,a y ) only intersects with the jth segment of the long strip solid, then its velocity U s (a x ,a y ,t) is:

[0114]

[0115] Among them, l P (a x ,a y ,j,t) is the grid at the current time t (a x ,a y ) points to the center of mass of the j-th long solid segment;

[0116] If grid(a x ,a y ) intersects with multiple long strip solid segments, then the velocity U at the intersection with each long strip solid segment is calculated separately. s (a x ,a y ,t), and then take the average value;

[0117] Step 6: Except for the boundary grid of the fluid calculation domain, perform collision operations on other grids. For each grid that performs collision operations, calculate the collision terms Ω in the nine directions. i (a x ,a y ,t); the collision term of the boundary grid of the fluid calculation domain is Ω i (a x ,a y ,t)=(0,0);

[0118] LBM describes the flow of fluid by solving the collision (evolution) and migration process of fluid particles in the computational domain. The collision process mainly evolves in the internal area of ​​the computational domain and does not involve the boundary; migration requires traversing all nodes in the cyclic computational domain, taking the current node as the center and migrating to the corresponding nodes in the surrounding 9 directions, replacing the distribution function of the corresponding direction of the corresponding node, such as Figure 5 and Figure 6 shown.

[0119]

[0120] in,

[0121]

[0122]

[0123] Among them, -i is the opposite direction of the i-th direction,

[0124] Step 7: Calculate the total force F exerted by the fluid on the long solid f (t) and the torque T on each long solid segment f (j,t);

[0125]

[0126] Step 8: Perform migration operations on each grid in the fluid calculation domain. For each grid, calculate the updated distribution vector f in the nine directions. i (a x ,a y ,t+dt);

[0127] f i (a x ,a y ,t+dt)=f i (a x ,a y ,t)+Ω i (a x ,a y ,t)

[0128] Step 9: Execute the DEM loop to update the velocity of each long solid segment Angular velocity ω dtime (j, t+dt), displacement in the x-axis direction Displacement in the y-axis direction and azimuth angle θ dtime (j,t+dt);

[0129] The DEM sub-cycle mainly characterizes the deformation and damage of solid materials. The main contact models are divided into two categories. One is the contact model that follows the simple Hooke's law between discrete blocks after the plate is broken; the other is the parallel bonding contact model between discrete blocks of the deformed plate. The force between discrete blocks after the plate is broken is related to the overlapping area between the blocks. After the above-mentioned bounding box search algorithm, whether the two blocks collide depends on whether they intersect, that is, whether there is an overlapping area, such as Figure 7 shown.

[0130] When the plate does not break, the parallel bond contact model is used between the blocks. Since there is no break, the positions of adjacent blocks are determined, and there is no need to perform adjacent search. It is only necessary to calculate the structural force and deformation of the given blocks. For example, Figure 8 As shown, there will be compressive stress and tensile stress at the contact surface between block i and block j at the same time, with compressive stress (compressive deformation) generated at the upper end and tensile stress (tensile deformation) generated at the lower end.

[0131] Step 9.1: Initialize t 2 = 0, set the iteration step size Δt of the DEM sub-cycle; ω 0 (j,t+dt)=ω dtime (j,t), θ 0 (j,t+dt)=θ dtime (j,t);

[0132] Step 9.2: Traverse each long strip solid segment, obtain the contact conditions between each long strip solid segment, and form a contact list; Figure 2 As shown, the contact determination uses a bounding box search algorithm to determine whether there is contact by determining the maximum and minimum coordinate relationship of each solid circumscribed rectangular bounding box.

[0133] For every combination of two long strip solid segments in the contact list, determine whether they are in a bonding state; if two long strip solid segments are determined to be out of the bonding state in a certain calculation, they will not be re-bonded with any long strip solid segment in subsequent calculations;

[0134] If we judge the jth 1 Segment long solid segment and the jth 2 If the segments of the long strip solid are in a bonding state, calculate the ultimate normal stress at the contact surface of the two segments and determine whether the ultimate normal stress is less than the ultimate tensile strength and ultimate compressive strength of the long strip solid material; if so, determine the jth 1 Segment long solid segment and the jth 2 If there is no break between the segments of the long strip solid, execute step 9.3; otherwise, determine the jth 1Segment long solid segment and the jth 2 There is a break between the long strip solid segments, and there is no force vector and torque between the two long strip solid segments, that is, the force vector and torque are both zero;

[0135] No. 1 Segment long solid segment and the jth 2 The deformation at the bonding surface of the long strip solid segment is expressed as The corresponding normal deformation component is The tangential deformation component is The limiting normal stresses are and At the upper and lower ends of the bonding surface, the normal stress The calculation method is:

[0136]

[0137] If the j 1 Segment long solid segment and the jth 2 There is overlapping area between the segments of a long solid Then determine the jth 1 Segment long solid segment and the jth 2 If there is a collision between the long strip solid segments, execute step 9.4;

[0138] Step 9.3: For the jth 1 Segment long solid segment and the jth 2 The force vectors on the two long strip solid segments are equal in magnitude and opposite in direction; 1 The force vector of a long solid segment The calculation method is:

[0139]

[0140] in, The jth 1 The normal elastic force, normal viscous force, tangential elastic force and tangential viscous force on the long strip solid segment;

[0141]

[0142] No. 1 The torque on a long solid segment The calculation method is:

[0143]

[0144] in, For the jth1 Segment long solid segment and the jth 2 The centroid of the bonding surface of the long strip solid segment points to the jth 1 Segment: The vector of the centroid of a segment of a long solid;

[0145] Step 9.4: For the jth collision 1 Segment long solid segment and the jth 2 The force vectors on the two long strip solid segments are equal in magnitude and opposite in direction; 1 The force vector of a long solid segment The calculation method is:

[0146]

[0147] in, The overlapping area The angle between the line connecting the center of mass of and the origin of the plane rectangular coordinate system of the fluid calculation domain and the x-axis; and The jth 1 Normal contact force and tangential contact force on the long strip solid segment;

[0148]

[0149] in, For the jth 1 Segment long solid segment and the jth 2 Normal vector of the contact surface of the segmented solid strip; For the jth 1 Segment long solid segment and the jth 2 The tangent vector of the contact surface of the segmented solid strip; For the jth 1 Segment long strip solid segment relative to the jth 2 The velocity vector of a segment of a long solid;

[0150]

[0151] in, The overlapping area The center of mass points to the jth 1 Segment: The vector of the centroid of a segment of a long solid; The overlapping area The center of mass points to the jth 2 Segment: The vector of the centroid of a segment of a long solid;

[0152]

[0153] No. 1The torque on a long solid segment The calculation method is:

[0154]

[0155] Step 9.5: Update the speed of each solid strip segment Angular velocity Displacement in the x-axis direction Displacement in the y-axis direction and azimuth

[0156]

[0157] Step 9.6: If t 2 <dtime, then let t 2 =t 2 +Δt, return to step 9.2; otherwise, output t 2 = the speed of each long solid segment at dtime The average angular velocity of rotation is ω dtime (j, t+dt), displacement in the x-axis direction Displacement in the y-axis direction and azimuth angle θ dtime (j,t+dt);

[0158] Step 10: If t<T, set t=t+dt and return to step 5; otherwise, output the calculation result when t=T to complete the prediction of the movement, collision and deformable breakage of the long solid in the liquid-filled tube.

[0159] Embodiment 1:

[0160] In this embodiment, if Fig.10 As shown in the figure, a thin plate perpendicular to the inflow is placed in a fluid pipe. The length of the fluid pipe is L = 4 cm and the height is H = 1 cm. The thickness of the thin plate is 0.04 cm, the length is b = 0.8 cm, and the aspect ratio is 1 / 20. The thin plate is placed 1.0 cm away from the velocity inlet of the pipe and the lower end is fixed. The fluid density is set to 1.0 g / cm 3 , the fluid viscosity is set to 0.1 g / (cm·s), and the sheet density is 7.8 g / cm 3 , Young's modulus is 10 5 g / (cm·s 2 ), Poisson's ratio is 0.3, and the ultimate compressive strength is set to 6000g / (cm·s 2 ), the ultimate tensile strength is set to -8000g / (cm·s 2 ), and the gravity and damping forces of the thin plate are neglected in this simulation.

[0161] Computational domain boundary conditions: Apply parabolic velocity inlet boundary condition V on the left side of the fluid domain x =1.5(-y 2 +2y)cm / s,V y = 0; the right boundary is set to the pressure outlet boundary condition p = 0; the upper and lower boundaries are no-slip fixed boundary conditions, that is, V x =0, V y =0.

[0162] The calculation results are as follows Fig.11 As shown, it can be seen that the top point of the thin plate is driven by the inlet fluid to move in the early stage, and then fluctuates slightly back and forth. At about 0.3s, the bottom end of the thin plate reaches the ultimate strength and breaks, and swims to the outlet with the fluid. Fig.12 The x-direction displacement and y-direction displacement of the top point of the thin plate drift over time.

[0163] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A research method for realizing solid deformable crushing based on LBM-DEM coupling algorithm, characterized in that: The following steps are involved: Step 1: For the problem of movement, collision and deformable breakage of a long strip solid in a liquid-filled tube, determine the total prediction time, obtain the parameter information of the fluid and the long strip solid, and construct a two-dimensional fluid calculation domain. The upper and lower boundaries of the fluid calculation domain are the upper and lower walls of the liquid-filled tube, and have a left boundary and a right boundary; determine the initial layout of the long strip solid in the fluid calculation domain; Step 2: Divide the long strip solid into multiple segments evenly in the length direction; divide the fluid calculation domain evenly into grids, each grid is the same square grid, and determine the calculation direction of each grid; Step 3: Calculate the time step of the LBM main cycle and the duration of the DEM secondary cycle based on the parameter information of the fluid and the long strip solid; Step 4: Initialize the segments of the long strip solid to be in a bonding state, and the speed, rotation angular velocity, displacement and azimuth of each long strip solid segment are all zero; initialize the distribution vector of each grid in each calculation direction in the fluid calculation domain; Step 5: Execute the LBM main loop and use the immersed moving boundary method to process the boundaries of the fluid and the solid. Except for the boundary grids of the fluid calculation domain, perform collision operations on other grids. For each grid that performs collision operations, calculate the collision term based on the distribution vectors of each calculation direction of the grid. The collision term of the boundary grid of the fluid calculation domain is a zero vector. Calculate the total force of the fluid on the long strip solid and the torque on each long strip solid segment based on the collision term. Finally, perform migration operations on each grid in the fluid calculation domain and update the distribution vectors of each calculation direction of each grid. Step 6: Execute the DEM sub-cycle until the duration of the DEM sub-cycle is reached; In each iteration of the DEM cycle, each long strip solid segment is traversed to obtain the contact conditions between each long strip solid segment and form a contact list; for each combination of two long strip solid segments in the contact list, it is determined whether they are in a bonding state; if two long strip solid segments are determined to be out of the bonding state in a certain calculation, then in subsequent calculations, the two will not be re-bonded with any long strip solid segment; If it is judged that the two long strip solid segments are in a bonding state, the limit normal stress at the contact surface of the two segments is calculated, and it is judged whether the limit normal stress is less than the limit tensile strength and the limit compressive strength of the long strip solid material; if it is less than, it is judged that no fracture occurs between the two long strip solid segments, and the force vector and torque received by each long strip solid segment are calculated; otherwise, it is judged that a fracture occurs between the two long strip solid segments, and there is no force vector and torque between the two; If there is an overlapping area between the two long strip solid segments, it is determined that a collision occurs between the two, and the force vector and torque received by each long strip solid segment are calculated; For each long strip solid segment, according to the total force of the fluid on the long strip solid, the torque of the fluid on the long strip solid segment, and the force vector and torque of the long strip solid segment from other long strip solid segments, update the speed, angular velocity of rotation, displacement in the x-axis direction, displacement in the y-axis direction and azimuth of the long strip solid segment; Step 7: Repeat steps 5 to 6 according to the time step of the LBM main loop until the total prediction time is reached, output the calculation results, and complete the prediction of the movement, collision and deformable breakage of the long solid in the liquid-filled tube.

2. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 1 is characterized in that: In the step 1, the total prediction time T is determined, and the parameter information of the fluid and the long strip solid is obtained, including the flow velocity u0, density ρ0, viscosity visco and relaxation factor τ of the fluid, the length L1, width L2, density ρ s , the elastic modulus E, tangential modulus G, Poisson's ratio poisson, and normal stiffness coefficient k of the long strip solid material n , tangential stiffness coefficient k t , normal viscosity coefficient k nv , tangential viscosity coefficient k tv , damping coefficient η, friction coefficient μ, ultimate tensile strength, ultimate compressive strength.

3. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 2 is characterized in that: In step 2, the long strip solid is evenly divided into N s The mass of each segment of the long solid strip is The moment of inertia is Take the intersection of the left boundary and the lower boundary of the fluid calculation domain as the origin, take the lower boundary as the x-axis, the positive direction of the x-axis is the flow direction of the fluid, take the y-axis of the left boundary, and establish a plane rectangular coordinate system; divide the fluid calculation domain into grids, each grid is a square grid with a side length of dx, and the index of each grid is (a x ,a y ), represents the ath x Row a y A grid of columns, a x =1,2,...,N x , a y =1,2,...,N y , N x With N y are the number of rows and columns of the grid respectively, and the total area of ​​the fluid calculation domain is N x ·dx·N y ·dx; For each grid, take its center point as the 0th direction, the vertical line from the center point to the right side of the grid as the 1st direction, the vertical line from the center point to the left side of the grid as the 2nd direction, the vertical line from the center point to the top of the grid as the 3rd direction, the vertical line from the center point to the bottom of the grid as the 4th direction, the vertical line from the center point to the upper right vertex of the grid as the 5th direction, the vertical line from the center point to the upper left vertex of the grid as the 6th direction, the vertical line from the center point to the lower left vertex of the grid as the 7th direction, and the vertical line from the center point to the lower right vertex of the grid as the 8th direction, and construct the vector e i for:

4. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 3 is characterized in that: The method for calculating the time step dt of the LBM main cycle and the time step dtime of the DEM secondary cycle in step 3 is specifically: in, Correct dtime0 so that the corrected dtime satisfies dt=N2·dtime, N2 is a positive integer, and dtime≤dtime0.

5. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 4 is characterized in that: In step 4, t=0 is initialized, and the segments of the long strip solid are initialized to be in a bonding state, and the speed, rotation angular velocity, displacement and azimuth of each long strip solid segment are all 0, that is, ω dtime (j,0)=0, θ 0 (j,0) = 0; Initialize the distribution vector f in 9 directions of each grid in the fluid calculation domain i (a x ,a y ,0) and the equilibrium distribution vector f i eq (a x ,a y ,0); in, 6. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 5 is characterized in that: The LBM main cycle is executed in step 5, which specifically includes the following steps: Step 5.1: Use the immersed moving boundary method to deal with the boundary between the fluid and the solid. For the grid at the interface between the fluid and the solid, calculate the proportion of the solid area in the grid ε(a x ,a y ,t) and speed U s (a x ,a y ,t); If grid(a x ,a y ) only intersects with the jth segment of the long strip solid, then its velocity U s (a x ,a y ,t) is: Among them, l P (a x ,a y ,j,t) is the grid at the current time t (a x ,a y ) points to the center of mass of the j-th long solid segment; If grid(a x ,a y ) intersects with multiple long strip solid segments, then the velocity U at the intersection with each long strip solid segment is calculated separately. s (a x ,a y ,t), and then take the average value; Step 5.2: Except for the boundary grid of the fluid calculation domain, perform collision operations on other grids. For each grid that performs collision operations, calculate the collision terms Ω in the nine directions. i (a x ,a y ,t); the collision term of the boundary grid of the fluid calculation domain is Ω i (a x ,a y ,t)=(0,0); in: Among them, -i is the opposite direction of the i-th direction, Step 5.3: Calculate the total force F exerted by the fluid on the long solid bar f (t) and the torque T on each long solid segment f (j,t); Step 5.4: Perform migration operations on each grid in the fluid computational domain. For each grid, calculate the updated distribution vector f in the nine directions. i (a x ,a y ,t+dt); f i (a x ,a y ,t+dt)=f i (a x ,a y ,t)+Ω i (a x ,a y ,t)。 7. The method for realizing solid deformable crushing based on LBM-DEM coupling algorithm according to claim 6 is characterized in that: The DEM cycle is performed in step 6, which specifically includes the following steps: Step 6.1: Initialize t2 = 0 and set the iteration step length Δt of the DEM sub-cycle; ω 0 (j,t+dt)=ω dtime (j,t), θ 0 (j,t+dt)=θ dtime (j,t); Step 6.2: Traverse each long strip solid segment, obtain the contact status between each long strip solid segment, and form a contact list; for each combination of two long strip solid segments in the contact list, determine whether they are in a bonding state; if two long strip solid segments are determined to be out of the bonding state in a certain calculation, then in subsequent calculations, the two will not be re-bonded with any long strip solid segment; If it is determined that the j1th long strip solid segment and the j2th long strip solid segment are in a bonding state, the ultimate normal stress at the contact surface of the two segments is calculated, and it is determined whether the ultimate normal stress is less than the ultimate tensile strength and the ultimate compressive strength of the long strip solid material; if it is less than, it is determined that no fracture occurs between the j1th long strip solid segment and the j2th long strip solid segment, and step 6.3 is executed; otherwise, it is determined that a fracture occurs between the j1th long strip solid segment and the j2th long strip solid segment, and there is no force vector and torque between the two long strip solid segments, that is, the force vector and the torque are both zero; The deformation at the bonding surface between the j1th long strip solid segment and the j2th long strip solid segment is expressed as The corresponding normal deformation component is The tangential deformation component is The limiting normal stresses are and At the upper and lower ends of the bonding surface, the normal stress The calculation method is: If there is an overlapping area between the j1th long strip solid segment and the j2th long strip solid segment It is determined that the j1th long strip solid segment collides with the j2th long strip solid segment, and step 6.4 is executed; Step 6.3: For the j1st long solid segment and the j2nd long solid segment that are not broken and in a bonded state, the force vectors on the two long solid segments are equal in magnitude and opposite in direction; among them, the force vector on the j1st long solid segment is The calculation method is: in, They are respectively the normal elastic force, normal viscous force, tangential elastic force, and tangential viscous force on the j1th long strip solid segment; The torque on the j1th long solid segment The calculation method is: in, is the vector from the center of mass of the bonding surface between the j1th long strip solid segment and the j2th long strip solid segment to the center of mass of the j1th long strip solid segment; Step 6.4: For the j1th long strip solid segment and the j2th long strip solid segment that collide, the force vectors on the two long strip solid segments are equal in magnitude and opposite in direction; among them, the force vector on the j1th long strip solid segment is The calculation method is: in, The overlapping area The angle between the line connecting the center of mass of and the origin of the plane rectangular coordinate system of the fluid calculation domain and the x-axis; and are the normal contact force and tangential contact force on the j1th long strip solid segment respectively; in, is the normal vector of the contact surface between the j1th long strip solid segment and the j2th long strip solid segment; is the tangent vector of the contact surface between the j1th long strip solid segment and the j2th long strip solid segment; is the velocity vector of the j1th long strip solid segment relative to the j2th long strip solid segment; in, The overlapping area The center of mass of the j1th long solid segment points to the center of mass of the j1th long solid segment; The overlapping area The center of mass of the j2-th long solid segment points to the vector of the center of mass of the j2-th long solid segment; The torque on the j1th long solid segment The calculation method is: Step 6.5: Update the speed of each solid strip segment Angular velocity Displacement in the x-axis direction Displacement in the y-axis direction and azimuth Step 6.6: If t2 < dtime, set t2 = t2 + Δt and return to step 6.2; otherwise, output the speed of each long solid segment when t2 = dtime The average angular velocity of rotation is ω dtime (j, t+dt), displacement in the x-axis direction Displacement in the y-axis direction and azimuth angle θ dtime (j,t+dt).

8. A computer device / equipment / system comprising a memory, a processor and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program / instruction stored thereon, characterized in that: When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

10. A computer program product comprising a computer program / instructions, characterized in that: When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 7 are implemented.

Images