A method for predicting particle motion trajectories in a DLD microfluidic sorting device

Abstract

The present invention discloses a method for predicting the motion trajectory of particles in a DLD microfluidic sorting device; the method is as follows: 1. Construct a DLD microfluidic sorting device model in multi-physics field finite element simulation software, and obtain the velocity field and pressure field of the fluid domain in the model. 2. Perform boundary extraction and reconstruction on the DLD microfluidic sorting device model constructed in step 1, and iterate the particle motion trajectory. The present invention proposes a method for predicting the Dc value of an actual device using 2D data, which expands the particle from a 2D circle to a 3D sphere in a 2D simulation and performs spherical segmentation, taking into account the influence of height on parameters such as flow rate, and optimizing the influence of particle velocity and fluid pressure on particles and iterative operations, so as to achieve more accurate and efficient prediction of the motion trajectory of particles in a deterministic lateral displacement device, and obtain algorithm prediction results that are almost consistent with the actual test results.

Description

Technical Field

[0001] The present invention belongs to the technical field of DLD microfluidic sorting device design, and in particular relates to a method for predicting particle motion trajectories in a DLD microfluidic sorting device. Background Art

[0002] Circulating tumor cells (CTCs) are tumor cells that enter the peripheral blood. Tumor cell metastasis is a major cause of high cancer mortality. The sizes of CTCs and blood cells are micrometer-sized, ranging from 5 to 35 μm, and CTCs are larger than most blood cells. Numerous microfluidic sorting technologies have been proposed based on these size differences, including surface acoustic wave sorting, inertial spiral sorting, and deterministic lateral displacement (DLD) sorting. Deterministic lateral displacement (DLD) sorting devices are passive microfluidic sorting technologies with advantages such as compact devices, simple testing systems, and high precision. Since its proposal by Lotien Richard Huang [Continuous particle separation through deterministic lateral displacement. Science 304, 987-990], numerous application studies have been conducted. Various micropillar structures, including triangular, circular, parallelogram, and L-shaped, have been designed for different samples, as well as their respective geometric arrangements. DLD sorting devices exhibit a critical diameter, or Dc value. When the particle size is larger than the Dc value, it exhibits a bump trajectory, causing the particle to deviate from the main flow line. When the particle size is smaller than the Dc value, it exhibits a zigzag trajectory, aligning with the main flow line. Ultimately, these different particle trajectories enable size-based sorting.

[0003] However, the Dc value of a DLD device depends on the geometry of its internal micropillar array. Typically, achieving the desired structural parameters requires multiple iterations of fabrication, which undoubtedly increases labor and time costs. To mitigate this shortcoming, researchers have developed empirical formulas based on simplified models. However, existing empirical formulas cannot cover all geometries. Wang Junchao et al. proposed the MOPSA algorithm [MOPSA: A microfluidics-optimized particle simulation algorithm. Biomicrofluidics 11,14]. However, this algorithm contains a hyperparameter factor β, which is related to the fluid properties, micropillar shape, and array parameters, requiring correction based on test results from actual fabricated devices. Therefore, the MOPSA algorithm cannot effectively predict the Dc value of newly designed DLD devices. The particle trajectory function in the multiphysics simulation software COMSOL treats particles as point masses, disregarding their size, and therefore cannot directly simulate the Dc value of DLD devices. Therefore, an algorithm that can quickly and easily calculate the Dc value of DLD devices would be of great value in reducing trial-and-error costs and accelerating device development. Summary of the Invention

[0004] The purpose of the present invention is to propose a method for predicting particle motion trajectories in a DLD microfluidic sorting device, so as to predict the Dc value (particle sorting critical diameter) of the DLD sorting device.

[0005] In a first aspect, the present invention provides a method for predicting particle motion trajectories in a DLD microfluidic sorting device, comprising the following steps:

[0006] Step 1: Construct a DLD microfluidic sorting device model in a multi-physics field finite element simulation software and obtain the velocity field and pressure field of the fluid domain in the model.

[0007] Step 2: Extract and reconstruct the boundaries of the DLD microfluidic sorting device model constructed in step 1, and iterate the particle motion trajectory.

[0008] 2-1. Obtain all boundary segment information in the model constructed in step 1.

[0009] 2-2. Initialize particle parameters.

[0010] 2-3. Divide the particle surface into multiple unit surfaces. Collect the coordinates of the center point of each unit surface.

[0011] 2-4. Initialize the particle position and velocity.

[0012] 2-5. Divide the DLD microfluidic sorting device model into multiple unit areas arranged in a matrix.

[0013] 2-6. Iterative calculation of particle trajectories.

[0014] a) Calculate the x-direction flow velocity v at the center point of each unit surface while taking into account the height information x and the y-direction velocity v y , and based on the flow velocity v at the center of each unit surface x and flow rate v y , obtain the x-direction flow velocity u at the particle location f and the y-direction velocity v f .

[0015] b) According to the flow velocity (uf, vf) and the current velocity of the particle (u current ,v current ) calculate the drag force (F D_x ,F D_y ).

[0016] c) Calculate the pressure value P at the center point of each unit surface p , and based on the pressure value P at the center point of each unit surface p Get the pressure resultant force (F p_x ,F p_y ).

[0017] d) According to the current speed (u current ,v current ), drag force (F D_x ,F D_y ), pressure force (F p_x ,F p_y ) predicts the displacement of the particle from the current moment to the next moment (Δs x ,Δs y ) and obtain the position at the next moment (P next_x ,P next_y ) and speed (u next ,v next ).

[0018] e) According to whether the particle position coincides with the microcolumn, the particle position (P next_x ,P next_y ) and particle velocity (u next ,v next ) for correction.

[0019] f) Repeat the iteration until the particle leaves the fluid domain and draw the particle trajectory.

[0020] Preferably, step 2 is repeated using particles of different diameters; and the predicted particle sorting critical diameter of the DLD microfluidic sorting device is obtained based on the trajectories of the particles of different diameters.

[0021] Preferably, the process of constructing the model in step 1 is as follows: first, the shape of the micropillars is drawn, and then an array is formed to obtain a micropillar array corresponding to the predicted parameters of the DLD microfluidic sorting device; then, a rectangle is drawn that partially or completely covers the micropillar array; the difference between the rectangle and the micropillar array is used as the fluid domain; and finally, the fluid domain is meshed.

[0022] Preferably, in step 2, only multiple unit areas around the location of the particle are extracted in each calculation.

[0023] Step 2 is completed in MATLAB software.

[0024] Preferably, in step 2-1, based on the contour's characteristic of being connected end to end, the originally randomly distributed boundary line segments are reordered by coordinate value. If the endpoints of two closed curves coincide, a gradient prediction method is used to determine the next coordinate point. The rearranged coordinate points can be connected to form a closed pattern.

[0025] Preferably, the particle parameters in step 2-2 include particle diameter, particle density, particle volume and particle mass.

[0026] Preferably, in step 2-3, the particles are divided into multiple pieces along the height direction, and then divided into multiple pieces along the circumference using a vertical surface.

[0027] As a preference, in step b) of steps 2-6, the x-direction velocity v of the midpoints of all unit surfaces is taken. x The mean of the flow velocity u at the particle location f ; The y-direction flow velocity v of all unit surface midpoints y The mean of the flow velocity v at the particle location f The flow velocity v at the midpoint of the unit surface x and flow rate v y The acquisition process is the same as that of , which is: obtain the coordinates of the projection point P of the center of each unit surface on the particle on the XOY plane of the particle local coordinate system. Determine the triangle mesh position of the projection point P in the fluid domain; calculate the velocity value V of the projection point P p =(V B -V A )·μ+(V C -V A )·η+V A ; Among them, V A , V B , V C are the velocity values ​​at the three vertices A, B, and C of the triangle mesh. μ and η are coefficients determined by the position of the projection point P in the triangle mesh. p , calculate the velocity value of each unit surface midpoint on the particle Wherein, H is the channel thickness of the DLD microfluidic sorting device; P z The Z-axis coordinate value of the center point of the unit surface in the particle's local coordinate system.

[0028] Preferably, in step b) of steps 2-6, the resultant pressure force (F p_x ,F p_y ) is:

[0029]

[0030] Among them, m represents the layer number along the height direction; n represents the segmentation number along the circumferential direction; M represents the number of equal divisions along the height direction; N represents the number of equal divisions along the circumferential direction. m(n+N / 2) 、p mn The pressure value of the center point of the nth unit surface of the mth layer from top to bottom;

[0031] As a preferred method, the process of obtaining the pressure value of the midpoint of the unit surface is as follows: obtain the coordinates of the projection point P of the center of each unit surface on the particle on the XOY plane of the particle local coordinate system. Determine the triangle mesh position of the projection point P in the fluid domain; calculate the pressure value P corresponding to the projection point P p =(P B -P A )·μ+(P C -P A )·η+P A ; where P A , P B , P C are the pressure values ​​at the three vertices A, B, and C of the triangle mesh where the projection point P is located. μ and η are coefficients determined by the position of the projection point P in the triangle mesh. The pressure value at the projection point P is used as the pressure value at the center of the corresponding unit surface.

[0032] In a second aspect, the present invention provides a DLD microfluidic sorting device design method, which comprises the following steps:

[0033] Step 1: Set the target particle sorting critical diameter and the initial parameters of the DLD microfluidic sorting device.

[0034] Step 2: Using particles corresponding to the target particle sorting critical diameter, the predicted particle sorting critical diameter of the currently set DLD microfluidic sorting device is obtained according to the aforementioned particle motion trajectory prediction method within the DLD microfluidic sorting device.

[0035] Step 3: According to the difference between the predicted particle sorting critical diameter and the target particle sorting critical diameter, adjust the DLD microfluidic sorting device and repeat step 2 until the predicted particle sorting critical diameter is within a preset range centered on the target particle sorting critical diameter.

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

[0037] 1. This paper proposes a method for predicting the Dc value of an actual device using 2D data. In a 2D simulation, particles are expanded from 2D circles into 3D spheres and then segmented. This method considers the effect of height on parameters such as flow velocity, the effects of particle velocity and fluid pressure on particle motion, and iterative calculation optimization. This method achieves a more accurate and efficient prediction of particle motion trajectories in a deterministic lateral displacement device (DLD), obtaining algorithmic prediction results that are nearly consistent with actual test results.

[0038] 2. Compared with the empirical formula and the original MOPSA algorithm, the present invention improves the prediction accuracy of particle trajectories in DLD devices, especially the prediction accuracy for new and unknown DLD devices, which helps to realize the digital design of DLD devices, thereby quickly obtaining DLD devices with target Dc values.

[0039] 3. The algorithm proposed in the present invention can not only predict the Dc value of the DLD device with high accuracy, but also significantly reduce the amount of calculation and accelerate the device finalization. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 Flowchart of Example 1 of the present invention

[0041] Figure 2 Flow field diagram of the velocity field u in the DLD fluid domain obtained in step 1 of Example 1 of the present invention;

[0042] Figure 3 Flow field diagram of the velocity field v in the DLD fluid domain obtained in step 1 of Example 1 of the present invention;

[0043] Figure 4 The pressure field diagram in the DLD fluid domain obtained in step 1 of Example 1 of the present invention;

[0044] Figure 5 Schematic diagram of expanding the 2D circle of the particle model into a 3D sphere and performing finite element segmentation in step 2 of Example 1 of the present invention;

[0045] Figure 6 1 is a geometric comparison diagram of the DLD device model and the actual device in Example 1 of the present invention;

[0046] Figure 7 The particle trajectory diagram outputted by Example 1 of the present invention;

[0047] Figure 8 This is a mask diagram for preparing a physical device in Example 1 of the present invention;

[0048] Figure 9 This is a photo of the inlet of the physical device prepared in Example 1 of the present invention, where particles with diameters of 10 μm and 11 μm are input;

[0049] Figure 10 This is the fluorescence image at the outlet of the physical device prepared in Example 1 of the present invention when inputting particles with diameters of 10 μm and 11 μm;

[0050] Figure 11 A schematic diagram of structural parameters of a DLD device designed for Example 2 of the present invention;

[0051] Figure 12 This is a simulation diagram showing a zigzag trajectory of the flow trajectory of particles with a diameter of 7.5 μm when the vertical gap is 30 μm in Example 2 of the present invention;

[0052] Figure 13 This is a simulation diagram of a bump trajectory of a particle with a diameter of 7.5 μm when the vertical gap is 20 μm in Example 2 of the present invention;

[0053] Figure 14 This is a simulation diagram of a zigzag trajectory of a particle with a diameter of 7.5 μm when the vertical gap is 25 μm in Example 2 of the present invention;

[0054] Figure 15 This is a simulation diagram showing that when the vertical gap is 22.5 μm in Example 2 of the present invention, the flow trajectory of particles with a diameter of 7.5 μm presents a bump trajectory, and the flow trajectory of particles with a diameter of 7.4 μm presents a zigzag trajectory. DETAILED DESCRIPTION

[0055] The present invention will be further described below with reference to the accompanying drawings.

[0056] In the following embodiments, the preferred software implementation platforms are COMSOL Multiphysics and MATLAB, and the two software are linked using the COMSOL Multiphysics with MATLAB program interface. The description process is divided into the COMSOL software process and the MATLAB software process, as shown in the flow chart. Figure 1 shown.

[0057] Example 1

[0058] A method for predicting particle motion trajectories in a DLD microfluidic sorting device comprises the following steps:

[0059] Step 1: Complete the model construction of the DLD microfluidic sorting device in the multi-physics field finite element simulation software to obtain the velocity field and pressure field of the flow domain. The specific process is as follows:

[0060] 1-1. In this embodiment, COMSOL Multiphysics software is used as the multi-physics finite element simulation software; in the modeling wizard of the multi-physics finite element simulation software, "2D model" is selected, "laminar flow" is selected as the physical field, and "steady-state study" is selected as the study.

[0061] 1-2. Create a model in the geometry node of the multi-physics finite element simulation software. The process is as follows:

[0062] a) Draw the shape of a single micropillar. The cross-section of the micropillar is an inverted triangle with a corner radius R of 2.5 μm and a height of 25 μm.

[0063] b) Obtain a tilted row of micropillars using a linear array in the array. The number of micropillars in the array is 20, and the x offset of the linear array is 60 μm and the y offset is 3.125 μm.

[0064] c) Using the matrix array in the array for the linear array pattern obtained in step b), a matrix array with multiple rows is obtained. The number of arrays in the X direction is 1, the number of arrays in the Y direction is 13, and the y offset is 50 μm.

[0065] d) Draw a rectangle covering part of the matrix array with a length of 1.53 mm and a width of 0.55 mm.

[0066] e) Use the Subtraction command to subtract the micropillar array from the rectangle drawn in step d) to obtain the shape of the fluid domain. The fluid domain is the area within the rectangle not occupied by the micropillar array.

[0067] 1-3. The material of the fluid domain uses water from the material library of the multi-physics finite element simulation software.

[0068] In the Laminar Flow module, set the boundaries to No Slip, add inlet and outlet boundaries, and set the leftmost boundary of the fluid domain as the inlet with an average flow velocity of 10 mm / s. Set the leftmost boundary of the fluid domain as the outlet with a pressure of 0 Pa.

[0069] 1-5. Select physical field control, free triangle mesh, and cell size as fine.

[0070] 1-6. Click Calculate in the study. After the calculation is completed, the velocity field u in the X direction of the fluid domain can be obtained, such as Figure 2 As shown, the velocity field v in the Y direction is Figure 3 As shown, the pressure field p is Figure 4 shown.

[0071] 1-7. Save the simulation project.

[0072] Step 2: Complete the boundary extraction and reconstruction of the DLD microfluidic sorting device model constructed in step 1, particle parameter initialization, variable initialization, iterative calculation of particle motion trajectory, and particle trajectory display in MATLAB software. The specific process is as follows:

[0073] 2-1. Use the mphopen function in MATLAB to allow MATLAB to access the model project file created in step 1.

[0074] 2-2. Use the mpheval function in MATLAB to obtain the velocity field u in the X direction, the velocity field v in the Y direction, and the pressure field p in the model. At this time, the corresponding values ​​are only available at the nodes of the triangular mesh.

[0075] 2-3. Use the mphmeshstats function in MATLAB to obtain information about all boundary segments in the model. Because contours are always connected end to end, the originally randomly distributed boundary segments are reordered by coordinate value and duplicate coordinate points are removed. If the endpoints of two closed curves coincide (share a coordinate point), the next coordinate point is determined using a gradient prediction method. The rearranged coordinate points can be connected to form a closed pattern. All coordinate points of a closed boundary are placed in a separate cell. Storing closed boundaries separately facilitates program access and improves calculation speed.

[0076] 2-4. Initialize particle parameters, including determining particle size information, particle density, particle volume (3D), and particle mass (3D). In this embodiment, the particle density is ρ = 1.05 g / cm 3 ; Particle volume calculation formula The particle mass calculation formula is m = ρ·V. D is the particle diameter. In this embodiment, the particle diameter D is an undetermined parameter, and the motion trajectory of particles with different diameters is ultimately obtained. The goal of this embodiment is to determine the critical particle diameter that can be separated by the model constructed in step 1.

[0077] 2-5. Particle surface division: parallel to the XOY plane, cut the particle into 8 slices of equal thickness; perpendicular to the XOY plane, divide each slice into 16 equal parts to obtain 8×16 unit surfaces; Figure 5 As shown. Take the particle center as the origin of the coordinate system and record the midpoint coordinates of each unit surface. The center coordinates of each unit surface (x m,n ,y m,n ,z m,n ) is as follows:

[0078]

[0079] Where m represents the number of slices from top to bottom (parallel to the XOY plane), ranging from 1 to 8. n represents the number of slices perpendicular to the XOY plane, ranging from 1 to 16. R in the formula represents the particle radius, R = D / 2.

[0080] 2-6. Position and velocity initialization: In the DLD microfluidic sorting device model, set the current midpoint position of the particle (P currrent_x ,P current_y ) is (20μm, 20μm); set the current speed of the particle (u current ,v current ) has an initial value of (0.01,0) m / s

[0081] 2-7. Data preprocessing: The entire DLD microfluidic sorting device model is divided into 31×11 unit regions according to the coordinate space. The corresponding velocity and pressure field values ​​are stored in each small region. During each calculation, only the data and boundaries within the unit region near the particle location (in this example, nine 3×3 unit regions around the particle location) are extracted for calculation. This reduces the calculation of invalid data and speeds up the calculation, which is particularly significant when the data volume is very large.

[0082] 2-8. Iterative Calculation of Particle Trajectories

[0083] a) Set the condition for stopping the calculation. The iteration is completed when the particle leaves the fluid domain.

[0084] b) Obtain the velocity at the particle position by taking the projection coordinates of the center of each unit surface on the particle on the XOY plane of the particle's local coordinate system, i.e., the x- and y-coordinates in steps 2-5. Use the vector method Determine the triangle mesh where the projection point P is located, where A, B, and C are the coordinate points of the triangle mesh. The necessary and sufficient conditions for the projection point P to be inside the triangle mesh ABC are μ∈[0.1], η∈[0.1], and μ+η≤1.

[0085] After that, calculate the x-direction flow velocity v of all unit surface midpoints on the particle x and the y-direction velocity v y ; Calculate the flow velocity v in the x direction x and the y-direction velocity v y The process is the same as the following: Calculate the velocity value V of the projection point P p =(V B -V A )·μ+(V C -V A )·η+V A , where V A , V B , V CThe velocity values ​​at the three vertices A, B, and C of the triangle mesh are respectively. Further considering the parabolic distribution of velocity in the Z direction; according to the velocity value V of the projection point P p , calculate the velocity value of each unit surface midpoint on the particle Wherein, H is the predicted channel thickness of the DLD microfluidic sorting device, which is set to 25 μm in this embodiment; P z The Z-axis coordinate value of the center point of the unit surface in the particle's local coordinate system.

[0086] From the spherical crown formula, we know that the area of ​​each layer from top to bottom is The area of ​​each surface is Therefore, take the x-direction velocity v of all the midpoints of the unit surface x The mean of the flow velocity u at the particle location f ; The y-direction flow velocity v of all unit surface midpoints y The mean of the flow velocity v at the particle location f .

[0087] c) Calculate the Stokes drag force (F D_x ,F D_y )=6πηR·Δv=6πηR·[(u f ,v f )-(u current ,v current )]. Where η is the dynamic viscosity of the fluid, R is the radius of the particle, and Δv is the velocity difference between the fluid velocity and the particle velocity. (F D_x ,F D_y ) represent the drag force in the X direction and the drag force in the Y direction, respectively.

[0088] d) Obtain the pressure on the particle by taking the projection point P of the center of each unit surface on the particle on the XOY plane of the particle's local coordinate system; calculate the pressure value P at the projection point P p =(P B -P A )·μ+(P C -P A )·η+P A ; where P A , P B , P C They are the pressure values ​​on the three vertices A, B, and C of the triangle mesh where the projection point P is located. Pressure value P p Only the pressure component parallel to the XOY plane is considered, and the pressure component in the Z direction is ignored.

[0089] Combined with the relationship of pressure = pressure × area, we can get the X-direction pressure force F on the particle: p_x , Y-direction pressure force F p_yas follows:

[0090]

[0091] Among them, m represents the slice number from top to bottom (parallel to the XOY plane), ranging from 1 to 8; n represents the slice number perpendicular to the XOY plane, ranging from 1 to 8; there are 16 slices perpendicular to the XOY plane, but (p m(n+8) -p mn ) item considers the difference of symmetrical points, so the data index value of the last 8 points can be obtained by adding 8 to the first 8 points. mn Indicates the pressure value on the nth unit surface in the mth slice perpendicular to the XOY plane. m represents the surface area of ​​the mth slice, which is a fraction of the particle's outer surface area.

[0092] e) Use Newton's second law and the acceleration formula to predict the displacement of the particle along the X and Y directions within Δt (Δs x ,Δs y )as follows:

[0093]

[0094] Where, mp is the mass of the particle; (u current ,v current ) represents the X and Y components of the particle's current velocity. (Δs x ,Δs y ) represents the displacement within the time Δt.

[0095] f) Add the current position of the particle to the displacement within Δt time and predict the position of the particle after Δt time (P next_x ,P next_y )as follows:

[0096] (P next_x ,P next_y )=(P currrent_x ,P current_y )+(Δs x ,Δs y ),

[0097] Among them, (P currrent_x ,P current_y ) represents the X and Y coordinates of the current particle center; (P next_x ,P next_y ) represents the X- and Y-coordinates of the particle center position after the predicted Δt time.

[0098] g) Use Newton's second law and the velocity formula to predict the velocity of the particle in the X and Y directions after Δt time (u next ,vnext )as follows:

[0099]

[0100] Where mp is the mass of the particle. (u next ,v next ) represents the velocity of the particle in the X direction and the velocity in the Y direction after Δt time.

[0101] h) Avoid the particle from coinciding with the wall. next_x ,P next_y ), the particle's outline is overlapped with the outline of the micropillar wall near the particle to determine whether the particle overlaps with the wall. If overlap occurs, the maximum distance the particle embeds into the wall is calculated. This maximum distance is then decomposed into dx in the X direction and dy in the Y direction. dx is positive if the particle moves in the positive X direction, and negative if it moves in the negative X direction. dy is positive if the particle moves in the positive Y direction, and negative if it moves in the negative Y direction.

[0102] Corrected particle position (P next_x ,P next_y ) is as follows:

[0103] (P next_x ,P next_y )=(P next_x,0 ,P next_y,0 )+(dx,dy)

[0104] Among them, (P next_x,0 ,P next_y,0 ) is the position of the particle before correction.

[0105] Correct the predicted position of the particle after Δt time, and halve the predicted velocity after Δt time to obtain the corrected particle velocity (u next ,v next ):

[0106]

[0107] Among them, (u next,0 ,v next,0 ) is the particle velocity before correction.

[0108] If no overlap occurs, no correction is required.

[0109] i) Variable preservation and update. Assign the original particle center coordinates after Δt to the variable representing the current particle center coordinates (P currrent_x ,P current_y )=(P next_x ,P next_y) and save the data of the current particle center coordinate position. Assign the particle velocity after the original Δt moment to the variable representing the current particle velocity (u current ,v current )=(u next ,v next )

[0110] j) Determine whether the particle area is still within the fluid domain. If the particle area is still within the fluid domain, repeat steps b) to i) to start a new round of iteration; otherwise, the iteration is completed.

[0111] 2-9. Draw the outline of the fluid domain and the trajectory of the particles.

[0112] Step 3: Simulate the trajectories of particles of different diameters in the DLD microfluidic sorting device to obtain the Dc value (particle sorting critical diameter) of the DLD device: Set multiple monotonically changing particle diameters D and execute step 2 separately to obtain the particle trajectories corresponding to different particle diameters. The particle diameter D at which the particle trajectory changes significantly is used as the Dc value of the DLD device. The specific process is as follows:

[0113] If the particle trajectory exhibits a zigzag trajectory in the DLD microfluidic sorting device model, increase the particle diameter D and repeat the trajectory algorithm in step 2 until the trajectory exhibits an altered or bump trajectory. The particle diameter D obtained in step 2 is considered the Dc value of the DLD device being tested.

[0114] If the particle trajectory shows an altered trajectory or a bump trajectory, reduce the D value and repeatedly run the trajectory algorithm in step 2 until the trajectory shows a Zigzag trajectory. The particle diameter D obtained from the last execution of the trajectory algorithm in step 2 is considered to be the Dc value of the DLD device under test.

[0115] This embodiment uses a DLD device using triangular micropillars as the prediction target, such as Figure 6 As shown in the figure, the fillet radius of the triangular microcolumn is 2.5μm and the height of the triangle is 25μm. The offset between adjacent rows is 3.125μm, the horizontal spacing between rows is 60μm, the column spacing is 50μm, and the pipe height is 25μm. The Dc value of the DLD device predicted by the algorithm provided in this embodiment is 10.6μm. The predicted trajectory is shown in Figure 7 The actual size is as shown. Figure 6 The actual object and design dimensions are similar. Figure 8 The mask diagram of the device is shown. The actual device tests the 10μm and 11μm phosphor particles. The 10μm and 11μm phosphor particles are not separated in the entrance area of ​​the device. Figure 9As shown. Two separate fluorescent bands appear in the exit area, as shown Figure 10 As shown, it is shown that the particle sorting of 10μm and 11μm is achieved.

[0116] Example 2

[0117] A DLD device design method comprises the following steps:

[0118] Step 1: Set the initial parameters of the DLD device and the target Dc value.

[0119] The design goal is to obtain a cylindrical microcolumn DLD device with a target Dc value of 7.5μm. The initial setting is 25μm for the channel height, 50μm for the diameter of the cylinder, 20μm for the horizontal microcolumn gap, 30μm for the vertical microcolumn gap, and 3.2μm for the row offset. Figure 11 As shown. The particle size D is set to 7.5 μm, and the trajectory is obtained using the algorithm in Implementation 1. Figure 12 As shown, a zigzag trajectory is observed. This indicates that the current geometric parameters do not meet the design requirements and need to be optimized.

[0120] Step 2: Optimize the DLD device process based on the target Dc value.

[0121] In this embodiment, the Dc value of the DLD device is adjusted by changing the gap distance in the vertical direction.

[0122] Change the vertical gap to 20 μm, start from modeling again, calculate and draw the particle trajectory with a diameter of 7.5, as shown in Figure 13 As shown, a bump trajectory is observed. Substituting the particle diameter D = 7.4 μm into the calculation, the trajectory also shows a bump trajectory. This indicates that the vertical clearance value obtained for Dc = 7.5 μm is between 20 μm and 30 μm.

[0123] Change the vertical gap to 25 μm, start from modeling again, calculate and draw the trajectory of particles with a diameter of 7.5, as shown in Figure 14 As shown, a zigzag trajectory is exhibited, indicating that the gap value obtained for Dc=7.5μm is between 20μm and 25μm.

[0124] Change the vertical gap to 22.5μm and start again from modeling. Calculate and draw the particle trajectory with a diameter of 7.5. The particle trajectory is a bump trajectory. Then, substitute the particle diameter D = 7.4μm into the calculation. The trajectory is a zigzag trajectory, as shown in the following example: Figure 15Thus, the geometric parameters of the DLD device with Dc=7.5μm were determined: the channel height is 25μm, the diameter of the cylinder is 50μm, the horizontal micro-pillar gap is 20μm, the vertical micro-pillar gap is 22.5μm, and the row offset is 3.2μm.

[0125] The above is only a preferred embodiment of the present invention. It should be noted that those skilled in the art can make several improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered as the scope of protection of the present invention.

Claims

1. A method for predicting particle motion trajectories in a DLD microfluidic sorting device, characterized by: The following steps are involved: Step 1: Build a DLD microfluidic sorting device model in multi-physics finite element simulation software and obtain the velocity field and pressure field of the fluid domain in the model; Step 2: Extract and reconstruct the boundaries of the DLD microfluidic sorting device model constructed in step 1, and iterate the particle motion trajectory; Step 2-1. Obtain all boundary segment information in the model constructed in step 1; Step 2-2. Initialize particle parameters; Step 2-3. Divide the particle surface into multiple unit surfaces; collect the coordinates of the center point of each unit surface; Step 2-4. Initialize the particle position and velocity; Step 2-5. Divide the DLD microfluidic sorting device model into multiple unit areas arranged in a matrix; Step 2-6. Iterative calculation of particle trajectories; a) Calculate the x-direction flow velocity v at the center point of each unit surface while taking into account the height information x and the y-direction velocity v y , and based on the flow velocity v at the center of each unit surface x and flow rate v y , obtain the x-direction flow velocity u at the particle location f and the y-direction velocity v f ; b) According to the flow velocity (uf, vf) and the current velocity of the particle (u current ,v current ) calculate the drag force (F D_x ,F D_y ); c) Calculate the pressure value P at the center point of each unit surface p , and based on the pressure value P at the center point of each unit surface p Get the pressure resultant force (F p_x ,F p_y ); d) According to the current speed (u current ,v current ), drag force (F D_x ,F D_y ), pressure force (F p_x ,F p_y ) predicts the displacement of the particle from the current moment to the next moment (Δs x ,Δs y ) and obtain the position at the next moment (P next_x ,P next_y ) and speed (u next ,v next ); e) According to whether the particle position coincides with the microcolumn, the particle position (P next_x ,P next_y ) and particle velocity (u next ,v next ) f) Repeat the iteration until the particle leaves the fluid domain and obtain the particle trajectory.

2. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 1, wherein: Repeat step 2 using particles of different diameters; obtain the predicted particle sorting critical diameter of the DLD microfluidic sorting device based on the trajectories of the particles of different diameters.

3. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 1, wherein: The process of building the model in step 1 is as follows: first draw the shape of the micropillars, then array them to obtain a micropillar array corresponding to the predicted DLD microfluidic sorting device parameters; then draw a rectangle that partially or completely covers the micropillar array to obtain the fluid domain; finally, mesh the fluid domain.

4. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 1, wherein: In step 2-1, based on the characteristic of the contour being connected end to end, the originally disordered boundary line segments are reordered according to their coordinate values. If the endpoints of two closed curves coincide, the gradient prediction method is used to determine the next coordinate point. The rearranged coordinate points can be connected to form a closed pattern.

5. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 1, wherein: The particle parameters described in step 2-2 include particle diameter, particle density, particle volume, and particle mass.

6. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 1, characterized in that: In step 2, only multiple unit areas around the particle location are extracted for each calculation; in steps 2-3, the particle is divided into multiple equal pieces along the height direction, and then divided into multiple equal parts along the circumference using the vertical plane.

7. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 6, characterized in that: In step b) of step 2-6, take the x-direction velocity v of all the midpoints of the unit surface x The mean of the flow velocity u at the particle location f ; The y-direction flow velocity v of all unit surface midpoints y The mean of the flow velocity v at the particle location f ; Flow velocity v at the midpoint of the unit surface x and flow rate v y The acquisition process is the same as that of , which is: obtain the coordinates of the projection point P of the center of each unit surface on the particle on the XOY plane of the particle local coordinate system; determine the triangle mesh position of the projection point P in the fluid domain; Calculate the velocity value V of the projection point P p =(V B -V A )·μ+(V C -V A )·η+V A ; Among them, V A , V B , V C They are the velocity values ​​on the three vertices A, B, and C of the triangular mesh respectively; μ and η are coefficients determined by the position of the projection point P in the triangular mesh; according to the velocity value V of the projection point P p , calculate the velocity value of each unit surface midpoint on the particle Wherein, H is the channel thickness of the DLD microfluidic sorting device; P z The Z-axis coordinate value of the center point of the unit surface in the particle's local coordinate system.

8. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 6, characterized in that: In step b) of step 2-6, the pressure resultant force (F p_x ,F p_y ) is: Where m represents the layer number along the height direction; n represents the segmentation number along the circumferential direction; M represents the number of equal divisions along the height direction; N represents the number of equal divisions along the circumferential direction; p m(n+N / 2) 、p mn The pressure value of the center point of the nth unit surface of the mth layer from top to bottom.

9. The method for predicting particle motion trajectories in a DLD microfluidic sorting device according to claim 6, wherein: The process of obtaining the pressure value of the midpoint of the unit surface is as follows: obtain the coordinates of the projection point P of the center of each unit surface on the particle on the XOY plane of the particle local coordinate system; determine the triangle mesh position of the projection point P in the fluid domain; calculate the pressure value P corresponding to the projection point P p =(P B -P A )·μ+(P C -P A )·η+P A ; where P A , P B , P C are the pressure values ​​on the three vertices A, B, and C of the triangle mesh where the projection point P is located; μ and η are coefficients determined by the position of the projection point P in the triangle mesh; The pressure value of the projection point P is used as the pressure value of the corresponding unit surface center.

10. A DLD microfluidic sorting device design method, characterized by: The following steps are involved: Step 1: Set the target particle sorting critical diameter and the initial parameters of the DLD microfluidic sorting device; Step 2: Using particles corresponding to the target particle sorting critical diameter, the predicted particle sorting critical diameter of the currently set DLD microfluidic sorting device is obtained according to the particle motion trajectory prediction method in the DLD microfluidic sorting device as described in claim 2; Step 3: According to the difference between the predicted particle sorting critical diameter and the target particle sorting critical diameter, adjust the DLD microfluidic sorting device and repeat step 2 until the predicted particle sorting critical diameter is within a preset range centered on the target particle sorting critical diameter.

Images