Method for calculating stabilized viscoelastic fluid-solid coupling problem based on overlapped grids

Abstract

The invention provides a method for calculating a stabilized viscoelastic fluid-solid coupling problem based on an overlapped grid, and relates to the technical field of fluid-solid coupling calculation, the method comprises the following steps: processing a flow field region and solid particles to obtain an interpolation boundary, a contribution unit and an interpolation unit, and carrying out flow field calculation on a background grid to obtain background grid flow field information; transmitting the background grid flow field information to the sub-grids to obtain the flow field information of the sub-grids; calculating by using the flow field information of the sub-grids through a finite element method to obtain updated flow field information of the sub-grids; transmitting the updated sub-grid flow field information to a background grid by using bilinear Lagrange interpolation to obtain updated background grid flow field information; based on the updated background grid flow field information, parameter iteration updating is carried out through calculation, and simulation data of the motion trail and the rotation state of the solid particles in the flow field area are obtained. According to the method, the problems of low precision and low efficiency of fluid-structure interaction calculation are solved.

Description

TECHNICAL FIELD

[0001] The present specification relates to the technical field of fluid-structure interaction calculation, in particular to a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids. BACKGROUND

[0002] Viscoelastic fluid-structure interaction (VFSI) problems are ubiquitous in many engineering and scientific applications. For example, in the biomedical field, blood is a typical non-Newtonian viscoelastic fluid, and the blood vessel wall is a viscoelastic soft tissue. The interaction between blood flow and the blood vessel wall is of great significance for the diagnosis and treatment of cardiovascular diseases. In the processing industry, simulating the coupling between polymer melt and the mold during processing such as extrusion, injection molding, and blow molding is crucial for optimizing the mold and improving product quality. In traditional methods, the structure is often treated as a pure elastic body, or the fluid is treated as a pure Newtonian fluid. With the increasing complexity of modern engineering structures, traditional analysis methods often fail to accurately capture the interaction between fluid and solid. Therefore, the study of fluid-structure interaction problems not only improves the safety and reliability of structures, but also provides theoretical support for the development of new materials and technologies.

[0003] In the traditional Arbitrary Lagrangian-Eulerian (ALE) method, the grid needs to be updated at each time step according to the motion of the object. If the quality of the updated grid is not good, grid reconstruction is needed. This not only increases the computational complexity and affects the efficiency of the calculation, but also introduces interpolation errors and affects the accuracy of the calculation. When using non-body-fitted grid methods to handle fluid-structure interaction problems, for example, the Immersed Boundary Method (IBM) uses Euler background grids and Lagrangian grids on the surface of the object to describe the fluid and solid, respectively. By treating the solid as an area immersed in the fluid, the motion of the fluid and the solid is coupled, simplifying the simulation process. However, in the case of strong fluid-structure interaction, there may be a problem of insufficient accuracy.

[0004] The overlapping grid method also has some limitations. The overlapping grid method relies on the interpolation process between grids. In areas with complex flow characteristics, interpolation errors will increase. In viscoelastic fluids, stress has a history-dependent and strong directional nature. Simple interpolation of variables between grids can easily smooth the fine stress field distribution, making it difficult to accurately represent the unique memory effect and normal stress difference phenomenon of viscoelastic fluids, affecting the physical authenticity of the simulation. When simulating complex flow fields such as flow separation and vortex shedding, the boundary conditions of the interpolation boundary will be complex. Inaccurate application of boundary conditions may lead to instability of the numerical solution, which affects its application ability in complex flow fields. SUMMARY

[0005] The present application provides a kind of based on the stabilized viscoelastic fluid-structure interaction problem calculation method of overlapping grid, solve the low precision and low efficiency problem existing in fluid-structure interaction calculation.

[0006] In order to achieve the above invention purposes, the technical scheme adopted by the present application is as follows: a kind of based on the stabilized viscoelastic fluid-structure interaction problem calculation method of overlapping grid, comprising: S1: the flow field region and solid particles are processed, and the overlapping background grid and sub-grid are obtained; S2: based on background grid and sub-grid, by analysis, interpolation boundary, contribution unit and interpolation unit are obtained; S3: based on interpolation boundary, contribution unit and interpolation unit, for viscoelastic fluid-structure interaction problem, by finite element method and stabilization method, flow field calculation is carried out on background grid, and background grid flow field information is obtained; S4: using bicubic spline interpolation, background grid flow field information is transmitted to sub-grid, and the flow field information of sub-grid is obtained; S5: using the flow field information of sub-grid, calculation is carried out by finite element method and stabilization method, and the updated flow field information of sub-grid is obtained; S6: using bilinear Lagrange interpolation, the updated flow field information of sub-grid is transmitted to background grid, and the updated background grid flow field information is obtained; S7: based on the updated background grid flow field information, the simulation data of the motion trajectory and the rotation state of solid particles in flow field region are obtained by parameter iterative updating through calculation, and the calculation of stabilized viscoelastic fluid-structure interaction problem is completed.

[0007] Further, the S1 includes: The flow field region is divided using grid spacing, and the background grid is obtained; Based on the geometric shape and spatial position of solid particles, sub-grid generation is carried out using sub-grid generation parameters, and the sub-grid is obtained.

[0008] Further, the S2 includes: Based on background grid and sub-grid, the outer boundary of sub-grid is directly defined as interpolation boundary, on the interpolation boundary, the element of sub-grid is used as interpolation unit, for receiving information from background grid, the element of background grid is used as contribution unit, for providing flow field data;Wherein, contribution unit judges whether the node is inside the object by calculating the distance from each node of background grid to the center of the object: if the distance from the node to the center of the object is less than or equal to the radius of the object, the background grid element where the node is located is regarded as inside the object, and is excluded, and then flow field data is provided.

[0009] Furthermore, the background grid flow field information is determined by a velocity correction equation, wherein the spatial discrete format expression of the velocity correction equation is: ; in, Represents the mass matrix, Indicates the first The velocity field updated at each time step This represents the intermediate velocity encountered during the solution process. Indicates the time step. Represents the Reynolds number. Represents the pressure gradient matrix. Indicates the first The pressure field after each time step update This represents the pressure field at the nth time step.

[0010] Furthermore, the flow field information of the sub-grid is determined by an interpolation formula, wherein the expression of the interpolation formula is: ; ; in, This represents the flow field information at the sub-mesh boundary nodes. This represents the x-axis coordinate of the node. This represents the y-coordinate of the node. This represents the flow field information of the node at the i-th row and j-th column in the contributing unit corresponding to this node. Let represent the improved form of the i-th uniform B-spline basis function. This represents the x-coordinate of the bottom-left corner node of the cell in the background mesh. This represents the y-coordinate of the bottom-left node of the cell in the background mesh. Indicates the grid size of the background grid.

[0011] Furthermore, the updated sub-grid flow field information is determined by the fluid control equations within the sub-grid, wherein the expression for the fluid control equations within the sub-grid is: ; in, Indicates speed, Indicates time, Indicates the velocity of solid particles. Indicates the pressure gradient. Indicates stress, Indicates viscosity ratio.

[0012] Furthermore, the expression for the bilinear Lagrange interpolation is: ; ; ; ; wherein, represents the updated flow field information of a node in the background mesh that needs to be updated by interpolation, represents the x-axis coordinate of the above-mentioned node, represents the y-axis coordinate of the above-mentioned node, represents the Lagrange interpolation coefficient, represents the flow field information of each node in the contribution cell corresponding to the above-mentioned node, corresponding to u, represents the x-axis coordinate of the i-th node in the contribution cell, represents the y-axis coordinate of the i-th node in the contribution cell, represents the Lagrange interpolation coefficient corresponding to the i-th node in the contribution cell.

[0013] Further, the expression of the simulation data is: ; ; ; ; ; ; ; wherein, represents mass, represents acceleration, represents gravity, represents the fluid force received by the solid particle, represents the stress tensor generated by the fluid at a point on the particle surface, represents the solid boundary, represents the outward normal vector of a point on the particle surface, represents the particle surface, represents the angular acceleration, J represents the moment of inertia, represents the resultant moment, represents the spatial coordinates of a point on the particle surface, represents the spatial coordinates of the center of gravity of the solid particle, P represents pressure, and I represents the unit tensor, represents the viscous stress tensor, represents the velocity of the solid particle at time t, denotes the velocity of the solid particle at time t, denotes denotes the acceleration of the solid particle at time t, denotes the acceleration of the solid particle at time t, denotes the time step, denotes denotes the spatial coordinates of the solid particle at time t, denotes denotes the of the solid particle at time t, denotes denotes the angular velocity of the solid particle at time t, denotes denotes the angular velocity of the solid particle at time t, denotes denotes the angular acceleration of the solid particle at time t, denotes denotes the angular acceleration of the solid particle at time t, denotes denotes the angle of rotation of the solid particle at time t, denotes denotes the angle of rotation of the solid particle at time t.

[0014] The present application has the beneficial effects that: the present application provides a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids, a complex flow region is divided into a plurality of sub-regions, the calculation grids of each sub-region exist independently, there is an overlapping relationship between the grids, and the flow field information is transmitted and coupled through interpolation on the boundary of the overlapping region. By using irregular grids as sub-grids around the object, on the one hand, the irregular grids can be very flexible to fit extremely complex geometric shapes, without the need for excessive simplification of the complex geometry, and the adaptability to complex geometric shapes is enhanced; on the other hand, the local encryption can be easily realized in the region around the object with complex flow characteristics, and higher local calculation precision is obtained with less total number of grids, and the precision and efficiency of the algorithm are improved. Directly aiming at the complexity and viscoelasticity of the material constitutive, the precision and efficiency problems existing in the existing fluid-structure interaction algorithm are solved, the applicability of the overlapping grid method under complex flow field conditions is improved, and it can be applied to a wider range of actual fluid-solid flow situations. (1) The irregular grids are used in the sub-grids, without the need for excessive simplification of the complex geometry, and the adaptability to complex geometric shapes is enhanced; and the local encryption can be easily realized in the region around the object with complex flow characteristics, and the precision and efficiency of the algorithm are improved. (2) The regular grids are used in the background grid, the grid nodes are arranged in order, the neighbor information of each node can be directly obtained through simple index calculation, without additional storage of topological relationship, and high-efficiency parallel computing can be easily realized, so that the solution of large-scale fluid-structure interaction problems becomes feasible. (3) The fluid-structure interaction model constructed by the present application can not only be applied to the fluid-structure interaction problem of Newtonian fluid, but also be extended to viscoelastic fluid, filling the gap of the overlapping grid method in practical application and popularization. BRIEF DESCRIPTION OF DRAWINGS

[0015] The present application will be further illustrated in the form of exemplary embodiments, which will be described in detail with reference to the accompanying drawings. These embodiments are not limiting, and in these embodiments, the same reference numbers represent the same structures, wherein: Figure 1 is an exemplary flow chart of a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids according to some embodiments of the present application; Figure 2 is an exemplary schematic diagram of bicubic spline interpolation according to some embodiments of the present application; Figure 3 is an exemplary schematic diagram of bilinear Lagrange interpolation according to some embodiments of the present application; Figure 4 is an exemplary schematic diagram of horizontal velocity cloud of Re 40 according to some embodiments of the present application; Figure 5 is an exemplary schematic diagram of horizontal velocity cloud of Re 200 according to some embodiments of the present application; Figure 6 is an exemplary schematic diagram of vorticity at different time of the circular particle settling process shown according to some embodiments of the present specification; Figure 7 is an exemplary schematic diagram of pressure cloud of two circular particles relative motion at t=24 shown according to some embodiments of the present specification; Figure 8 is an exemplary schematic diagram of pressure cloud of two circular particles relative motion at t=32 shown according to some embodiments of the present specification; Figure 9 is an exemplary schematic diagram of vorticity of two circular particles relative motion at t=24 shown according to some embodiments of the present specification; Figure 10 is an exemplary schematic diagram of vorticity of two circular particles relative motion at t=32 shown according to some embodiments of the present specification; Figure 11 is an exemplary schematic diagram of lift and drag coefficients variation of upper layer particle shown according to some embodiments of the present specification; Figure 12 is an exemplary schematic diagram of lift and drag coefficients variation of lower layer particle shown according to some embodiments of the present specification; Figure 13 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle 0° at t=24 shown according to some embodiments of the present specification; Figure 14 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle 0° at t=32 shown according to some embodiments of the present specification; Figure 15 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle 90° at t=24 shown according to some embodiments of the present specification; Figure 16 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle 90° at t=32 shown according to some embodiments of the present specification; Figure 17 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle -90° at t=24 shown according to some embodiments of the present specification; Figure 18 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle -90° at t=32 shown according to some embodiments of the present specification; Figure 19 is an exemplary schematic diagram of vorticity of two triangular particles with rotation angle 180° at t=24 shown according to some embodiments of the present specification; Figure 20 is an exemplary schematic diagram of vorticity of a triangular particle with two rotation angles of 180° at t = 32 according to some embodiments of the present specification; Figure 21 is an exemplary schematic diagram of horizontal velocity contours of a viscoelastic Poiseuille flow (Re = 0.01, Wi = 0.1, β = 0.59) according to some embodiments of the present specification; Figure 22 is an exemplary schematic diagram of vertical velocity contours of a viscoelastic Poiseuille flow (Re = 0.01, Wi = 0.1, β = 0.59) according to some embodiments of the present specification; Figure 23 is an exemplary schematic diagram of pressure contours of a viscoelastic Poiseuille flow (Re = 0.01, Wi = 0.1, β = 0.59) according to some embodiments of the present specification; Figure 24 is an exemplary schematic diagram of b xx contours of a viscoelastic Poiseuille flow (Re = 0.01, Wi = 0.1, β = 0.59) according to some embodiments of the present specification; Figure 25 is an exemplary schematic diagram of horizontal velocity contours of a viscoelastic cylinder flow (Re = 0.01, Wi = 0.1, β = 0.59) according to some embodiments of the present specification. DETAILED DESCRIPTION

[0016] The specific embodiments of the present application are described herein. It should be apparent, however, to those skilled in the art that many modifications might be made to the specific embodiments without departing from the spirit and scope of the application as set forth in the claims. Accordingly, the techniques described herein should be understood to be merely illustrative of specific applications and should not be construed as limiting.

[0017] Embodiment One Figure 1 is an exemplary flow chart of a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids according to some embodiments of the present specification. As shown in Figure 1 , the flow includes the following steps. In some embodiments, the flow can be performed by a processor.

[0018] S1: process the flow field region and the solid particle to obtain an overlapping background grid and a sub-grid.

[0019] The flow field region is a region around an object with flow characteristics. For example, the flow field region can include the flow field region around the inside of an engine, an airplane as a whole, a biological organ, and the like, a Poiseuille of a viscoelastic fluid, and a flow around a cylinder of a viscoelastic fluid, and the like. Specifically, the flow field region can include the range of the flow field region (for example, inlet, outlet, wall), and the grid spacing d.

[0020] The solid particle is a solid particle moving in the flow field region. For example, the solid particle can include a circular particle, a triangular particle, and the like. Specifically, the solid particle can include the geometric shape and the spatial position of the solid object (for example, the center coordinates (x cen, y cen) and the radius r of a circular object, the coordinates of three vertices of a triangle, or a more complex shape description), and the like.

[0021] The background grid is a uniformly distributed node generated in the entire calculation domain with a fixed grid spacing d. For example, the background grid can include the node coordinates of the background grid, the element connection information, and the boundary node coordinates of the inlet, outlet, and wall. The node coordinates node is a matrix, the element connection information element is a matrix (when Taylor Hood element is used, element1 (en x 3) and element2 (en x 6) two element connection information matrices are output respectively corresponding to P1 and P2 elements), and the boundary node coordinates

[0022] The sub-grid is a non-regular triangular grid generated by using the PDEToolbox tool in MATLAB based on the geometric shape of the object. For example, the sub-grid can include the sub-grid generation parameters (such as the shape size of the sub-grid, the grid density, and the like), the node coordinates of the sub-grid, the element connection information, and the solid boundary node coordinates. The node coordinates p is a n_sub x 2 matrix, the element connection information t is a en_sub x 3 matrix (when Taylor Hood element is used, t1 (en_sub x 3) and t2 (en_sub x 6) two element connection information matrices are output respectively corresponding to P1 and P2 elements), and the boundary node coordinates cy is a n_cy x 2 matrix.

[0023] In some embodiments, the processor can divide the flow field region by using the grid spacing to obtain the background grid; and generate the sub-grid by using the sub-grid generation parameters based on the geometric shape and the spatial position of the solid particle.

[0024] ​​​​​​In some embodiments, the processor can arrange the nodes in the background grid in a two-dimensional rectangular coordinate system, and form regular rectangular triangular grid cells by connecting adjacent nodes. The background grid covers the entire computational domain, including the area where the solid object can move, ensuring the integrity of the computational domain. The grid spacing d is pre-set according to the accuracy and efficiency requirements of the calculation, and is usually determined based on the characteristic scale of the flow field.

[0025] In some embodiments, the sub-grid closely fits the surface of the object and can adaptively refine the area near the object to capture complex flow details. The generation range of the sub-grid is extended to a large enough area to ensure that its outer boundary has sufficient overlap with the background grid. The width of the overlapping area is usually set to at least 5 grid cell sizes to ensure interpolation accuracy.

[0026] S2: Based on the background grid and the sub-grid, the interpolation boundary, the contribution cell and the interpolation cell are obtained through analysis.

[0027] The interpolation boundary is a boundary area for interpolation.

[0028] The contribution cell is a cell for information exchange.

[0029] The interpolation cell is a cell for interpolation.

[0030] In some embodiments, the processor can define the outer boundary of the sub-grid as the interpolation boundary based on the background grid and the sub-grid. On the interpolation boundary, the cells of the sub-grid are used as interpolation cells to receive information from the background grid, and the cells of the background grid are used as contribution cells to provide flow field data. The contribution cell determines whether a node is inside the object by calculating the distance from each node of the background grid to the center of the object: if the distance from the node to the center of the object is less than or equal to the radius of the object, the background grid cell where the node is located is considered to be inside the object and is excluded, and then the flow field data is provided.

[0031] In some embodiments, the processor can define the outer boundary of the sub-grid as an interpolation boundary directly based on the node coordinates node and the element connection information element of the background grid, the node coordinates p and the element connection information t of the sub-grid, that is, the interface of the sub-grid and the background grid for exchanging flow field information. On the interpolation boundary, the elements of the sub-grid serve as interpolation elements for receiving information from the background grid, and the elements of the background grid serve as contribution elements for providing flow field data. This method of directly selecting the boundary of the sub-grid simplifies the identification process of the interpolation boundary and improves the calculation efficiency; in the background grid, the elements occupied by the solid are not involved in the flow field calculation, so as to avoid invalid flow field calculation in these areas. Whether a node is inside the object is determined by calculating the distance from each node of the background grid to the center of the object: if the distance from the node to the center of the object is less than or equal to the radius (or equivalent radius) of the object, the background grid element where the node is located is considered to be inside the object and is excluded. This method is suitable for spherical or approximately spherical objects. For non-spherical objects, a more accurate judgment method based on the geometric shape (such as the ray casting method) can be used to obtain the outer boundary boundary of the sub-grid (the interpolation boundary of the sub-grid), the calculation boundary hole of the background grid, the contribution element information p_id of the boundary nodes of the sub-grid in the background grid, and the contribution element information R of the boundary nodes of the background grid in the sub-grid.

[0032] S3: Based on the interpolation boundary, the contribution element and the interpolation element, the flow field of the background grid is calculated by the finite element method and the stabilization method for the viscoelastic fluid-structure coupling problem to obtain the background grid flow field information.

[0033] The background grid flow field information is the updated flow field information (u, v, P, b, ξ) of the background grid.

[0034] The stabilization method includes SRCR, SUPG, DEVSS, PGP, etc.

[0035] In some embodiments, the background grid flow field information is determined by a velocity correction equation, wherein the spatial discrete format expression of the velocity correction equation is: ; wherein, represents a mass matrix, represents the updated velocity field at the th time step, represents an intermediate velocity appearing in the solving process, represents a time step, represents a Reynolds number, represents a pressure gradient matrix, represents the updated pressure field at the th time step, Pn represents the pressure field at the nth time step.

[0036] In some embodiments, the processor can solve the incompressible viscous flow problem based on the background mesh boundary conditions (e.g., as shown in FIG. 1, inlet uses VelocityInlet to give parabolic velocity, outlet uses PressureOutlet to specify the static pressure at the outlet, wall uses No-SlipWall to assume the fluid velocity close to the wall is the same as the wall velocity, hole uses the matching boundary condition, and the boundary nodes use the information from the previous time step), and initial conditions, such as u, v, P, b, Figure 4 , Figure 5 and Figure 25 , Re, dt (time step), Wi, , the processor can solve the NS equations using the finite element method to obtain the flow field information of the background mesh.

[0037] In some embodiments, the processor can solve the viscoelastic fluid problem to obtain the spatial discrete format of the velocity correction equation by the following way: The mass and momentum conservation equations of incompressible viscoelastic flow are given in the following dimensionless form:

[0038]

[0039] where u is the velocity, p is the pressure, is the stress, is the solvent viscosity coefficient, and Re is the Reynolds number. The unknown quantities we need to solve are u, p, and . The stress satisfies

[0040] where Wi is the Weissenberg number, c represents the conformation tensor, which in the Oldroyd-B model takes the following form:

[0041] The dimensionless parameters Re and Wi are defined as:

[0042] where is the fluid density, is the total viscosity, is the relaxation time, U is the characteristic velocity, L is the characteristic length.​

[0043] To ensure the conformation tensor c is always positive definite, introduce c the square root of the conformation tensor

[0044]

[0045]

[0046]

[0047]

[0048]

[0049] The non-zero components of the anti-symmetric matrix a have the following form:

[0050] The momentum equation is time-discretized by BDF2, and the constitutive equation is solved by second-order Runge-Kutta method. The flow field information at time step n is denoted by , , , and . The square root of the conformation tensor is first updated, and is defined as b

[0051] The intermediate value

[0052]

[0053] is then calculated. The intermediate stress is obtained according to .

[0054] The momentum conservation equation is discretized into BDF2 format:

[0055] The intermediate velocity is calculated according to the following format:

[0056] The velocity at time step n+1 is obtained by subtracting the two equations:

[0057] The mass conservation equation is given by: ​

[0058] Thus obtained Pressure at time and velocity Finally, the stress is updated using the second-order Runge-Kutta scheme:

[0059]

[0060] In the process of solving fluid-structure interaction problems of viscoelastic fluid, numerical oscillation often occurs, which leads to the failure of convergence. In this case, the SUPG formula is used to stabilize the dominant convection, and the test function is written as:

[0061] where is the stabilization parameter, denoted as

[0062] where represents the average velocity of a certain element, represents the local element size in the streamline direction. The following auxiliary matrix is defined:

[0063] We can get the matrix form of the discrete format (13) of

[0064] where

[0065]

[0066]

[0067] The DEVSS method is used to avoid the LBB constraint of velocity-stress interpolation, which maintains the elliptic property of the momentum equation. The discrete approximation of the deformation tensor rate is introduced:

[0068]

[0069] The approximation space of is usually the same as Since the same linear interpolation is used, is also approximated as a continuous first-order polynomial with other variables.

[0070] ; The right-hand side term has a stabilizing effect on the discrete level, while it is zero on the continuous level. To avoid iterations, we approximate the stabilizing term at time t by the values at the two previous time steps: ; The equation for the intermediate velocity ; We approximate D by

[0071] We introduce two auxiliary matrices:

[0072] The weak form of the above equation is obtained as

[0073]

[0074] where

[0075]

[0076]

[0077] To overcome the poor pressure stability in the incremental fractional scheme, we add an extra stabilizing term to the weak form of the pressure Poisson equation. Let denote the introduced auxiliary variable, which is approximated by

[0078] The matrix form of the pressure Poisson equation is now written as

[0079]

[0080]

[0081]

[0082] In the matrix form of the pressure Poisson equation, is a stabilizing parameter depending on the local element size and velocity. If , it is the standard Galerkin discretization. In general, it can be defined as ​​

[0083] In addition, corresponding to the two cases of formula (33) , Both explicit and implicit processing can be performed. Obviously, implicit format requires more calculation, but since the explicit format and implicit format have the same stability effect. Therefore, the explicit format is adopted.

[0084] Wherein, the parameter meanings of the above formula include: u Velocity, p represents pressure, Stress, Solvent viscosity coefficient, Re Reynolds number, Wi Weissenberg number, c Conformation tensor, Fluid density, Total viscosity, Relaxation time, U Characteristic velocity, L Characteristic length, Stability parameter, Average velocity of a certain unit, Local cell size in the streamline direction, Introduced auxiliary variable, A stability parameter depending on the local cell size and velocity, M, K, G, B represent matrix form expression.

[0085] S4: using bicubic spline interpolation, the background grid flow field information is transmitted to the sub-grid, and the flow field information of the sub-grid is obtained.

[0086] The flow field information of the sub-grid includes ; wherein, Horizontal velocity of the sub-grid node, Vertical velocity of the sub-grid node, Pressure of the sub-grid node, Square root of the stress tensor of the sub-grid node, Artificial variable of the sub-grid node for stabilizing the pressure term.

[0087] In some embodiments, the flow field information of the sub-grid is determined by an interpolation formula, wherein the expression of the interpolation formula is: ; ; Wherein, Flow field information (such as velocity field, pressure field, stress field) of the sub-grid boundary node,​ x-coordinate of the node, y-coordinate of the node, flow field information of the node at the i-th row and j-th column in the contribution cell corresponding to the node, improved form of the i-th uniform B-spline basis function, x-coordinate of the bottom-left node of the cell in the background mesh where the node is located, y-coordinate of the bottom-left node of the cell in the background mesh where the node is located, grid size of the background mesh.

[0088] In some embodiments, the processor can use bicubic spline interpolation to transfer the flow field information of the background mesh to the sub-mesh based on the flow field information (u, v, P, b, ), the contribution cell information of the sub-mesh outer boundary boundary (sub-mesh interpolation boundary) and the sub-mesh boundary node in the background mesh Bicubic spline interpolation is a high-order interpolation method based on uniform B-spline basis functions, which can construct a smooth surface in a two-dimensional plane, ensuring the accuracy and continuity of the flow field information during the transfer process, and obtaining the flow field information of the sub-mesh.

[0089] These basis functions ensure the smoothness and cubic precision (truncation error O(h^4), C^2 continuity) of the interpolation surface, effectively reducing numerical dissipation and oscillation in information transfer, and ensuring the stability and accuracy of flow field coupling. See the accompanying drawings for an interpolation diagram Figure 2 .

[0090] S5: Using the flow field information of the sub-mesh, the updated flow field information of the sub-mesh is obtained by finite element method and stabilization method.

[0091] The updated flow field information of the sub-mesh includes ; wherein, updated horizontal velocity of the sub-mesh node, updated vertical velocity of the sub-mesh node, updated pressure of the sub-mesh node, updated square root of stress tensor of the sub-mesh node, updated artificial variable of the sub-mesh node for stabilizing the pressure term.

[0092] In some embodiments, the updated flow field information of the sub-mesh is determined by the fluid control equation in the sub-mesh, wherein the expression of the fluid control equation in the sub-mesh is: ; wherein, Indicates speed, Indicates time, Indicates the velocity of solid particles. Indicates the pressure gradient. Indicates stress, Indicates viscosity ratio.

[0093] In some embodiments, such as Figures 6-24 The processor can base its work on various flow field information from the sub-mesh. The subgrid boundary conditions (the outer boundary uses a butt joint boundary condition to ensure the flow field information is maintained by the interpolation obtained in step 4, and the solid boundary cy uses a No-SlipWall) are also considered. In the bidirectional coupled model, when solid particles translate and rotate, the subgrid must also translate and rotate accordingly. Therefore, the updated subgrid flow field information is obtained using the fluid control equations described by ALE within the subgrid.

[0094] S6: Use bilinear Lagrange interpolation to transfer the updated subgrid flow field information to the background grid, and obtain the updated background grid flow field information.

[0095] The updated background mesh flow field information includes (u, v, P, b, ... );in, This indicates the horizontal velocity of the sub-mesh node after update. This represents the vertical velocity of the sub-mesh node after update. denoted by , where 'b' represents the pressure after the sub-mesh node update, and 'b' represents the square root of the stress tensor after the sub-mesh node update. This represents the artificial variable used to stabilize the pressure term after the subgrid nodes are updated.

[0096] In some embodiments, the expression for bilinear Lagrange interpolation is: ; ; ; ; in, This represents the updated flow field information of the nodes in the background mesh that need to have their flow field information updated via interpolation. This represents the x-axis coordinate of the above nodes. This represents the y-axis coordinate of the above nodes. Represents the Lagrange interpolation coefficients. This represents the flow field information of each node in the contributing unit corresponding to the above node, corresponding to u. This represents the x-axis coordinate of the i-th node in the contributing unit. y-coordinate of the i-th node in the contribution cell, Lagrange interpolation coefficient corresponding to the i-th node in the contribution cell.

[0097] In some embodiments, the processor can calculate the flow field information of the sub-grid based on the flow field information of the background grid, The background grid calculates the boundary hole and the contribution cell information R of the background grid boundary node in the sub-grid, and transmits the flow field information of the sub-grid to the background grid by bilinear Lagrange interpolation to obtain the updated flow field information of the background grid.

[0098] The bilinear Lagrange interpolation has a first order of accuracy (truncation error and continuity), and is simple and efficient in calculation, suitable for large-scale calculation and dynamic grid update. Although the accuracy is low, it can meet the stability and conservation requirements of the overall calculation when used for information feedback from the sub-grid to the background grid. The interpolation schematic diagram is shown in the accompanying drawings Figure 3 .

[0099] S7: Based on the updated flow field information of the background grid, the parameter iteration is updated by calculation to obtain the simulation data of the motion trajectory and rotation state of the solid particles in the flow field region, and the calculation of the stabilized viscoelastic fluid-structure coupling problem is completed.

[0100] The simulation data of the motion trajectory and rotation state of the solid particles in the flow field region includes the updated geometry and spatial position of the solid object, the node coordinates p of the sub-grid, the total force and the total moment exerted on the solid object at the current time step .

[0101] In some embodiments, Figures 6-10 The particle sedimentation shown in the figure needs to consider the gravity G of the particle itself; while Figures 11-24 The relative motion is discussed, which is a planar motion, and G does not need to be considered.

[0102] In some embodiments, the expression of the simulation data is: ; ; ; ; ; ; ; Among them, mass, acceleration, gravity, fluid force received by a solid particle, stress tensor generated by fluid at a point on the surface of a particle, solid boundary, outward normal vector at a point on the surface of a particle, surface of a particle, angular acceleration of rotation, J represents moment of inertia, resultant moment, spatial coordinates of a point on the surface of a particle, spatial coordinates of the center of gravity of a solid particle, P represents pressure, I represents unit tensor, viscous stress tensor, velocity of a solid particle at time t, velocity of a solid particle at time t, acceleration of a solid particle at time t, acceleration of a solid particle at time t, time step, spatial coordinates of a solid particle at time t, angular velocity of rotation of a solid particle at time t, angular velocity of rotation of a solid particle at time t, angular acceleration of rotation of a solid particle at time t, angular acceleration of rotation of a solid particle at time t, angle of rotation of a solid particle at time t, angle of rotation of a solid particle at time t.

[0103] In some embodiments, the processor can determine the flow field information of the sub-grid based on the flow field information of the sub-grid, the geometric shape and spatial position of the solid object, and the solid boundary cy, by assuming that the surface of the particle is a curved surface , the outward normal vector of which is . The stress tensor generated by the fluid at any point on the surface of the particle is , which is the joint action of pressure (normal) and shear stress (tangential), and can be represented as wherein,​​​​​​​​​​​ is a unit tensor, is pressure, is viscous stress tensor, the update of translational and rotational motion uses trapezoidal formula which shows time integration. The sub-grid is also updated synchronously with the movement and rotation of solid particles to get the updated geometry and spatial position of solid objects, the node coordinates p of the sub-grid, the total force and total moment experienced by the solid object in the current time step, and the velocity of the solid object in the current time step.

[0104] In some embodiments, the processor can determine whether a termination condition is reached. As Figure 4 , the termination condition is that the flow field reaches a steady state, i.e., the flow field information changes by less than a threshold value after a number of time steps; as Figures 6-10 , the termination condition is that the circular particle falls to the lower wall, i.e., the spatial coordinates of the particle reach the target position; as Figures 11-24 , the termination condition is that the cumulative time step reaches the requirement, in the embodiment shown in Figures 11-24 , when the number of simulation steps makes the physical time reach 32, the termination condition is reached. If the termination condition is reached, the program is ended; if the termination condition is not reached, the calculation is returned to step 2 until the result converges or the target time is reached.

[0105] In some embodiments of the present specification, a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids is provided, a complex flow region is divided into a plurality of sub-regions, the calculation grid of each sub-region independently exists, and there is an overlapping relationship between the grids. The flow field information is transmitted and coupled through interpolation on the boundary of the overlapping region. By using irregular grids as sub-grids around the object, on the one hand, the extremely complex geometric shape can be very flexibly fitted without excessive simplification of the complex geometry, and the adaptability to the complex geometric shape is enhanced. On the other hand, the local encryption in the region around the object with complex flow characteristics can be easily realized, and higher local calculation precision is obtained with less total number of grids, and the precision and efficiency of the algorithm are improved. Directly aiming at the complexity and viscoelasticity of the material constitutive, the precision and efficiency problems existing in the existing fluid-structure interaction algorithm are solved, the applicability of the overlapping grid method under complex flow field conditions is improved, and it can be applied to a wider range of actual fluid-solid flow situations. (1) The irregular grid is used in the sub-grid, the complex geometry does not need to be excessively simplified, and the adaptability to the complex geometric shape is enhanced. Local encryption in the region around the object with complex flow characteristics can be easily realized, and the precision and efficiency of the algorithm are improved. (2) The regular grid is used in the background grid, the grid nodes are arranged in order, and the neighbor information of each node can be directly obtained by simple index calculation without additional storage of topological relationship. Efficient parallel computing can be easily realized, and the solution of large-scale fluid-structure interaction problems becomes feasible. (3) The fluid-structure interaction model constructed by the present application can not only be applied to the fluid-structure interaction problem of Newtonian fluid, but also be extended to viscoelastic fluid, filling the gap of overlapping grid method in practical application and popularization.

[0106] Embodiment two In the field of biomedicine, blood is a typical non-Newtonian viscoelastic fluid, and the blood vessel wall is a viscoelastic soft tissue. The interaction between blood flow and blood vessel wall is of great significance for the diagnosis and treatment of cardiovascular diseases. However, in traditional methods, the structure is often regarded as a pure elastic body, or the fluid is regarded as a pure Newtonian fluid, which often cannot accurately capture the interaction between fluid and solid.

[0107] In some embodiments, a stabilized viscoelastic fluid-structure interaction problem calculation method based on overlapping grids is provided, comprising: S1: as Figures 21-24 , simulate the movement of the high polymer melt in the mold; as Figure 25 , simulate the microcellular injection molding scene of the high polymer melt, process the flow field region and the solid particles, and obtain the overlapping background grid and sub-grid; S2: based on the background grid and the sub-grid, through analysis, the interpolation boundary, the contribution element and the interpolation element are obtained; S3: Based on the interpolation boundary, the contribution unit and the interpolation unit, for the viscoelastic fluid-structure coupling problem, the flow field information of the background grid is calculated by the finite element method and the stabilization method, to obtain the background grid flow field information; S4: The background grid flow field information is transmitted to the sub-grid by using the bicubic spline interpolation to obtain the flow field information of the sub-grid; S5: The flow field information of the sub-grid is calculated by using the finite element method and the stabilization method to obtain the updated flow field information of the sub-grid; S6: The updated flow field information of the sub-grid is transmitted to the background grid by using the bilinear Lagrange interpolation to obtain the updated flow field information of the background grid; S7: Based on the updated flow field information of the background grid, the parameter iteration is updated by calculation to obtain the simulation data of the motion trajectory and the rotation state of the solid particles in the flow field region, and the calculation of the stabilized viscoelastic fluid-structure coupling problem is completed.

[0108] In some embodiments, the S1 comprises: The flow field region is divided by using the grid spacing to obtain the background grid; Based on the geometric shape and the spatial position of the solid particles, the sub-grid is generated by using the sub-grid generation parameters to obtain the sub-grid.

[0109] In some embodiments, the S2 comprises: Based on the background grid and the sub-grid, the outer boundary of the sub-grid is directly defined as the interpolation boundary, on which the unit of the sub-grid is used as the interpolation unit to receive the information from the background grid, and the unit of the background grid is used as the contribution unit to provide the flow field data; wherein the contribution unit judges whether the node is inside the object by calculating the distance from each node of the background grid to the center of the object: if the distance from the node to the center of the object is less than or equal to the radius of the object, the background grid unit where the node is located is considered to be inside the object, and is excluded to provide the flow field data.

[0110] In some embodiments, the background grid flow field information is determined by a velocity correction equation, wherein the spatial discrete format expression of the velocity correction equation is: ; Wherein, represents the mass matrix, represents the updated velocity field of the first time step, represents the intermediate velocity appearing in the solving process, represents the time step, represents the Reynolds number, represents the pressure gradient matrix, represents the first the updated pressure field at the time step, denotes the pressure field at the nth time step.

[0111] In some embodiments, the flow field information of the sub-grid is determined by an interpolation formula, where the expression of the interpolation formula is: ; ; wherein, denotes the flow field information (such as velocity field, pressure field, stress field) of the boundary node of the sub-grid, denotes the x-axis coordinate of the node, denotes the y-axis coordinate of the node, denotes the flow field information of the node at the i-th row and j-th column in the contribution element corresponding to the node, denotes the modified form of the i-th uniform B-spline basis function, denotes the x-axis coordinate of the lower left corner node of the element in the background grid where the node is located, denotes the y-axis coordinate of the lower left corner node of the element in the background grid where the node is located, denotes the grid size of the background grid.

[0112] In some embodiments, the updated flow field information of the sub-grid is determined by the fluid control equation in the sub-grid, where the expression of the fluid control equation in the sub-grid is: ; wherein, denotes the velocity, denotes the time, denotes the velocity of the solid particles, denotes the pressure gradient, denotes the stress, denotes the viscosity ratio.

[0113] In some embodiments, the expression of the bilinear Lagrange interpolation is: ; ; ; ; wherein, denotes the updated flow field information of the node in the background grid that needs to be updated by interpolation, denotes the x-axis coordinate of the node, denotes the y-axis coordinate of the node, denotes the Lagrange interpolation coefficient, represents the flow field information of each node in the contribution unit corresponding to the node, and corresponds to u, represents the x-axis coordinate of the i-th node in the contribution unit, represents the y-axis coordinate of the i-th node in the contribution unit, represents the Lagrange interpolation coefficient corresponding to the i-th node in the contribution unit.

[0114] In some embodiments, the expression of the simulation data is: ; ; ; ; ; ; ; wherein, represents mass, represents acceleration, represents gravity, represents the fluid force received by the solid particle, represents the stress tensor generated by the fluid at a point on the surface of the particle, represents the solid boundary, represents the outward normal vector of a point on the surface of the particle, represents the surface of the particle, represents the angular acceleration, and J represents the moment of inertia, represents the resultant moment, represents the spatial coordinates of a point on the surface of the particle, represents the spatial coordinates of the center of gravity of the solid particle, P represents pressure, and I represents the unit tensor, represents the viscous stress tensor, represents the velocity of the solid particle at time t, represents the velocity of the solid particle at time t, represents the acceleration of the solid particle at time t, represents the acceleration of the solid particle at time t, represents the time step, represents the spatial coordinates of the solid particle at time t, represents the of the solid particle at time t, represents the angular velocity of the solid particle at time t, represents a rotational angular velocity of the solid particle at the moment, denotes a rotational angular acceleration of the solid particle at the moment, denotes a rotational angular acceleration of the solid particle at the moment, denotes a rotational angle of the solid particle at the moment, denotes a rotational angle of the solid particle at the moment.

[0115] The present application directly aims at the complexity (viscoelasticity) of material constitutive, solves the precision and efficiency problems existing in the current fluid-structure coupling algorithm, can not only improve the safety and reliability of engineering structure, but also provide theoretical support for the development of new materials and new technologies.

Claims

1. A computational method for a stabilized viscoelastic-fluid-structure interaction problem based on overlapping meshes, characterized in that, include: S1: Process the flow field region and solid particles to obtain overlapping background mesh and sub-mesh; S2: Based on the background grid and subgrid, the interpolation boundary, contribution unit, and interpolation unit are obtained through analysis; S3: Based on interpolation boundaries, contribution elements, and interpolation elements, for viscoelastic fluid-structure interaction problems, the flow field of the background mesh is calculated using the finite element method and stabilization method to obtain the flow field information of the background mesh; S4: Using bicubic spline interpolation, the flow field information of the background grid is transferred to the sub-grid to obtain the flow field information of the sub-grid; S5: Using the flow field information of the subgrid, the updated flow field information of the subgrid is obtained by calculation using the finite element method and the stabilization method; S6: Use bilinear Lagrange interpolation to transfer the updated subgrid flow field information to the background grid to obtain the updated background grid flow field information; S7: Based on the updated background grid flow field information, parameters are iteratively updated through calculation to obtain simulation data of the motion trajectory and rotation state of solid particles in the flow field region, thus completing the calculation of the stabilized viscoelastic fluid-structure interaction problem.

2. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, S1 includes: The flow field region is divided using the grid spacing to obtain the background grid; Sub-mesh is generated based on the geometry and spatial position of solid particles using sub-mesh generation parameters.

3. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, S2 includes: Based on the background mesh and sub-mesh, the outer boundary of the sub-mesh is directly defined as the interpolation boundary. On the interpolation boundary, the cells of the sub-mesh serve as interpolation cells to receive information from the background mesh, while the cells of the background mesh serve as contributing cells to provide flow field data. The contributing cells determine whether a node is inside the object by calculating the distance from each node of the background mesh to the center of the object: if the distance from the node to the center of the object is less than or equal to the radius of the object, the background mesh cell containing that node is considered to be inside the object and is discarded, thus providing flow field data.

4. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, The background grid flow field information is determined by a velocity correction equation, wherein the spatial discrete format expression of the velocity correction equation is: ； in, Represents the mass matrix, Indicates the first The velocity field updated at each time step This represents the intermediate velocity encountered during the solution process. Indicates the time step. Represents the Reynolds number. Represents the pressure gradient matrix. Indicates the first The pressure field after each time step update This represents the pressure field at the nth time step.

5. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, The flow field information of the sub-grid is determined by an interpolation formula, wherein the expression of the interpolation formula is: ； ； in, This represents the flow field information at the sub-mesh boundary nodes. Represents the x-axis coordinate of the node. Represents the y-axis coordinate of the node. This represents the flow field information of the node at the i-th row and j-th column in the contributing unit corresponding to this node. Let represent the improved form of the i-th uniform B-spline basis function. This represents the x-coordinate of the bottom-left corner node of the cell in the background mesh. This represents the y-coordinate of the bottom-left node of the cell in the background mesh. Indicates the grid size of the background grid.

6. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, The updated subgrid flow field information is determined by the fluid control equations within the subgrid, wherein the expression for the fluid control equations within the subgrid is: ； in, Indicates speed, Indicates time, Indicates the velocity of solid particles. Indicates the pressure gradient. Indicates stress, Indicates viscosity ratio.

7. The calculation method for stabilized viscoelastic-fluid-structure interaction problems based on overlapping meshes according to claim 1, characterized in that, The expression for the bilinear Lagrange interpolation is: ； ； ； ； in, This represents the updated flow field information of the nodes in the background mesh that need to have their flow field information updated via interpolation. This represents the x-axis coordinate of the above nodes. This represents the y-axis coordinate of the above nodes. Represents the Lagrange interpolation coefficients. This represents the flow field information of each node in the contributing unit corresponding to the above node, corresponding to u. This represents the x-axis coordinate of the i-th node in the contributing unit. This represents the y-coordinate of the i-th node in the contributing unit. This represents the Lagrange interpolation coefficient corresponding to the i-th node in the contributing unit.

8. The calculation method for the stabilized viscoelastic-fluid-structure interaction problem based on overlapping meshes according to claim 1, characterized in that, The expression for the simulated data is: ； ； ； ； ； ； ； in, Indicates quality, Indicates acceleration. Represents gravity. This represents the fluid force received by a solid particle. The stress tensor represents the stress tensor generated by a fluid at a point on the particle surface. Indicates solid boundary, This represents the outward normal vector at a point on the particle surface. Indicates the surface of particles, J represents the rotational angular acceleration, and J represents the moment of inertia. Indicates the resultant torque. This represents the spatial coordinates of a point on the particle surface. The coordinates of the centroid of the solid particle are represented by P, where P represents pressure and I represents the unit tensor. Represents the viscous stress tensor. express The velocity of the solid particles at any given moment This represents the velocity of the solid particle at time t. express The acceleration of the solid particles at any given moment. This represents the acceleration of the solid particle at time t. Indicates the time step. express The spatial coordinates of the solid particles at any given time. express Solid particles at all times, express The angular velocity of the solid particle at any given moment. express The angular velocity of the solid particle at any given moment. express The rotational angular acceleration of the solid particle at any given time. express The rotational angular acceleration of the solid particle at any given time. express The rotation angle of the solid particles at any given time. express The rotation angle of the solid particle at any given time.

Images