Numerical method for engulfment behavior of immiscible droplets in confined shear flow

Abstract

The present application provides a kind of numerical method of limited shear flow in immiscible droplet engulfment behavior, it is related to the technical field of immiscible droplet engulfment behavior, the simulation task based on phase field model and Navier-Stokes equation is established in the present application, the calculation domain is dispersed by finite element method and the physical field is initialized, in each time step, the updated velocity field and phase field data are obtained by coupling solution, then the interface curvature change rate and the vorticity intensity of flow field are comprehensively calculated, the global index representing the dynamic intensity of system is generated, the index is compared with the preset threshold value, the adjustment strategy of next time step is dynamically selected: when the index is higher than the threshold value, smaller step and higher grid resolution are used, otherwise, larger step and lower resolution are used, finally, the solver configuration is automatically reconstructed according to the selected parameters, the simulation is continued under the new adjustment strategy, forming the adaptive numerical simulation cycle, until the simulation is completed and the complete field data of droplet engulfment process is output.

Description

Technical Field

[0001] This invention relates to the field of the field of droplet swallowing behavior of immiscible solutions, specifically a numerical method for the droplet swallowing behavior of immiscible solutions in confined shear flow. Background Technology

[0002] In fields such as microfluidics, materials synthesis, and chemical separation, the engulfment behavior of immiscible droplets in confined shear flow is a key interfacial phenomenon affecting the mass transfer efficiency and product structure uniformity of the system. Accurate simulation of such processes is of great value for optimizing equipment design and process parameters. Especially in applications involving the preparation of multiple emulsions and the controlled release of drug carriers, where precise control of interface morphology is required, numerical methods have become an indispensable tool for studying microscopic interface dynamics. Although traditional computational fluid dynamics methods can describe basic flows, they still have significant limitations in the dynamic response to droplet interface topological changes and flow field coupling effects.

[0003] Existing technologies mainly advance numerical solutions through fixed grid systems and constant time steps, or by pre-setting the grid refinement region based on experience at specific times. These methods rely on prior knowledge to set discrete parameters and adopt a uniform numerical resolution in different stages such as droplet approach, film drainage, and final fusion. Some improved schemes attempt to introduce interface tracking or static adaptive grids, but their parameter adjustments are often out of sync with the real-time physical field evolution and cannot autonomously optimize the allocation of computing resources based on the dynamic characteristics of the system.

[0004] The shortcomings of existing technologies lie in their lack of real-time perception and response capabilities to abrupt changes in physical field characteristics. When the curvature of the droplet interface changes drastically or the vorticity of the flow field increases instantaneously, the fixed parameter system cannot synchronously improve the local calculation accuracy, resulting in blurred interface capture or loss of key dynamic features. On the other hand, maintaining high-resolution calculations during the smooth flow field phase leads to the unnecessary consumption of a large amount of computing resources. This contradiction is essentially due to the fact that existing methods decouple the numerical discretization strategy from the physical field evolution process, making it impossible to establish an autonomous decision-making mechanism based on real-time dynamic characteristics, thus making it difficult to achieve a dynamic balance between computational efficiency and numerical accuracy.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a numerical method for the engulfment behavior of immiscible droplets in confined shear flow, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A numerical method for the co-occurrence behavior of immiscible droplets in confined shear flow, comprising the following steps:

[0009] Step 1: Establish a simulation task of confined shear flow in the computer system. Use the phase-field equation coupled with the Navier-Stokes equation as the governing equation. Use the finite element method to spatially discretize the computational domain to generate a computational grid composed of several grid elements. Assign initial states to all grid elements according to the initial physical conditions to generate the initial phase field and velocity field. Set the initial time step and start the time step iteration process.

[0010] Step 2: Perform coupled numerical solution based on the current time step and current grid resolution to obtain the updated phase field numerical solution and velocity field numerical solution for each grid cell. Calculate the interface curvature change rate and flow field vorticity intensity of each grid cell, and generate a global index characterizing the dynamic intensity of the entire system through weighted aggregation.

[0011] Step 3: Compare the global metrics with the preset thresholds. Based on the comparison results, dynamically select a set of adjustment strategies for the next time step. The adjustment strategy is to use a smaller time step and a higher grid resolution level when the global metrics are higher than the preset thresholds, and a larger time step and a lower grid resolution level when they are lower.

[0012] Step 4: Automatically reconstruct the discrete configuration of the numerical solver according to the selected adjustment strategy, that is, divide the computational grid by applying the time step and grid resolution level. The iterative process returns to Step 2, and performs coupled numerical solution for the next time step according to the selected adjustment strategy. This loop continues until the simulation task is completed, and finally outputs serialized field data for analyzing the co-occurrence behavior of immiscible droplets.

[0013] Furthermore, establishing the simulation task specifically includes:

[0014] Establishing the phase-field equation coupled with the Navier-Stokes equation as the governing equation means using the Cahn-Hilliard equation to describe the evolution process of the immiscible fluid interface and the transient Navier-Stokes equation to describe the motion process of the background fluid. The two are coupled through the interfacial tension term. Spatial discretization adopts the finite element method, which means dividing the continuous computational domain into multiple non-overlapping grid elements to form an initial computational grid for numerical calculation.

[0015] The state includes a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The logic for generating the phase field and velocity field is as follows: based on the initial physical conditions, each grid cell is assigned a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The scalar values ​​representing the fluid phase identity in all grid cells are aggregated to generate the phase field, and the vector values ​​representing the fluid velocity in all grid cells are aggregated to generate the velocity field. The initial physical conditions include the initial flow field velocity distribution, the initial droplet position and size, and the physical property parameters of each fluid. The physical property parameters include density, viscosity, interfacial tension coefficient, mobility, and interfacial thickness parameters.

[0016] Furthermore, for coupled numerical solutions, a step-by-step solution algorithm is used to solve the phase field control equations and the fluid dynamics equations sequentially within each time step, thereby obtaining the updated phase field numerical solution and velocity field numerical solution defined on all grid cells.

[0017] The specific process of generating global metrics includes:

[0018] Based on the numerical solution of the phase field at the current time step, the curvature field is obtained by calculating its second spatial derivative, and the curvature rate of change field is obtained by subtracting it from the curvature field at the previous time step; based on the numerical solution of the velocity field at the current time step, the vorticity field is obtained by calculating its curl.

[0019] Calculate the L2 norm of the curvature change rate field over the entire computational domain and weight it using a preset first weighting coefficient. Calculate the L2 norm of the vorticity field over the entire computational domain and weight it using a preset second weighting coefficient. Sum the two weighted results to obtain the global index.

[0020] The first weighting coefficient and the second weighting coefficient are determined based on the sensitivity analysis of each physical attribute parameter, and satisfy the condition that the sum of the first weighting coefficient and the second weighting coefficient is 1.

[0021] Furthermore, the preset threshold includes a preset high threshold and a preset low threshold, and the adjustment strategy includes a high-precision parameter strategy, a high-efficiency parameter strategy, and a maintenance parameter strategy. The logic for dynamically selecting the adjustment strategy is as follows:

[0022] The global index is compared and analyzed with the preset high and low thresholds: when the value of the global index is greater than the high threshold, a high-precision parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset refinement coefficient less than 1, and the grid resolution level is adjusted by adding a positive integer number of level adjustment amounts to the current level.

[0023] When the value of the global metric is less than the low threshold, a high-efficiency parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset coarsening coefficient greater than 1, and the grid resolution level is adjusted by reducing the current level by a positive integer number of level adjustments.

[0024] When the value of the global metric is between the high and low thresholds, a parameter maintenance strategy is selected for the next time step to maintain the current time step size and grid resolution level unchanged.

[0025] The preset refinement coefficient, preset coarsening coefficient, and level adjustment amount are preset based on numerical stability requirements during simulation initialization.

[0026] Furthermore, for the discrete configuration of the automatically reconfigured numerical solver, the following logical steps are executed according to the selected adjustment strategy:

[0027] 1) By calling the parameter setting interface provided by the numerical solver, the time step size of the next time step is set to a new time step size determined based on the selected adjustment strategy;

[0028] 2) Trigger the mesh adaptation function of the numerical solver to identify regions of interest based on the numerical solutions of the velocity field and phase field at the current time step. Specifically: based on the numerical solution of the phase field, calculate its spatial gradient magnitude distribution and mark all mesh elements with gradient magnitudes greater than a preset interface gradient threshold as high-dynamic regions of interest; based on the numerical solution of the velocity field, calculate its vorticity magnitude distribution and mark all mesh elements with vorticity magnitudes less than a preset vorticity intensity threshold as low-dynamic regions of interest.

[0029] 3) Based on the selected adjustment strategy, perform corresponding mesh topology transformations on the region of interest: when a high-precision parameter strategy is selected, refine the mesh cells in the high-dynamic region of interest by performing a mesh resolution level refinement operation; when a high-efficiency parameter strategy is selected, coarsen the mesh cells in the low-dynamic region of interest by performing a mesh resolution level coarsening operation; when a maintenance parameter strategy is selected, keep all regions of the current computational mesh unchanged.

[0030] Furthermore, the termination condition of the simulation task is determined by the following logic: the simulation task is considered complete when any of the following conditions are met:

[0031] 1) The accumulated physical simulation time since the start of the simulation reaches the preset simulation duration limit;

[0032] 2) The following three sub-conditions must be met simultaneously:

[0033] 2.1) Droplet agglomeration process completed: The droplet agglomeration process is considered complete when only one connected droplet phase region exists within the computational domain, determined by performing connectivity analysis on the current phase field numerical solution. 2.2) Droplet interface motion convergence: The droplet interface motion convergence is considered complete when the maximum change in the droplet centroid position over N consecutive time steps is less than a preset position convergence threshold. 2.3) System dynamic convergence: The system dynamic convergence is considered complete when the difference between the maximum and minimum values ​​of the global index over N consecutive time steps is less than a preset index convergence threshold.

[0034] Where N is the preset number of time steps greater than 1.

[0035] Furthermore, the iterative process returns to step 2, whose operational logic is as follows:

[0036] After completing the discrete configuration reconstruction in step 4, the system automatically advances the computation state to the next time step, and based on the reconstructed time step and grid resolution, executes the coupled numerical solution process in step 2 again, thus forming a closed loop of time stepping.

[0037] When the simulation task is completed, the system decodes and encodes the velocity field and phase field values ​​corresponding to all time steps into binary data streams according to the time step sequence, and writes them to the storage device to form a single data file.

[0038] Compared with the prior art, the beneficial effects of the present invention are:

[0039] This invention effectively solves the core problem of the disconnect between computational resource allocation and physical field evolution in existing technologies by establishing a real-time feedback mechanism between the dynamic intensity index of the physical field and the numerical discrete parameters. Specifically, based on the global index calculated by combining the interface curvature change rate and the vorticity intensity of the flow field, the invention achieves accurate quantification of the dynamic intensity of the system. According to the comparison result of this index with the preset threshold, the invention dynamically selects the adjustment strategy, so that the numerical method automatically enhances the spatiotemporal resolution during the rapid deformation stage of the interface, and reasonably reduces the computational accuracy during the stable flow field stage, thereby realizing the on-demand allocation of computational resources.

[0040] In critical stages such as droplet approach, film drainage, and interface fusion, this invention uses a solver reconfiguration mechanism to promptly employ smaller time steps and higher mesh resolution, ensuring accurate capture of interface topology changes. In slow evolution stages such as droplet rotation and equilibrium, a larger time step and a coarser mesh system are used, significantly reducing unnecessary computational consumption.

[0041] This invention forms a complete technical closed loop of solution-analysis-decision-reconstruction, breaking through the technical bottleneck of traditional methods that make it difficult to balance computational efficiency and numerical accuracy. This method can not only accurately reproduce the dynamic behavior of the entire droplet engulfment process, especially transient phenomena such as interface instability and film rupture, but also significantly improve the utilization efficiency of computing resources by autonomously optimizing discrete parameters, providing a more reliable numerical analysis tool for the design of microfluidic devices and the optimization of chemical processes. Attached Figure Description

[0042] Figure 1 This is a schematic diagram of the overall method flow of the present invention;

[0043] Figure 2 This is a bar chart of the global index of the present invention - the field norm of the rate of change of curvature;

[0044] Figure 3 This is a bar chart showing the field norm of curvature change rate versus vorticity field norm in this invention.

[0045] Figure 4 This is a curve showing the fitting of the field norm of the rate of change of curvature to the global index in this invention.

[0046] Figure 5 This is a fitting curve of the vorticity field norm-global index of the present invention;

[0047] Figure 6 This is the bubble color mapping diagram of the first weighting coefficient and the second weighting coefficient of the present invention. Detailed Implementation

[0048] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0049] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0050] Example:

[0051] Please see Figures 1-6 The present invention provides a technical solution:

[0052] A numerical method for the co-occurrence behavior of immiscible droplets in confined shear flow, comprising the following steps:

[0053] Step 1: Establish a simulation task of confined shear flow in the computer system. Use the phase-field equation coupled with the Navier-Stokes equation as the governing equation. Use the finite element method to spatially discretize the computational domain to generate a computational grid composed of several grid elements. Assign initial states to all grid elements according to the initial physical conditions to generate the initial phase field and velocity field. Set the initial time step and start the time step iteration process.

[0054] Establishing a simulation task specifically includes:

[0055] Establishing the phase-field equation coupled with the Navier-Stokes equation as the governing equation means using the Cahn-Hilliard equation to describe the evolution process of the immiscible fluid interface and the transient Navier-Stokes equation to describe the motion process of the background fluid. The two are coupled through the interfacial tension term. Spatial discretization adopts the finite element method, which means dividing the continuous computational domain into multiple non-overlapping grid elements to form an initial computational grid for numerical calculation.

[0056] The Cahn-Hilliard equation is the core of the phase-field method. Instead of directly tracking sharp interfaces, it introduces a spatially continuously varying phase-field variable (or order parameter), which has a value of 1 inside a droplet, -1 in the surrounding fluid, and a smooth transition at the interface. The equation describes the evolution of this phase-field variable over time. The properties of its higher-order differential terms allow the interface to naturally maintain a fixed and very small thickness and drive the interface to move in the direction of minimizing surface energy. It describes the evolution dynamics of the interface morphology.

[0057] The Navier-Stokes equations describe the conservation of momentum and mass of Newtonian fluids. They solve for the velocity vector and pressure scalar of the flow field, which determine the motion of the fluid under external forces such as shear, pressure, and interfacial tension. They describe the overall dynamics of the fluid.

[0058] Interfacial tension, as an internal force, is manifested as a volume force source term in the Navier-Stokes equations. The specific technique is as follows: First, the interfacial tension volume force acting on the fluid is obtained by calculating the gradient of the chemical potential (the core variable of the Cahn-Hilliard equation) from the current phase field distribution. Then, this volume force term is added to the right side of the momentum conservation equation of the Navier-Stokes equations. Subsequently, the velocity obtained from the flow field solution will, in turn, affect the transport of the phase field through the convection term. This is a typical two-way coupling mechanism. In numerical implementation, methods such as the continuous surface force model are used to smoothly apply this volume force to the computational unit near the interface.

[0059] Spatial discretization uses the finite element method to transform continuous partial differential equations (governing equations) into a system of algebraic equations that can be solved by a computer. This is because computers cannot directly process infinite continuous domains and must divide them into a finite number of small, simple parts.

[0060] The computational domain is filled with elements of triangular (2D) or tetrahedral (3D) shapes to ensure that all elements do not overlap and cover the entire computational domain. This type of element is highly adaptable to complex geometric boundaries. Mesh generation is completed by preprocessing software (such as Gmsh) or the built-in function of the solver. It requires presetting the geometry of the computational domain, such as the length, width, and height of the channels and mesh size control parameters, such as the global maximum / minimum element size, or local refinement at specific boundaries or regions. The software automatically generates a mesh file based on these parameters. The generated mesh file contains all coordinate information and the connection relationships of the elements, which constitutes the geometric basis of numerical computation.

[0061] Within each grid cell, the real physical fields (such as velocity, pressure, and phase field variables) are approximated using simple polynomial functions (shape functions). The governing equations are weighted and residual-processed on each cell using the Galerkin method or similar variational principles. This set of equations is then solved by a numerical solver.

[0062] The state includes a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The logic for generating the phase field and velocity field is as follows: based on the initial physical conditions, each grid cell is assigned a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The scalar values ​​representing the fluid phase identity in all grid cells are summed to generate the phase field, and the vector values ​​representing the fluid velocity in all grid cells are summed to generate the velocity field. The initial physical conditions include the initial flow field velocity distribution, the initial droplet position and size, and the physical property parameters of each fluid. The physical property parameters include density, viscosity, interfacial tension coefficient, mobility, and interfacial thickness parameters.

[0063] Density refers to the mass per unit volume of a fluid. In confined shear flow, the density difference between the droplet phase and the surrounding continuous phase is a key factor driving or influencing fluid motion (such as sedimentation and buoyancy effects) and interfacial instability. Its value is obtained by consulting physical chemistry handbooks or by experimental measurement (such as a densitometer), and is read as a constant in the input file of the program.

[0064] Dynamic viscosity is a measure of a fluid’s ability to resist flow (shear deformation). Viscosity ratio greatly affects the droplet’s ability to deform in shear flow, the intensity of internal circulation, and the dynamics of the engulfment process. High viscosity will inhibit interfacial deformation and flow instability. Its value can be obtained by consulting physical chemistry handbooks or by rheometer experiments.

[0065] The interfacial tension coefficient refers to the energy per unit length of the interface separating two immiscible fluids. It is the most direct driving force for interface evolution, tending to minimize the interfacial area, thereby keeping the droplets spherical or promoting fusion. This coefficient is the core bridge connecting the phase-field model and physical reality. Its value can be obtained by consulting physical chemistry handbooks or by measuring it with special experimental equipment (such as the pendant drop method or the rotating drop method).

[0066] In the phase-field model, the interface is treated as a diffusion layer with a finite thickness. This parameter controls the thickness of the interface. It must be much smaller than all relevant physical length scales (such as the droplet radius), but large enough to be resolved on the computational grid. It is set based on numerical analysis, requiring the finite thickness to be greater than the grid size to ensure interface resolution, but at the same time, the finite thickness is much smaller than R (where R is the droplet radius) to ensure that the interface is physically sharp. Its specific value is related to the interfacial tension and the energy barrier height in the phase-field model.

[0067] In the Cahn-Hilliard equations, mobility is a kinetic coefficient that controls the rate at which the phase field variables tend toward equilibrium. It affects the dynamic relaxation process of the interface motion and is set based on a trade-off between numerical stability and physical accuracy. Too small a mobility will result in slow interface dynamics, deviating from physical reality, while too large a mobility will increase the rigidity of the equations, requiring extremely small and long synchronization, which is computationally expensive. Its setting depends on asymptotic analysis to ensure that the correct interface tension effect and interface dynamics can be recovered when the interface thickness tends to zero.

[0068] For the phase field, based on the initial droplet position and size, a signed distance function or a hyperbolic tangent function is defined on the computational grid. Specifically, for a droplet located at... A circular droplet of radius R is located at point R, and each coordinate is a point on the circle. calculate ,in It is a parameter for controlling the interface thickness, so that inside the droplet... ,external The interface has a smooth transition, which is the scalar value representing the fluid phase identity;

[0069] For the velocity field, it is set according to the initial flow field velocity distribution. For example, for a simple shear flow, it is assigned... ,in It is the preset shear rate. It is a vertical coordinate, which is the vector value representing the velocity of the fluid.

[0070] The values ​​assigned above directly constitute the initial solution vector of the system of algebraic equations to be solved. The so-called phase field and velocity field are two (or one containing all degrees of freedom) large arrays in computer memory, which store the phase field scalar value and the velocity vector value respectively.

[0071] The initial flow field velocity distribution, as described above for shear flow distribution, requires a preset shear rate. The initial droplet position and size require pre-setting the coordinates of the droplet center. And radius R; the physical property parameters of each fluid need to be preset, such as density, viscosity, interfacial tension coefficient (used to calculate the magnitude of interfacial tension volume force), and phase field model related parameters such as mobility and interfacial thickness parameters.

[0072] Set the initial time step, start the time step iteration process, that is, start the simulation clock, and begin solving in chronological order; the initial time step is estimated based on numerical stability criteria such as CFL conditions.

[0073] Step 2: Perform coupled numerical solution based on the current time step and current grid resolution to obtain the updated phase field numerical solution and velocity field numerical solution for each grid cell. Calculate the interface curvature change rate and flow field vorticity intensity of each grid cell, and generate a global index characterizing the dynamic intensity of the entire system through weighted aggregation.

[0074] The current time step determines the accuracy of time advancement. The smaller the step, the more accurately it captures fast physical processes, but the higher the computational cost. The current grid resolution determines the ability to capture spatial details. The denser the grid, the higher the resolution for fine interface structures and flow field gradients, but the greater the computational cost.

[0075] For coupled numerical solutions, a step-by-step solution algorithm is used to solve the phase field control equations and the fluid dynamics equations sequentially within each time step, so as to obtain the updated phase field numerical solution and velocity field numerical solution defined on all grid cells.

[0076] The phase field equations (Cahn-Hilliard equations, fourth-order nonlinearity) and Navier-Stokes equations (second-order nonlinearity) are coupled together, making direct simultaneous solution very difficult and computationally expensive. Step-by-step solution is an efficient numerical strategy that decomposes the complex coupled problem into several relatively simple subproblems and solves them sequentially in each time step.

[0077] The process of solving the phase field equations involves inputting the velocity field and the phase field. With the velocity field kept constant, the new phase field is calculated by solving the Cahn-Hilliard equations. This equation determines how the interface is stretched and deformed by the flow field, as well as the diffusion and sharpening tendencies of the interface itself.

[0078] The process of solving the fluid dynamics equations involves inputting the velocity field and the newly calculated phase field, using the new phase field to calculate the interfacial tension, adding it as a source term to the momentum equation, and then solving the Navier-Stokes equations to obtain the new velocity field.

[0079] The output defines the new velocity field and the new phase field updated on all grid cells, which serve as the data basis for subsequent analysis.

[0080] The specific process of generating global metrics includes:

[0081] Based on the numerical solution of the phase field at the current time step, the curvature field is obtained by calculating its second spatial derivative, and the curvature rate of change field is obtained by subtracting it from the curvature field at the previous time step; based on the numerical solution of the velocity field at the current time step, the vorticity field is obtained by calculating its curl.

[0082] The formula for calculating the curvature field is as follows: The formula first calculates the phase field. Spatial gradient Then, the gradient vector is normalized to obtain the interface unit normal vector field. Finally, the divergence of the unit normal vector field is calculated. The result is the curvature field. ; It reflects the curvature of the interface in space; for example, a smaller, more rounded droplet will have a greater curvature value. The rate of spatial change, i.e., the gradient The larger and more concentrated the curvature, the higher the calculated curvature. The larger the absolute value, the better; positive curvature indicates a convex phase 1 at the interface (e.g., a droplet), while negative curvature indicates a concave phase 1 at the interface. The larger the absolute value, the more severe the interface curvature and the smaller the local radius of curvature. For example, the curvature can be very large at the neck of a droplet or at a filament that is about to break.

[0083] The formula for calculating the rate of change of curvature field is: , set the current time step Calculated curvature field Compared with the previous time step Stored curvature field Subtract, then divide by the time step between the two time steps. , where n is the index of the time step;

[0084] The instantaneous rate of interface shape change was quantified. A larger value indicates that the interface is undergoing rapid deformation, stretching, contraction, or wrinkling; for example, this value increases sharply at the moment a droplet is about to break up or merge. (Curvature field) The faster the change in a quantity over time t, the higher the calculated rate of change. The larger it is, the more it directly captures the dynamic characteristics of the interface motion, rather than the static geometric characteristics;

[0085] The formula for calculating the vorticity field is: For the velocity field To find the curl in three dimensions, It is also a vector; in two dimensions, it degenerates into a scalar (the component perpendicular to the plane). The rotational intensity and direction of the fluid micro-element were quantified. The larger the value, the more intense the fluid rotation at that location. In shear flow, high vorticity regions are located where the velocity gradient is large; in droplet wakes, vortices will form; velocity field Spatial variation (i.e., velocity gradient) The larger and more uneven the vorticity, the more accurate the calculated vorticity. The larger the vortex, the better; in a uniform flow, the vortex is zero.

[0086] The L2 norm of the curvature change rate field over the entire computational domain is calculated and weighted using a preset first weighting coefficient. The L2 norm of the vorticity field over the entire computational domain is also calculated and weighted using a preset second weighting coefficient. The two weighted results are then summed to obtain the global index. ;

[0087] The first weighting coefficient and the second weighting coefficient are determined based on the sensitivity analysis of each physical property parameter, and satisfy the condition that the sum of the first weighting coefficient and the second weighting coefficient is 1. The specific calculation formula is as follows:

[0088] ;

[0089] in, Let L2 norm be the field of curvature change rate. Let L2 norm be the vorticity field in the computational domain. The square root is taken after integrating the square of each physical quantity. In the discrete finite element model, this integration is achieved by summing the contributions of all mesh elements.

[0090] The L2 norm reflects how a local physical quantity field distributed throughout the entire field is aggregated into a single scalar value that characterizes the overall strength of that physical quantity across the entire field. A large value indicates that the entire system's interface is undergoing drastic and widespread deformation. A large value indicates that the entire flow field is under intense rotational or shearing motion. The L2 norm is mathematically similar to calculating the energy of a field, taking into account both the magnitude of the intensity and the range of its influence.

[0091] Global Indicators It characterizes the intensity of the dynamic evolution of the entire system (interface + flow field) at the current moment. The larger the value, the more dynamic, unstable, and rapidly changing the system is, requiring high computational precision to capture the moments of physical detail.

[0092] The two were determined through parameter sensitivity analysis. The specific method was to run a series of test cases under typical operating conditions and observe them individually. and The magnitude and importance of changes before and after critical events (such as droplet breakage) can be adjusted through trial and error or algorithm optimization. and This makes the global indicators It can most sensitively and accurately indicate the occurrence of these key events in advance or simultaneously; for example, if interface deformation is found to be the dominant mechanism of the consuming process, then it can assign... Larger values ​​(e.g., 0.7); if the vorticity variation of the background shear flow is also important, then set... , .

[0093] Specific data for partial time steps and global metrics are shown in Table 1.

[0094] Table 1 Global Indicator Statistics Chart

[0095]

[0096] Through data analysis, it was observed that there is a clear synergistic relationship between different characteristic parameters. For example, the data shows that the magnitude of the global index is simultaneously affected by the field norm of the rate of curvature change and the field norm of vorticity, and the two contribute to different degrees through weighting coefficients.

[0097] When analyzing the relationship between the field norm of the rate of change of curvature and the global index, the global index tends to show a higher value when the field norm of the rate of change of curvature is large and its corresponding weight coefficient is high. For example, in time step 7, the field norm of the rate of change of curvature reaches 3.0, and the first weight coefficient is 0.84, which makes the global index rise to a relatively high level of 3.19. Conversely, in time step 3, although the field norm of the rate of change of curvature reaches a relatively large value of 3.1, the global index is only 1.42 because its first weight coefficient is only 0.18.

[0098] Meanwhile, the influence of the vorticity field norm on the global index also shows a similar pattern. In time step 15, the vorticity field norm reaches a relatively large value of 4.5, and its second weighting coefficient is 0.33, which causes the global index to rise to 2.48. In time step 11, although the second weighting coefficient reaches a relatively high level of 0.77, the global index is only 1.48 because the vorticity field norm is only a relatively small value of 1.1.

[0099] It is particularly noteworthy that when either of the two norms reaches an extreme value and has a high corresponding weight, it will significantly affect the numerical trend of the global index. For example, in time step 4, although the vorticity field norm reaches a relatively large value of 3.5, its second weight coefficient is only 0.23, while the curvature change rate field norm is at a low level of 0.8, resulting in a global index of only 1.41. This characteristic indicates that when assessing the dynamic intensity of the system, it is necessary to comprehensively consider the balance between the actual values ​​of each physical quantity and their relative importance.

[0100] Through this multi-parameter collaborative analysis, we can more accurately grasp the dynamic characteristics of the system and provide a reliable basis for subsequent grid adaptive strategies. In practical applications, the weight coefficients should be dynamically adjusted according to the changing characteristics of each physical quantity under different working conditions in order to achieve the optimal allocation of computing resources.

[0101] Step 3: Compare the global index with the preset threshold. Based on the comparison result, dynamically select a set of adjustment strategies for the next time step. The adjustment strategy is to use a smaller time step and a higher grid resolution level when the global index is higher than the preset threshold, and a larger time step and a lower grid resolution level when the global index is lower.

[0102] The preset thresholds include preset high and low thresholds, and the adjustment strategies include high-precision parameter strategies, high-efficiency parameter strategies, and maintenance parameter strategies. The logic for dynamically selecting the adjustment strategy is as follows:

[0103] High threshold and low threshold This creates a hysteresis interval. This design avoids frequent policy oscillations caused by small fluctuations of the global metric around a single threshold, thus enhancing the stability of the system. For example, if there is only one threshold, when the global metric fluctuates slightly above and below it, the time step and grid will keep switching, generating numerical noise.

[0104] The adjustment strategy reflects the system's decision on the allocation of computational resources for the next time step. It directly indicates which operating mode the solver should operate in: whether to enter fine capture mode, fast advance mode, or maintain the current mode.

[0105] The global index is compared and analyzed with the preset high and low thresholds: when the value of the global index is greater than the high threshold, a high-precision parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset refinement coefficient less than 1, and the grid resolution level is adjusted by adding a positive integer number of level adjustment amounts to the current level.

[0106] When the value of the global metric is less than the low threshold, a high-efficiency parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset coarsening coefficient greater than 1, and the grid resolution level is adjusted by reducing the current level by a positive integer number of level adjustments.

[0107] When the value of the global metric is between the high and low thresholds, a parameter maintenance strategy is selected for the next time step to maintain the current time step size and grid resolution level unchanged.

[0108] The formula for adjusting the time step under the high-precision adjustment strategy is as follows: ,in The formula for adjusting the time step under the high-efficiency adjustment strategy is: ,in The formula for adjusting the time step under the maintenance adjustment strategy is: ;

[0109] The new time step size set for the next time step. Smaller size means that simulations will use smaller time fragments to advance, more accurately capturing fast transient processes and meeting the requirements of numerical stability conditions (such as CFL conditions), but the computational cost is higher. Increasing the size means that the simulation will use larger step sizes to accelerate time progression and improve computational efficiency, but it will lose rapidly changing physical details.

[0110] For example, the current step size It forms the basis for adjustments, with preset refinement coefficients. It is a multiplier less than 1, such as 0.5 or 0.8. The smaller its value, the more drastically the time step is contracted, resulting in a more aggressive improvement in accuracy, but also a greater sacrifice in efficiency; preset coarsening factor. It is a multiplier greater than 1, such as 1.5 or 2.0. The larger its value, the more the time step expands, the more significant the improvement in efficiency, but the higher the risk of decreased accuracy.

[0111] The formula for adjusting the grid resolution level under the high-precision adjustment strategy is as follows: ,in The formula for adjusting the grid resolution level under the high-efficiency adjustment strategy is as follows: ,in The adjustment formula for the grid resolution level under the maintenance adjustment strategy is as follows: ;

[0112] This is the new grid resolution level set for the next time step. It is an integer index; the higher the level, the denser the corresponding grid. Increasing the size means that a finer mesh will be generated overall, which can distinguish smaller flow structures and interface features in space and has higher computational accuracy, but the amount of computation (memory and CPU time) will increase exponentially. Reducing the size means using a sparser grid overall, sacrificing spatial detail in exchange for a significant increase in computational speed;

[0113] The value is determined by the current grid resolution level. With level adjustment amount The addition and subtraction operations determine the level adjustment amount. It is a positive integer that determines the step size of the grid resolution change. This indicates that each adjustment is made to one level, resulting in a smooth and stable change. (If set...) If the adjustment is more radical, it will cause drastic changes in the numerical solution;

[0114] The preset refinement coefficient, preset coarsening coefficient, and level adjustment amount are preset based on numerical stability requirements during simulation initialization.

[0115] The preset refinement and coarsening coefficients must satisfy the stability constraints of the explicit or semi-implicit time integration method; for example, for convection-diffusion problems, the CFL condition must be met. High threshold values ​​are also required. and low threshold It needs to be calibrated through benchmarking, for example, by running a typical example with known results, such as droplet equilibrium and engulfment states, and observing the numerical range of the global exponent. Set to occur before a drastic change event (such as membrane rupture). When the system clearly enters a stable state, the level adjustment amount is set to 1 to ensure that the grid changes are not too drastic and to avoid large interpolation errors caused by abrupt changes in grid topology.

[0116] Step 4: Automatically reconstruct the discrete configuration of the numerical solver according to the selected adjustment strategy, that is, divide the computational grid by applying the time step and grid resolution level. The iterative process returns to Step 2, and performs coupled numerical solution for the next time step according to the selected adjustment strategy. This loop continues until the simulation task is completed, and finally outputs serialized field data for analyzing the swallowing behavior of immiscible droplets.

[0117] For the discrete configuration of the automatic reconfiguration numerical solver, the following logical steps are performed according to the selected tuning strategy:

[0118] 1) By calling the parameter setting interface provided by the numerical solver, the time step size of the next time step is set to a new time step size determined based on the selected adjustment strategy;

[0119] In the program, this is manifested as modifying a variable that controls the progression of time. For example, if there is a variable `double dt;` in the code, this step is to execute the operation `dt = new_dt;`, where `new_dt` is calculated in step 3. ;

[0120] If a commercial or open-source solver is used, such as ANSYS Fluent, OpenFOAM, or FEniCS, it provides a set of application programming interfaces (APIs) for program calls to dynamically modify parameters at runtime. This ensures that the solver uses a dynamically adjusted step size when advancing to the next time step, which is a direct means of controlling the accuracy of time advancement.

[0121] 2) Trigger the mesh adaptation function of the numerical solver to identify regions of interest based on the numerical solutions of the velocity field and phase field at the current time step. Specifically: based on the numerical solution of the phase field, calculate its spatial gradient magnitude distribution and mark all mesh elements with gradient magnitudes greater than a preset interface gradient threshold as high-dynamic regions of interest; based on the numerical solution of the velocity field, calculate its vorticity magnitude distribution and mark all mesh elements with vorticity magnitudes less than a preset vorticity intensity threshold as low-dynamic regions of interest.

[0122] The numerical solution of the phase field at the current time step is calculated. For each grid cell, the spatial gradient of its phase field variables is calculated, and the magnitude of this gradient is further calculated. This magnitude reaches its maximum in the interface region and approaches zero in the pure liquid phase and the pure background fluid phase. The larger the gradient magnitude, the more drastic the interface change at that location, and the more likely it is to be the core region of the interface. The gradient magnitude of each grid cell is compared with a preset interface gradient threshold. When the gradient magnitude of a cell is greater than the threshold, it is marked as a region of high dynamic interest. This interface gradient threshold is mainly used to filter numerical noise and accurately capture the core grid cells that constitute the interface, ensuring that high computational resolution is only allocated to the regions where physical quantities change most drastically.

[0123] The numerical solution of the velocity field at the current time step is calculated. For each grid cell, the vorticity of its velocity field is calculated, and the magnitude of the vorticity is further calculated. The vorticity magnitude reflects the region of near-irrotational, steady flow in the flow field. It indicates where computational resources can be saved because the flow structure in these regions is simple and does not require high resolution to characterize. The smaller the vorticity magnitude, the weaker the fluid rotation and the smoother the flow at that location. The vorticity magnitude of each grid cell is compared with a preset vorticity intensity threshold. When the vorticity magnitude of a cell is less than the threshold, it is marked as a low-dynamic region of interest. By setting an appropriate vorticity intensity threshold, regions with smooth flow can be reliably identified. These regions have simple flow structures and lower requirements for grid resolution, and can be marked as low-dynamic regions of interest where computational resources can be optimized.

[0124] 3) Based on the selected adjustment strategy, perform corresponding mesh topology transformations on the region of interest: when a high-precision parameter strategy is selected, refine the mesh cells in the high-dynamic region of interest by performing a mesh resolution level refinement operation; when a high-efficiency parameter strategy is selected, coarsen the mesh cells in the low-dynamic region of interest by performing a mesh resolution level coarsening operation; when a maintenance parameter strategy is selected, keep all regions of the current computational mesh unchanged.

[0125] High dynamic interest regions are identified by calculating the phase field gradient magnitude of each grid cell and comparing it with a preset interface gradient threshold. That is, all grid cells with gradient magnitudes greater than the preset interface gradient threshold are marked as high dynamic interest regions. Low dynamic interest regions are identified by calculating the vorticity magnitude of each grid cell and comparing it with a preset vorticity intensity threshold. That is, all grid cells with vorticity magnitudes less than the preset vorticity intensity threshold are marked as low dynamic interest regions.

[0126] The refinement operation under the high-precision parameter strategy involves the following steps when the high-precision parameter strategy is selected: The system automatically reads the preset grid resolution level increase from the high-precision parameter strategy, and within the high-dynamic-range interest area, adjusts the current grid resolution level... Upgraded to (Addition) According to the new resolution level requirements, the mesh cells in the high dynamic interest area are refined. Specifically, the parent cell is divided into multiple sub-cells with smaller geometric sizes, and the physical field data on the original coarse mesh is transferred to the newly generated fine mesh through a conservation interpolation algorithm.

[0127] The coarsening operation under the high-efficiency parameter strategy involves the following steps when the strategy is selected: The system reads the preset mesh resolution reduction amount from the high-efficiency parameter strategy and, within the low dynamic interest area, reduces the current mesh resolution level. Reduce to (Subtraction) According to the new resolution level requirements, the mesh cells in the low dynamic interest area are coarsened. Specifically, multiple adjacent small cells are merged into a parent cell with a larger geometric size, and the physical field data on the atomic cells are integrated into the newly generated coarse mesh through a volume weighted average algorithm.

[0128] The mesh adaptation module that triggers the solver involves refinement, coarsening, and field variable interpolation. The refinement operation divides a parent element within a marked region into multiple smaller sub-elements, such as dividing a triangle into four smaller triangles (regular subdivision) or a tetrahedron into eight smaller tetrahedrons. The coarsening operation merges multiple adjacent small elements within a marked region into a larger element. After mesh transformation, field variable interpolation maps the velocity field and phase field on the original mesh to the new mesh using a conservation interpolation algorithm to ensure computational continuity.

[0129] The high-precision parameter strategy refines the highly dynamic areas of interest, accurately allocating computing resources to the interface area and improving the accuracy of interface capture; the high-efficiency parameter strategy coarsens the low-dynamic areas of interest, saving computing resources in unimportant areas and improving overall efficiency; the maintenance parameter strategy has no operation, maintaining the stability of the current computing state.

[0130] Refining the high dynamic interest region can provide higher spatial resolution in the interface region, better capture key physical quantities such as curvature and interface tension, reduce the numerical dissipation of interface capture, avoid non-physical interface diffusion, and ensure accurate simulation of interface dynamics and morphological evolution.

[0131] Coarsening the region of interest with low dynamics can significantly reduce the total number of grid cells, reduce the computation time per time step, reduce computer memory requirements, and greatly improve computational efficiency while ensuring overall accuracy.

[0132] The simulation task is terminated based on the following logic. The simulation task is considered complete when any of the following conditions are met:

[0133] 1) The accumulated physical simulation time since the start of the simulation reaches the preset simulation duration limit;

[0134] This is a safety feature to prevent the simulation from looping indefinitely; the upper limit is preset based on the timescale of the physical process being studied.

[0135] 2) The following three sub-conditions must be met simultaneously:

[0136] 2.1) The droplet agglomeration process is completed. By performing connectivity analysis on the current phase field numerical solution, the droplet agglomeration process is considered complete when it is determined that there is only one connected droplet phase region in the computational domain.

[0137] The phase field data is binarized, for example, phase field data greater than 0.5 is regarded as droplets. Then, a connected component labeling algorithm, such as depth-first search or disjoint-set data structure algorithm, is used to traverse all grid cells and count the number of connected droplet regions. If the count is 1, it indicates that multiple droplets have merged into one, and the merging is completed.

[0138] 2.2) Droplet interface motion convergence: The droplet interface motion is considered to have converged when the maximum change in the position of the droplet centroid over N consecutive time steps is less than a preset position convergence threshold.

[0139] A centroid position queue of length N is maintained in memory. At each step, the current centroid is calculated, and the difference between it and the previous N-1 centroids is calculated. The maximum value is taken. The maximum position change reflects whether the overall spatial motion of the droplet has stopped. The smaller the value (less than the preset position convergence threshold), the more stable the droplet position is, indicating that the large-scale macroscopic motion has stopped. The centroid sequence is the independent variable. The greater the difference between these positions, the greater the calculated maximum change.

[0140] 2.3) System dynamic convergence: The system is considered to have reached dynamic convergence when the difference between the maximum and minimum values ​​of the global index over N consecutive time steps is less than a preset index convergence threshold.

[0141] Maintain a queue of global metrics of length N in memory, calculate the difference between the maximum and minimum values ​​in the queue, with the independent variable being the global metrics over N consecutive time steps, and the dependent variable being the range of these N metric values ​​(the difference between the maximum and minimum values). The larger the range, the more dynamic the system is still experiencing fluctuations; the smaller the range (less than the metric convergence threshold), the more stable the overall dynamic level of the system has become, indicating that the drastic changes that drive the mesh adaptation have disappeared.

[0142] Where N is the preset number of time steps greater than 1; N is a statistical window that uses multiple consecutive time steps for judgment to avoid misjudgment due to small fluctuations in single-step calculation; N of 1 will misjudge due to single-step fluctuations, and the larger N is (such as 5 or 10), the more reliable the judgment, but it will also delay termination.

[0143] The iterative process returns to step 2, whose operation logic is as follows:

[0144] After completing the discrete configuration reconstruction in step 4, the system automatically advances the computation state to the next time step, and based on the reconstructed time step and grid resolution, executes the coupled numerical solution process in step 2 again, thus forming a closed loop of time stepping.

[0145] When the simulation task is completed, the system decodes and encodes the velocity field and phase field values ​​corresponding to all time steps into binary data streams according to the time step sequence, and writes them to the storage device to form a single data file.

[0146] Serialization refers to converting complex data structures in memory (such as grid information and field data at each time step) into a stable byte stream for storage or network transmission; binary data streams use binary format instead of text format, which has the advantages of fast read and write speed and small file size, and is crucial for massive simulation data;

[0147] The entire field data at each time step is organized chronologically and encoded into an efficient binary format (such as HDF5 or VTK) to generate a complete database that can be used for post-processing. This allows users to replay the entire merging process or extract flow field and interface information at any time for quantitative analysis.

[0148] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0149] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0150] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0151] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A numerical method for the co-occurrence behavior of immiscible droplets in confined shear flow, characterized in that, The specific steps include: Step 1: Establish a simulation task of confined shear flow in the computer system. Use the phase-field equation coupled with the Navier-Stokes equation as the governing equation. Use the finite element method to spatially discretize the computational domain to generate a computational grid composed of several grid elements. Assign initial states to all grid elements according to the initial physical conditions to generate the initial phase field and velocity field. Set the initial time step and start the time step iteration process. Step 2: Perform coupled numerical solution based on the current time step and current grid resolution to obtain the updated phase field numerical solution and velocity field numerical solution for each grid cell. Calculate the interface curvature change rate and flow field vorticity intensity of each grid cell, and generate a global index characterizing the dynamic intensity of the entire system through weighted aggregation. Step 3: Compare the global metrics with the preset thresholds. Based on the comparison results, dynamically select a set of adjustment strategies for the next time step. The adjustment strategy is to use a smaller time step and a higher grid resolution level when the global metrics are higher than the preset thresholds, and a larger time step and a lower grid resolution level when they are lower. Step 4: Automatically reconstruct the discrete configuration of the numerical solver according to the selected adjustment strategy, that is, divide the computational grid by applying the time step and grid resolution level. The iterative process returns to Step 2, and performs coupled numerical solution for the next time step according to the selected adjustment strategy. This loop continues until the simulation task is completed, and finally outputs serialized field data for analyzing the swallowing behavior of immiscible droplets. For coupled numerical solutions, a step-by-step solution algorithm is used to solve the phase field control equations and the fluid dynamics equations sequentially within each time step, so as to obtain the updated phase field numerical solution and velocity field numerical solution defined on all grid cells. The specific process of generating global metrics includes: Based on the numerical solution of the phase field at the current time step, the curvature field is obtained by calculating its second spatial derivative, and the curvature rate of change field is obtained by subtracting it from the curvature field at the previous time step; based on the numerical solution of the velocity field at the current time step, the vorticity field is obtained by calculating its curl. Calculate the L2 norm of the curvature change rate field over the entire computational domain and weight it using a preset first weighting coefficient. Calculate the L2 norm of the vorticity field over the entire computational domain and weight it using a preset second weighting coefficient. Sum the two weighted results to obtain the global index. The first weighting coefficient and the second weighting coefficient are determined based on the sensitivity analysis of each physical attribute parameter, and satisfy the condition that the sum of the first weighting coefficient and the second weighting coefficient is 1.

2. The numerical method for the co-occurrence behavior of immiscible droplets in a confined shear flow according to claim 1, characterized in that: Establishing a simulation task specifically includes: Establishing the phase-field equation coupled with the Navier-Stokes equation as the governing equation means using the Cahn-Hilliard equation to describe the evolution process of the immiscible fluid interface and the transient Navier-Stokes equation to describe the motion process of the background fluid. The two are coupled through the interfacial tension term. Spatial discretization adopts the finite element method, which means dividing the continuous computational domain into multiple non-overlapping grid elements to form an initial computational grid for numerical calculation. The state includes a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The logic for generating the phase field and velocity field is as follows: based on the initial physical conditions, each grid cell is assigned a scalar value representing the fluid phase identity and a vector value representing the fluid velocity. The scalar values ​​representing the fluid phase identity in all grid cells are aggregated to generate the phase field, and the vector values ​​representing the fluid velocity in all grid cells are aggregated to generate the velocity field. The initial physical conditions include the initial flow field velocity distribution, the initial droplet position and size, and the physical property parameters of each fluid. The physical property parameters include density, viscosity, interfacial tension coefficient, mobility, and interfacial thickness parameters.

3. The numerical method for the co-occurrence behavior of immiscible droplets in a confined shear flow according to claim 2, characterized in that: The preset thresholds include preset high and low thresholds, and the adjustment strategies include high-precision parameter strategies, high-efficiency parameter strategies, and maintenance parameter strategies. The logic for dynamically selecting the adjustment strategy is as follows: The global index is compared and analyzed with the preset high and low thresholds: when the value of the global index is greater than the high threshold, a high-precision parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset refinement coefficient less than 1, and the grid resolution level is adjusted by adding a positive integer number of level adjustment amounts to the current level. When the value of the global metric is less than the low threshold, a high-efficiency parameter strategy is selected for the next time step to adjust the time step size and grid resolution. The time step size is adjusted by multiplying the current time step size by a preset coarsening coefficient greater than 1, and the grid resolution level is adjusted by reducing the current level by a positive integer number of level adjustments. When the value of the global metric is between the high and low thresholds, a parameter maintenance strategy is selected for the next time step to maintain the current time step size and grid resolution level unchanged. The preset refinement coefficient, preset coarsening coefficient, and level adjustment amount are preset based on numerical stability requirements during simulation initialization.

4. The numerical method for the co-occurrence behavior of immiscible droplets in a confined shear flow according to claim 3, characterized in that: For the discrete configuration of the automatic reconfiguration numerical solver, the following logical steps are performed according to the selected tuning strategy: 1) By calling the parameter setting interface provided by the numerical solver, the time step size of the next time step is set to a new time step size determined based on the selected adjustment strategy; 2) Trigger the mesh adaptation function of the numerical solver to identify the region of interest based on the velocity field numerical solution and phase field numerical solution of the current time step. Specifically, based on the phase field numerical solution, calculate its spatial gradient magnitude distribution and mark all mesh cells with gradient magnitude greater than the preset interface gradient threshold as high dynamic regions of interest. Based on the numerical solution of the velocity field, the distribution of its vorticity modulus is calculated, and all grid elements with vorticity modulus values ​​less than a preset vorticity intensity threshold are marked as low dynamic interest regions. 3) Based on the selected adjustment strategy, perform corresponding mesh topology transformations on the region of interest: when the high-precision parameter strategy is selected, refine the mesh cells in the high-dynamic region of interest by performing a mesh resolution level refinement operation; when the high-efficiency parameter strategy is selected, coarsen the mesh cells in the low-dynamic region of interest by performing a mesh resolution level coarsening operation; when the maintenance parameter strategy is selected, keep all regions of the current computational mesh unchanged.

5. The numerical method for the co-occurrence behavior of immiscible droplets in a confined shear flow according to claim 4, characterized in that: The simulation task is terminated based on the following logic. The simulation task is considered complete when any of the following conditions are met: 1) The accumulated physical simulation time since the start of the simulation has reached the preset simulation duration limit; 2) The following three sub-conditions must be met simultaneously: 2.1) The droplet agglomeration process is completed. By performing connectivity analysis on the current phase field numerical solution, the droplet agglomeration process is considered complete when it is determined that there is only one connected droplet phase region in the computational domain. 2.2) Droplet interface motion convergence: The droplet interface motion is considered to have converged when the maximum change in the position of the droplet centroid over N consecutive time steps is less than a preset position convergence threshold. 2.3) System dynamic convergence: The system is considered to have reached dynamic convergence when the difference between the maximum and minimum values ​​of the global index over N consecutive time steps is less than a preset index convergence threshold. Where N is the preset number of time steps greater than 1.

6. The numerical method for the co-occurrence behavior of immiscible droplets in a confined shear flow according to claim 5, characterized in that: The iterative process returns to step 2, whose operation logic is as follows: After completing the discrete configuration reconstruction in step 4, the system automatically advances the computation state to the next time step, and based on the reconstructed time step and grid resolution, executes the coupled numerical solution process in step 2 again, thus forming a closed loop of time stepping. When the simulation task is completed, the system decodes and encodes the velocity field and phase field values ​​corresponding to all time steps into binary data streams according to the time step sequence, and writes them to the storage device to form a single data file.

Images