Topological optimization method, device and equipment for material distribution based on particle flow

Abstract

The invention discloses a topological optimization method, device and equipment for material distribution based on particle flow, and relates to the technical field of material structure optimization. According to the method, isogeometric analysis is introduced into the material point method in topological optimization of material distribution, the isogeometric material point method is provided, consistency between geometric modeling and numerical analysis is achieved through a unified primary function of isogeometric analysis, benefit from high-order continuity of the primary function of isogeometric analysis, and the method has the advantages that the method is simple and convenient to implement. The introduction of isogeometric analysis can effectively eliminate and suppress numerical artifacts in topological optimization and unit cross instability in a traditional substance point method. According to the isogeometric material point method, movable particles of material distribution in a physical domain serve as design variables, a set of particle-grid physical information and sensitivity transmission model is established in combination with Lagrange-Euler mixed description, corresponding accurate geometric topology is output while particle evolution is achieved, and the method has the advantages of being high in accuracy and high in reliability. And the accuracy and the flexibility of topological optimization of the material distribution structure are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of material structure optimization, and particularly relates to a topology optimization method, device, and equipment for material distribution based on particle flow. Background Art

[0002] Currently, as an advanced design tool, topology optimization has shown its great potential in applications in multiple fields in recent years, especially in the fields of aerospace, automotive, civil engineering, biomedical engineering, and material design. The goal of topology optimization is to find the optimal material distribution with expected structural performance in the design space. With the rapid development of additive manufacturing and next-generation computing infrastructure, the application of topology optimization has also gradually expanded and become a key tool for solving complex engineering problems.

[0003] In topology optimization, how to effectively establish a material description model to express the evolution of the structural topology within the design domain has become a key issue affecting the optimization effect and computational efficiency. Currently, the description of the structural physical field in topology optimization is mainly divided into two types of methods: Eulerian description and Lagrangian description. The Eulerian description method defines the topological domain on a fixed grid or nodes and tracks the structural changes during the optimization process. The variable density method, the Solid Isotropic Material with Penalization (SIMP) method, the Evolutionary Structural Optimization (ESO) method, etc. are widely applied to optimization problems. However, the Eulerian description method has certain limitations in capturing complex structural boundaries and adapting to significant shape and topology changes, which may lead to a decrease in the numerical accuracy of the optimization results. To address these limitations, the Lagrangian method, as a method based on dynamic grids or particle motion, provides a more flexible way of topological modeling. The Lagrangian method is applicable to large deformation and fracture scenarios and provides new ideas for topology optimization through methods such as the Element-Free Galerkin (EFG), Reproducing Kernel (RK), Meshless Local Petrov-Galerkin Method (MLPG), and Smoothed Particle Hydrodynamics (SPH). However, isolated mass nodes may appear during the optimization process in the Lagrangian method, and these nodes are disconnected from the main load-bearing structure, and the computational cost is relatively high. To overcome the respective defects of the Eulerian method and the Lagrangian method, researchers have proposed various hybrid Lagrangian-Eulerian methods, such as the Particle-in-Cell (PIC) method and the Fluid-Implicit Particle (FLIP) method. Among these methods, the Material Point Method (MPM), as a Lagrangian-Eulerian hybrid method, has been widely applied in fields such as computer graphics and mechanical engineering. The MPM discretizes the design domain on a fixed Eulerian grid and solves the momentum balance equation, which has the advantages of improving computational efficiency and scalability.

[0004] In recent years, the Lagrangian-Eulerian Topology Optimization (LETO) method combines the MPM with fixed Eulerian grid points and uses a smooth radial SPH kernel to achieve information transfer between Lagrangian particles and Eulerian grids, thus effectively combining the flexibility of the Lagrangian method with the computational efficiency of the Eulerian method.

[0005] However, the LETO method still faces numerical defects brought about by the classical Finite Element Method (FEM), such as inconsistencies between Computer Aided Design (CAD) models and Computer Aided Engineering (CAE) models, and low-order continuity of basis functions of adjacent elements. To a certain extent, these problems lead to cell-cross errors in the Material Point Method, thereby affecting the accuracy of optimization results, causing connectivity problems and isolated elements in the optimization results, and reducing the topological optimization accuracy and flexibility of material distribution. Summary of the Invention

[0006] Based on this, it is necessary to provide a topological optimization method, device, and equipment for material distribution based on particle flow to address the above technical problems.

[0007] The present invention adopts the following technical solutions: The present invention provides a topological optimization method for material distribution based on particle flow, including: Initializing a particle design model according to design parameters to initialize the physical information corresponding to the design variables of particles in the physical space based on isogeometric analysis meshes; Mapping the physical information of particles to the control points of isogeometric analysis meshes through a smoothed radial SPH kernel, and performing numerical projection through a control density projection function to obtain the physical density of the control points; Constructing a density distribution function based on the physical density of the control points to represent the material distribution within the design domain, and calculating the density values at Gauss points; Determining the sensitivities of the objective function and constraint function of topological optimization to Gauss points, and based on the sensitivity chain derivative and the density values at Gauss points, transferring the sensitivities of the objective function and constraint function to Gauss points to the control points, and then transferring them to the design variables of the particles according to the sensitivity chain derivative to determine the sensitivities of the objective function and constraint function to the design variables of the particles; Iteratively updating the design variables of the particles by the method of moving asymptotes according to the sensitivities of the objective function and constraint function to the design variables of the particles to achieve particle flow topological optimization of the material distribution.

[0008] The present invention provides a topological optimization device for material distribution based on particle flow, including: A particle model establishment module for initializing a particle design model according to design parameters to initialize the physical information corresponding to the design variables of particles in the physical space based on isogeometric analysis meshes; A physical information mapping module, configured to map the physical information of particles to the control points of the isogeometric analysis grid through a smooth radial SPH kernel, and perform numerical projection through a control density projection function to obtain the physical density of the control points; A physical density analysis module, configured to represent the material distribution within the design domain through a density distribution function constructed based on the physical density of the control points, and calculate the density values at the Gauss points; A sensitivity analysis transfer module, configured to determine the sensitivities of the objective function and constraint function of the topology optimization with respect to the Gauss points, and based on the sensitivity chain derivative and the density values of the Gauss points, transfer the sensitivities of the objective function and constraint function with respect to the Gauss points to the control points, and then transfer them to the design variables of the particles according to the sensitivity chain derivative, to determine the sensitivities of the objective function and constraint function with respect to the design variables of the particles; A topology optimization module, configured to iteratively update the design variables of the particles by the moving asymptotes method according to the sensitivities of the objective function and constraint function with respect to the design variables of the particles, to achieve the particle flow topology optimization of the material distribution.

[0009] The present invention provides a computer-readable storage medium, where the storage medium stores a computer program, and when the computer program is executed by a processor, it implements the above-mentioned topology optimization method for material distribution based on particle flow.

[0010] The present invention provides a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, it implements the above-mentioned topology optimization method for material distribution based on particle flow.

[0011] The above-mentioned at least one technical solution adopted by the present invention can achieve the following beneficial effects: In the topology optimization of material distribution, the present invention introduces isogeometric analysis into the material point method, and proposes an isogeometric material point method. The isogeometric material point method of the present invention uses the movable particles of the material distribution in the physical domain as design variables, combines the Lagrangian-Eulerian mixed description, establishes a set of particle-grid physical information and sensitivity transfer models, constructs a set of particle flow topology optimization systems, and ensures the consistency of the particle-topology-analysis model. Through the unified basis function of isogeometric analysis, consistency is achieved between geometric modeling and numerical analysis. Thanks to the high-order continuity of the isogeometric analysis basis function, the introduction of isogeometric analysis can effectively eliminate and suppress numerical artifacts in topology optimization and the element crossing instability phenomenon in the traditional material point method. The isogeometric material point method of the present invention can output the corresponding accurate geometric topology while realizing particle evolution, significantly improving the accuracy and flexibility of the topology optimization of the material distribution structure. Description of the Drawings

[0012] The accompanying drawings described herein are used to provide a further understanding of the present invention, and constitute a part of the present invention. The schematic embodiments of the present invention and their descriptions are used to explain the present invention, and do not constitute an improper limitation of the present invention. In the drawings: Figure 1 It is a schematic flow chart of a topology optimization method for material distribution based on particle flow provided by the present invention; Figure 2 It is a schematic flow chart of a particle flow topology optimization framework under a Lagrangian-Euler hybrid system provided by the present invention; Figure 3 It is a schematic diagram of the mapping relationship among particles, control points and Gaussian points in the design domain provided by the present invention; Figure 4 It is a schematic diagram of three different density projection functions provided by the present invention. (a1) and (a2) correspond to the classical density projection function; (b1) and (b2) correspond to the density projection function in the negative derivative interval; (c1) and (c2) correspond to the smooth transition curve projection function based on the improved Heaviside function; Figure 5 It is a schematic diagram of the optimization results of three different density projection functions provided by the present invention; Figure 6 It is a schematic diagram of the optimization process and result 1 under an extreme initial distribution provided by the present invention; Figure 7 It is a schematic diagram of the optimization process and result 2 under an extreme initial distribution provided by the present invention; Figure 8 It is a schematic diagram of the optimization process and result of a 1 / 4 circular ring provided by the present invention; Figure 9 It is a schematic diagram of the optimization process and result of a three-dimensional michel structure provided by the present invention; Figure 10 It is a schematic diagram of a topology optimization device for material distribution based on particle flow provided by the present invention. Detailed implementation manners

[0013] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below in conjunction with specific embodiments of the present invention and the corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0014] Currently, in the topological optimization of material distribution, there are still some technical problems in the existing technical solutions. The Eulerian method has limitations in capturing complex geometries and mesh adaptability, and it is difficult to handle the accuracy problems brought by topological structure changes. Although the Lagrangian method is flexible, it has a high computational cost and may generate isolated mass nodes during the optimization process, affecting the reliability of the results. The hybrid Lagrangian-Eulerian method (such as the LETO method) combines the advantages of both, but still faces problems such as numerical accuracy, stability, and numerical artifacts. Isogeometric analysis (IGA) has high accuracy in topological optimization, but the combination with the Lagrangian-Eulerian method lacks research and its practical application is relatively limited. In addition, the numerical instability of topological optimization and the need for post-processing steps also make practical applications complicated, especially in multi-disciplinary design problems, how to balance the requirements of different disciplines. Future research needs to improve these methods to solve problems such as numerical accuracy, computational efficiency, stability, and post-processing, in order to improve the practical application effect and applicability of topological optimization.

[0015] The following will, in conjunction with the accompanying drawings, elaborate on the technical solutions provided by each embodiment of the present invention.

[0016] Figure 1 It is a schematic flow diagram of a topological optimization method for material distribution based on particle flow in the present invention, specifically including the following steps: S101: Initialize the particle design model according to the design parameters, and initialize the physical information corresponding to the design variables of the particles in the physical space based on the isogeometric analysis mesh.

[0017] S102: Map the physical information of the particles to the control points of the isogeometric analysis mesh through a smooth radial SPH kernel, and perform numerical projection through the control density projection function to obtain the physical density of the control points.

[0018] S103: Construct a density distribution function according to the physical density of the control points to represent the material distribution in the design domain, and calculate the density values at the Gauss points.

[0019] S104: Determine the sensitivities of the objective function and constraint function of the topological optimization to the Gauss points, and according to the sensitivity chain derivative and the density values at the Gauss points, transfer the sensitivities of the objective function and constraint function to the Gauss points to the control points, and then transfer them to the design variables of the particles according to the sensitivity chain derivative, and determine the sensitivities of the objective function and constraint function to the design variables of the particles.

[0020] S105: According to the sensitivities of the objective function and constraint function to the design variables of the particles, use the moving asymptotes method to iteratively update the design variables of the particles to achieve the particle flow topological optimization of the material distribution.

[0021] The server mentioned in the present invention can be a server set up in a service platform, or a device such as a desktop computer or a laptop computer that can execute the solution of the present invention. For the convenience of description, only the server will be used as the execution entity for description below.

[0022] Figure 2 It is a schematic flow diagram of a particle flow topology optimization framework under a Lagrangian-Euler hybrid system in the present invention. Figure 2 The topology optimization process is illustrated in the figure, and the following will specifically describe each part thereof.

[0023] First, the server can initialize the particle design model with design parameters, uniformly configure the physical coordinates of all particles in the parameter space based on the IGA background grid, and inversely solve the initial particle material consumption according to the volume constraint, and map the particles to the physical space using NURBS basis functions to obtain the initial design domain.

[0024] Specifically, in one or more embodiments of the present invention, the present invention is to construct an IGA design domain based on NURBS basis functions in the physical space. As is well known, NURBS is constructed by B-splines, so the definition of B-splines will be discussed first. In a one-dimensional parameter space, a set of non-decreasing coordinates are divided and represented as . is the knot, is the order of the polynomial, is the number of basis functions used to construct the B-spline curve. The parameter space is divided into a series of sub-interval elements by the knots. The B-spline basis functions are recursively defined by the Cox-de Boor formula, starting from a piecewise constant function . At this time, .

[0025] For the basis functions are recursively defined as follows: .

[0026] This is the Cox–de Boor recurrence formula. It can be seen that the basis function with order p = 0 is the same as the standard piecewise constant basis function, and the basis function with order p = 1 is the same as the standard linear finite element basis function. However, the basis functions with higher orders ( p ≥2) are completely different.

[0027] Although B-splines have great flexibility in modeling various geometries, some important shapes, such as circles and ellipsoids, cannot be accurately represented. Therefore, the NURBS basis functions are obtained by a projective transformation of the B-spline basis functions, and the specific form is: 。

[0028] Among them, the weight function is defined as: 。

[0029] Among them, is the weight of each basis function. The NURBS curve is constructed by the linear combination of NURBS basis functions, and the expression is: 。

[0030] Among them, are the control points. Given the additional knot vectors and , the definition of NURBS basis functions in two dimensions is as follows: 。

[0031] Similarly, for three dimensions, there is: 。

[0032] Among them, and are the B-spline basis functions in the other two parameter directions, and are the polynomial orders in the two parameter directions respectively, and are the corresponding numbers of control points. The NURBS basis function can inherit all the favorable characteristics of the B-spline basis function.

[0033] In the parameter space - physical space mapping, the single-index notation is used, and the index is used to label the NURBS basis function. The geometric mapping in the setting can be expressed as: 。

[0034] In the formula, is the parameter coordinate corresponding to the number of control points in the influence domain, is the th control point. represents the physical geometry. For the sake of illustration, assume that the IGA solution in the parameter domain is a scalar value and is controlled by the same NURBS basis function, and can be expressed as: 。

[0035] Among them, represents variables, such as the degrees of freedom (DOFs) of the control points. The IGA solution in the physical domain is expressed as the advancement of its parameter corresponding solution. The implementation of this transformation's geometric mapping can be expressed as: 。

[0036] Among them, 。

[0037] After obtaining the coordinate distribution of particles in the physical space and the material consumption, the preliminary setting of the design domain is completed at this time. It is necessary to inversely calculate the average material consumption of the particles according to the volume constraint in the design parameters. The ultimate goal of this step is to achieve a relative error of less than 0.1% between the initial volume fraction and the set volume constraint.

[0038] Then, the server can perform particle-control point information transfer (P2C), transfer the particle physical information to the control points using the SPH kernel, and obtain the physical density of the control points through further projection.

[0039] Specifically, in one or more embodiments of the present invention, the present invention uses a cubic B-spline basis function as the kernel function to complete the particle-control point mapping. Due to the gradient discontinuity of the linear basis function at the mesh boundary in the traditional finite element method (FEM), when a particle crosses the element boundary, an "element crossing error" will occur. This error will cause sudden changes in the internal force and velocity gradient, thereby causing sensitivity deviation in the topology optimization process and further reducing the numerical accuracy of the result. In addition, in IGA, if the inverse mapping of the NURBS basis function is to be realized, the material points must be inversely mapped from the physical space to the parametric space. This transformation involves solving complex non-linear equations and it is difficult to effectively ensure the accuracy and efficiency of the calculation. To address the above two problems, the present invention proposes a smooth radial SPH kernel, using a cubic B-spline kernel function to replace the original basis function, thereby reducing the "element crossing error". After mapping the particle physical information to the IGA mesh, the control equation is solved on the mesh. Due to the high-order continuity of the NURBS basis function, the updated control point information and the sensitivity information mapped to the particles will be smoother.

[0040] In the Material Point Method (MPM), the coordinate basis is determined by the position vector of the particle in the initial undeformed structure state, which is called the material coordinate or Lagrangian coordinate. The position of the same particle in the deformed configuration is represented by the spatial coordinate or Euler coordinate . By using a function to describe the motion and deformation of the solid, the relationship between the spatial coordinate and the material coordinate can be established as follows: .

[0041] The material point 's velocity is represented by and is defined as the rate of change of the position of the material point: .

[0042] This is the Lagrangian velocity field. Since MPM uses the Lagrangian description, the Euler form of this velocity field is not discussed. The acceleration of the material point , which is the derivative of the velocity with respect to the material time, is expressed as: .

[0043] Consider the general non - linear elasticity problem. For simplicity, neglecting the source terms and boundary conditions, the momentum balance problem can be expressed in the strong form in Eulerian description as follows: .

[0044] where is the mass density, is the velocity of the material point. is the material time derivative, is the Cauchy stress, is the time - varying region deformed with the velocity . Meanwhile, the law of conservation of mass is described by the mass - conservation equation, i.e., the continuity equation, which has the form: . Considering the mass conservation in the above equation, for the intermediate configuration , expressing it at the fixed material point gives: .

[0045] For simplicity, the time parameter of the function is represented by a subscript, e.g., . is the deformed position of the material point at time instant is the displacement from to , is the unit tensor, represents the local volume change. Following the convention of using an over - bar to denote the Lagrangian containing the reference configuration . Next, the concept of the material time derivative is introduced. Consider a field (which can be a scalar, vector, or tensor) defined throughout the body. The change of this field at a specific material point is determined by calculating the material time derivative of Lagrangian description: In the Lagrangian description, the two independent variables, the material coordinates X and time t, can be directly defined as the partial derivative of the given field with respect to time, expressed as: .

[0046] Eulerian description: In this description, the field is defined as . At this time, not only the change in time but also the change in the particle's spatial position X has to be considered. Therefore, the partial derivative of the material description of the given field with respect to time has to be calculated while keeping X fixed. Using the chain rule, the formula for the material time derivative of the Eulerian scalar field is obtained: .

[0047] According to the above definition, the weak form of the equation can be expressed as: find , such that for any there is: .

[0048] where is the incremental velocity. Let the initial condition at time be: and , where is the Euler velocity field. The initial data satisfies: .

[0049] In the formula, is the Lagrangian velocity field or material velocity field, which maps the material points in the domain to their corresponding velocities at time t. This velocity field satisfies the following evolution equation: .

[0050] Time domain discretization: Define , and use single-step finite differences to discretize the interval . To clarify the idea and establish the notation convention, describe a simple explicit Euler integration method: .

[0051] where satisfies .

[0052] The deformation evolves according to the following discrete form: .

[0053] where . Using the finite difference to approximate the acceleration, according to the time-discrete problem can be obtained: where, The elimination of indicates the role of as an approximation of the material time derivative. The n-level stress is determined by a material-specific constitutive model. The determinant is given by det , and the evolution of is as follows:

[0054] Using the finite difference to approximate the velocity field evolution equation, that is .

[0055] This can ensure that in the case of , the time-discrete Lagrangian velocity remains unchanged for each material point throughout the motion.

[0056] The standard MPM discretizes such semi - discrete problems by selecting a finite - dimensional subspace where and are discrete representations in space, and for material points (particles) within the domain, a finite quadrature rule is used for integration: .

[0057] In the formula represents the particle, is the weight of the particle or the particle volume.

[0058] In this framework, the control - point information is obtained by weighted summation of the physical quantities of particles in the neighboring domain , and the specific SPH kernel calculation method is as follows: .

[0059] where is the position of the control point, is the position of the particle within the neighborhood of the control point, represents the material consumption of the th particle, is the kernel function, is the kernel size, is the volume of each control point, which depends on the average proportion of each control point in the element. In this paper, the cubic - spline kernel function is adopted: .

[0060] where is a constant, and its values in 2D and 3D are and respectively.

[0061] The aggregated variable at the control point is initially obtained through SPH kernel mapping. At this time, the physical density of the control point is evaluated from the aggregated variable of the control point through a numerical projection scheme. The mathematical formula of the aggregated - variable projection scheme is expressed as: .

[0062] where and represent the coefficients of two functions. Represents the interval points of a piecewise function. The first segment of the numerical projection function can be regarded as a smooth variant of the Heaviside function, and the latter segment is a cubic function characterized by a negative derivative. The latter half of the projection function in the above formula has an important impact on the effectiveness of the projection scheme, especially in realizing a structural topology with smooth and continuous boundaries. This is mainly because the negative derivative can not only maintain the smoothness of the sensitivity information but also appropriately correct the excessive aggregation of particles at adjacent control points.

[0063] The above process is called the "particle to control point (P2C)" stage. This stage provides control information for the subsequent Eulerian grid in I-MPM. Subsequently, the control point information is updated on the grid and the information is retransmitted to the particles through NURBS mapping.

[0064] After that, the server can perform IGA static solution, including: control point - Gauss point information transfer (C2G), representing the material distribution in the design domain using the Density Distribution Function (DDF), and calculating the density values at Gauss points; solving the control equations based on the IGA background grid to obtain the sensitivities of the objective function and constraint function to Gauss points; Gauss point - control point information transfer (G2C), according to the chain derivative of sensitivity, using NURBS basis functions to transfer the Gauss point sensitivity information to control points. Figure 3 It is a schematic diagram of the mapping relationship among particles, control points, and Gauss points in a design domain of the present invention.

[0065] Specifically, in one or more embodiments of the present invention, the physical density of the control point, i.e., the control density can be linearly combined with NURBS basis functions to form a Density Distribution Function (DDF), denoted as , which is used to represent the structural topology of the entire design domain. The specific form is as follows: .

[0066] In the formula is the physical density of the control point, corresponding to a strict limit range , is the minimum value of the control density, which is used to eliminate numerical singularities in the optimization process. , and They are the total number of control points in the directions of relevant parameters respectively. During the optimization process, the density distribution function (DDF) can be regarded as a high-order density response surface. Therefore, the contour lines of the DDF are used to represent the structural boundary. The DDF values above the contour lines represent the entities in the structural domain, and the DDF values below the contour lines represent the voids in the design domain. In fact, the entire representation mechanism is directly inherited from the implicit description model in the classical level set method, and the differences have also been fully studied in previous works on isogeometric topology optimization. Moreover, NURBS is also used to construct the DDF. The detailed mechanism of the DDF representing the structural topology is given by the following formula: 。

[0067] In the formula is the value of the contour line (2D) or the isosurface (3D) of the DDF, which is set to 0.5 in the subsequent numerical examples. In the present invention, the DDF is constructed from the physical densities of the control points, and the physical densities of the control points are determined by the corresponding particle design variables. The density distribution function can be expressed in more detail as: 。

[0068] That is, the material distribution within the entire design domain. Calculate the density value of the Gauss point, that is: 。

[0069] NURBS has several important properties: (1) strong convex hull property (2) differentiability (3) local modification property (4) variation diminishing property. These properties are crucial for the construction of the Gauss density. This construction form can perform local representation, and it depends on the local support domain of the current node and there is 。

[0070] When solving the Gauss point density, the properties of the NURBS basis functions indicate that they are affected by the surrounding control points, and the number of them is where is the order of NURBS, represents the dimension of NURBS. The final Gauss density is achieved through the linear combination of the NURBS basis functions and the control density.

[0071] Next, the displacement field satisfying the Dirichlet boundary conditions can be solved on the IGA mesh is the known displacement vector on the Dirichlet boundary . In the present invention, for linearly elastic materials, only a single linear system needs to be solved.

[0072] In one or more embodiments of the present invention, taking the classical optimization design problem of maximizing the structural load-bearing capacity as an example, the effectiveness and advantages of the Lagrangian-Eulerian hybrid topology optimization description are studied. In the proposed particle flow topology optimization method, Lagrangian description particles are used as design variables, and the evolution process of the structural topology is derived, which can effectively optimize all physical information such as the position of particles and material consumption, and seek a reasonable material distribution with expected structural performance. In the Euler-based description grid, the isogeometric material point method is used to solve the linear elastic equilibrium equation, which ensures the high-order continuity of NURBS basis functions and the consistency of geometric models with high numerical accuracy, and successfully solves the key problem of the jump of particle stresses and other historical variables in adjacent elements. The general formulation of the above optimization form is to search for the minimum structural compliance of the displacement field by iteratively optimizing the design variables of the particles to search for the displacement field of the minimum structural compliance, where the interaction of the physical information of particles, control points, and Gauss points is described by the Lagrangian-Eulerian particle flow model, and its mathematical equation is: .

[0073] Where is the set of design variables, which consists of the positions of all particles and the material consumption of the particles within the influence domain. represents the design domain and is the total number of particles in the design domain. is the optimization objective characterizing the structural compliance, which can reflect the load-bearing capacity of the structure. represents the displacement field of the structure in the design domain, which is calculated by the linear elastic equivalent equation. is the volume constraint, where is the maximum value of the material volume constraint. is the prescribed displacement vector on the Dirichlet boundary,

[0074] is the applied load on the Neumann boundary. The sensitivities of the objective function and constraint function to the Gauss points can be calculated using the obtained displacement field .

[0075] Based on this, the control point-particle information transfer (C2P) can be performed. Based on the KD-Tree spatial neighborhood search to accelerate the search for control points in the particle influence domain, the inverse operation of the SPH kernel can be efficiently performed, and the sensitivities of the physical coordinates and material consumption design variables of the particles can be solved based on the sensitivity chain derivative.

[0076] Specifically, in the optimization design framework of the present invention, the physical information of particles is the design variable for optimizing the structural topology. To calculate the sensitivities of the objective function and the constraint function with respect to the design variables, that is and , an inverse projection calculation needs to be introduced in step S102, where the sensitivity information of the objective function and the constraint function with respect to the control points will be passed to the particles (C2P) using the SPH kernel function and the numerical projection scheme. The specific derivation is as follows:

[0077] In the linear elastic system, the sensitivity analysis of the objective function adopts the exact form obtained by the adjoint method in traditional linear topology optimization. Taking the physical information of the particles , that is, the position and the material consumption as the design variables, the evolution of the structural topology is derived. Using the physical information that controls the density, the corresponding sensitivity analysis of the objective function with respect to the position and the material consumption is calculated using the chain rule, expressed as: .

[0078] In the formula, the parameter is the physical density of the control point, and its derivative can be calculated by the first derivative of the SPH kernel. The above equation can be expressed in a more detailed multi-variable form as follows: .

[0079] Among them, represents the aggregation variable of the control point, which is defined by the SPH kernel function and transfers the physical information of the particles to the control point. The derivative of can be derived through the following equation, specifically expressed as:

[0080] Among them, and represent the material consumption and the position of the particle respectively. The derivative of the aggregation variable of the control point with respect to the physical variables of the particle contains two components. The aggregation variable of the control point is obtained by weighted averaging the particles within the influence domain through the SPH kernel function. Similarly, the sensitivities of these particles are also calculated according to this equation. First, find the partial derivative of the particle material consumption , and the first term of the design variable sensitivity can be calculated through the following equation: .

[0081] Secondly, the derivative of the relative coordinate offset of the particle needs to be calculated. Let represent the relative coordinates between the control points within the current particle influence domain, and represent the relative distance between the particle and the control point within the current influence domain. The results are as follows: .

[0082] Where is the unit direction vector between the current particle and the control point. Using the aggregation variable information of the control points, a numerical projection scheme is adopted to evaluate the physical density of the control points, so that the sensitivity analysis of the control density to the aggregation variables of the control points can be derived through the numerical projection scheme, as follows: .

[0083] To establish a connection between the objective function or volume constraint and the control density, it is necessary to introduce the Gaussian points for calculating the control equations as intermediate variables in the chain rule, as follows: .

[0084] The first-order sensitivity of the DDF to the control density can be derived as: .

[0085] The sensitivity analysis of the corresponding objective function and constraint function to the control density can be calculated as: .

[0086] Where represents the physical density at the th control point. The final form of the first-order derivative of the objective function with respect to the design variables (i.e., the physical position and material consumption of the particles) can be described as: .

[0087] In the formula, is the strain tensor, is the elastic tensor matrix, is the material consumption of the particles within the influence domain, is the penalty term.

[0088] In the above formula, it is not difficult to find that the adjoint formula of the design sensitivity uses the position and density information of the particles, as well as the density information of the Gaussian integration points related to the control points in the SPH kernel influence domain for calculation. In optimization, the physical information of the particles, that is, the position and material consumption of the particles, plays an important role in promoting the structure topology to seek the optimal material layout. Similarly, the sensitivity analysis of the constraint function to the two physical variables of the particles can also be calculated as: .

[0089] Finally, the server can update the particle design variables based on the Method of Moving Asymptotes (MMA) optimizer, drive the evolution of the material consumption and physical coordinates of the particles, and form the structure topology and performance optimization characterized by the particle flow.

[0090] Specifically, the Moving Asymptote Method (MMA) is used to solve the optimization problem that jointly optimizes the positions of particles and material consumption. This optimizer is designed for general structural optimization problems with inequality constraints and box constraints. The algorithm approximates the original problem through a series of separable convex optimizations. In each iteration, it sets two asymptotes for each variable to constrain the search interval. These asymptotes will be updated according to each sub-optimal value.

[0091] To obtain high-quality results from MMA, careful parameter tuning is required, and different examples may have different optimal sets of MMA parameters. There are three parameters to be tuned: asyinit; asyincr, and asydeck. In this paper, these parameters are set to 0.02, 1.05, and 0.65 respectively, the same as the original MMC method. These parameters are used in all numerical examples. To further stabilize and accelerate the optimization in a consistent manner, the following regularization is also applied in MMA:

[0092] 1. Control the step size of each variable by modifying its box constraint at each iteration. The single-step change in particle density is controlled to be below 0.5, and the position change in each dimension is also controlled to be below twice the background grid spacing. This step size control can stabilize the optimization.

[0093] 2. Following the MMC method, scale the objective gradient and volume constraint so that their norms are both 1. In addition, the objective and volume constraints are scaled accordingly to ensure that the scaled gradients are consistent. Scaling significantly accelerates the optimization.

[0094] Based on Figure 1 the topology optimization method of particle flow-based material distribution as shown, in the topology optimization of material distribution, the present invention introduces isogeometric analysis into the material point method, and proposes an isogeometric material point method. The isogeometric material point method of the present invention takes the movable particles of material distribution in the physical domain as design variables, combines the Lagrangian-Eulerian mixed description, establishes a set of particle-grid physical information and sensitivity transfer models, and constructs a set of particle flow topology optimization systems to ensure the consistency of the particle-topology-analysis model. Through the unified basis function of isogeometric analysis, the consistency between geometric modeling and numerical analysis is achieved. Benefiting from the high-order continuity of the isogeometric analysis basis function, the introduction of isogeometric analysis can effectively eliminate and suppress numerical artifacts in topology optimization and the element crossing instability phenomenon in the traditional material point method. The isogeometric material point method of the present invention can output the corresponding precise geometric topology while realizing particle evolution, significantly improving the accuracy and flexibility of the topology optimization of material distribution structures.

[0095] (1) The Lagrangian-Eulerian hybrid description particle flow topology optimization method provided by the present invention aims at the problems existing in the traditional finite element-based material point method and material point method topology optimization, such as element crossing errors, insufficient accuracy of topological field description, and problems such as QR mode, connectivity problems, and isolated elements in structures. An I-MPM analysis framework is proposed and integrated with topology optimization to form a particle flow topology optimization (PFTO) model. Numerical discussions clearly reveal the effectiveness and indispensability of I-MPM in solving the state equation with higher numerical accuracy, as well as the advantages of PFTO topology results in terms of description accuracy, solution efficiency, etc.

[0096] (2) The Lagrangian-Eulerian hybrid description particle flow topology optimization method provided by the present invention proposes a smooth control density projection curve with a negative derivative interval constructed based on the Heaviside function, aiming to eliminate the phenomenon of excessive particle aggregation existing in the traditional sensitivity-driven material point method topology optimization, which usually leads to numerical problems such as boundary waviness of the structural topology. By introducing the segmented interval points with negative derivatives in the projection function and applying reverse sensitivity to the particles in the over-aggregation region during the optimization process, it can effectively avoid the infinite increase of the control point aggregation variables in the optimization process and greatly reduce the numerical problems in the optimization solution.

[0097] Figure 4 It is a schematic diagram of three different density projection functions in the present invention. Among them, (a1) and (a2) correspond to the classical density projection function; (b1) and (b2) correspond to the density projection function with a negative derivative interval; (c1) and (c2) correspond to the smooth transition curve projection function based on the improved Heaviside function.

[0098] Figure 5 It is a schematic diagram of the optimization results of three different density projection functions in the present invention. (a1) corresponds to the density representation diagram of particles under the classical density projection function, and (a2) corresponds to the structure representation diagram obtained by optimization under the classical density projection function; (b1) corresponds to the density representation diagram of particles under the density projection function with a negative derivative interval, and (b2) corresponds to the structure representation diagram obtained by optimization under the density projection function with a negative derivative interval; (c1) corresponds to the density representation diagram of particles under the smooth transition curve projection function based on the improved Heaviside function, and (c2) corresponds to the structure representation diagram obtained by optimization under the smooth transition curve projection function based on the improved Heaviside function.

[0099] (3) The Lagrangian-Eulerian hybrid description particle flow topology optimization method provided by the present invention has strong applicability and can be applied to optimization under extreme initial design domain conditions. In the classical SIMP method and isogeometric method, when considering the non-uniform layout of the initial design domain density, the optimizer will lose the effectiveness of topological search. For PLSM, the optimization mechanism is to evolve the structure boundary by solving the hamili-jacobi partial differential equation to optimize the structure topology. The distribution of the initial design variables requires an uneven layout, and holes are set in the initial layout, which results in the inapplicability of the PLSM method to the setting where some regions are blank in extreme cases. Therefore, the optimizer can only find a local optimal solution. In the PFTO method proposed in the present invention, even if the particles are placed in extreme positions, the PFTO method can find an optimized topology with a more reasonable material layout. The main reason is that the movement of the particles can quickly and effectively compensate for the defects of the initial distribution through the sensitivity analysis of the objective function and constraint function, and transfer the physical information of the particles to the blank areas of the design region to ensure the generation of the main load transfer path. Therefore, in the PFTO method, the particle flow topology representation model based on the Lagrangian description shows good results in structural topology optimization and can effectively overcome the numerical defects and limitations of multiple traditional methods. At the same time, the solution of the static equation of the design domain adopts the Eulerian description, which makes the optimization process have higher computational efficiency at high numerical accuracy because the same NURBS basis functions are considered in the development of the analysis model. In summary, the PFTO method based on the Lagrangian-Eulerian description demonstrates its powerful effects, superiority, and indispensability, and can effectively find a reasonable material distribution under the desired performance.

[0100] Figure 6 Schematic diagram 1 of the optimization process and results under an extreme initial distribution in the present invention; Figure 7 Schematic diagram 2 of the optimization process and results under an extreme initial distribution in the present invention. As can be seen from Figure 6 and Figure 7 it can be seen that even if the particles are placed in extreme positions, the PFTO method proposed in the present invention can find an optimized topology with a more reasonable material layout.

[0101] Figure 8 Schematic diagram of the optimization process and results of a 1 / 4 circular ring in the present invention. As can be seen from Figure 8 the corresponding optimization example can illustrate the flexibility of the design domain after introducing isogeometric analysis; Figure 9 Schematic diagram of the optimization process and results of a three-dimensional michel structure in the present invention. As can be seen from Figure 9 the corresponding optimization example can illustrate the generality of the topology optimization method for material distribution based on particle flow in the present invention.

[0102] The Lagrangian-Eulerian hybrid description particle flow topology optimization method of the present invention can significantly improve the design flexibility of the topology optimization structure. By introducing isogeometric analysis and the particle flow optimization framework, it overcomes the numerical accuracy and computational efficiency problems faced by traditional methods in the topology optimization process. The implementation of this technology is expected to broaden its application prospects in high-tech fields such as industry, medicine, and aerospace. The commercial value is reflected in the following aspects: Multi-physics and intelligent optimization: The material point method combines the high efficiency of the Eulerian method in solving problems and the flexibility of the Lagrangian method in dealing with large deformations and other issues. In topology optimization, it has extremely high application value, especially in fluid-structure interaction simulation and intelligent design of soft robots.

[0103] Efficient design and manufacturing: Through optimized design methods, material waste can be reduced, manufacturing efficiency can be improved, and production costs can be lowered.

[0104] Cross-field applications: It has broad application potential in fields such as aerospace, medicine, automotive, and architecture, promoting technological innovation and product upgrading in these industries.

[0105] When applying the topology optimization method for material distribution based on particle flow provided by the present invention, it is not necessary to execute according to the order of the steps Figure 1 shown. The execution order of the specific steps can be determined according to needs, and the present invention does not limit this.

[0106] The above is the topology optimization method for material distribution based on particle flow provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding topology optimization device for material distribution based on particle flow, as Figure 10 shown.

[0107] Figure 10 The following is a schematic diagram of a topology optimization device for material distribution based on particle flow provided by the present invention, including: A particle model establishment module 201, configured to initialize a particle design model according to design parameters, so as to initialize the physical information corresponding to the design variables of the particles in the physical space based on the isogeometric analysis grid; A physical information mapping module 202, configured to map the physical information of the particles to the control points of the isogeometric analysis grid through a smooth radial SPH kernel, and perform numerical projection through a control density projection function to obtain the physical density of the control points; A physical density analysis module 203, configured to represent the material distribution in the design domain according to the density distribution function constructed by the control point physical density, and calculate the density value of the Gaussian point; The sensitivity analysis transfer module 204 is used to determine the sensitivities of the objective function and constraint function of the topology optimization to the Gauss points, and transfer the sensitivities of the objective function and constraint function to the Gauss points to the control points according to the sensitivity chain derivative and the density values of the Gauss points, and then transfer them to the design variables of the particles according to the sensitivity chain derivative to determine the sensitivities of the objective function and constraint function to the design variables of the particles; The topology optimization module 205 is used to iteratively update the design variables of the particles by the moving asymptote method according to the sensitivities of the objective function and constraint function to the design variables of the particles, so as to realize the particle flow topology optimization of the material distribution.

[0108] For the specific limitations of the topology optimization device for material distribution based on particle flow, reference can be made to the limitations of the topology optimization method for material distribution based on particle flow in the above text, which will not be elaborated here. Each module in the above topology optimization device for material distribution based on particle flow can be implemented in whole or in part by software, hardware and their combination. The above modules can be embedded in the processor in the computer device in the form of hardware or independent of it, or stored in the memory in the computer device in the form of software, so that the processor can call and execute the operations corresponding to the above modules.

[0109] The present invention also provides a computer-readable storage medium, which stores a computer program, and the computer program can be used to execute the above Figure 1 Provided topology optimization method for material distribution based on particle flow.

[0110] The present invention also provides a computer device. At the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory and a non-volatile memory. Of course, other hardware required for other services may also be included. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it to implement the above Figure 1 Provided topology optimization method for material distribution based on particle flow.

[0111] Those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above methods. Among them, any reference to a memory, storage, database, or other medium used in the various embodiments provided by the present invention can include at least one of non-volatile and volatile memories. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical memory, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0112] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in the present invention.

Claims

1. A topological optimization method for material distribution based on particle flow, characterized in that, Comprising: Initializing a particle design model according to design parameters to initialize physical information corresponding to design variables of particles in the physical space based on an isogeometric analysis grid; Mapping the physical information of the particles to the control points of the isogeometric analysis grid through a smoothed radial SPH kernel, and numerically projecting through a control density projection function to obtain the physical density of the control points; Constructing a density distribution function according to the physical density of the control points to represent the material distribution in the design domain, and calculating the density values of the Gauss points; Determining the sensitivities of the objective function and constraint function of the topology optimization to the Gauss points, and according to the sensitivity chain derivative and the density values of the Gauss points, transferring the sensitivities of the objective function and constraint function to the Gauss points to the control points, and then transferring to the design variables of the particles according to the sensitivity chain derivative to determine the sensitivities of the objective function and constraint function to the design variables of the particles; According to the sensitivities of the objective function and constraint function to the design variables of the particles, iteratively updating the design variables of the particles by the moving asymptote method to realize the particle flow topology optimization of the material distribution.

2. The topological optimization method for material distribution based on particle flow according to claim 1, characterized in that The initializing the particle design model according to the design parameters to initialize the physical information corresponding to the design variables of the particles in the physical space based on the isogeometric analysis grid specifically includes: Initializing the particle design model according to the design parameters to initialize the parameter coordinates of all particles in the parameter space based on the isogeometric analysis grid; Inversely solving the material consumption of the particles according to the volume constraint in the constraint function of the topology optimization; Mapping the parameter coordinates and material consumption of the particles to the physical space through NURBS basis functions.

3. The topology optimization method for material distribution based on particle flow according to claim 1, wherein The mapping the physical information of the particles to the control points of the isogeometric analysis grid through a smoothed radial SPH kernel, and numerically projecting through a control density projection function to obtain the physical density of the control points specifically includes: Using a smoothed radial SPH kernel based on a cubic B-spline kernel function to map the physical information of the particles to the isogeometric analysis grid, solving the control equation on the grid, and obtaining the aggregation variables of the control points through the following formula: ； The cubic B-spline kernel function is: ； The aggregation variables based on the control points are numerically projected by the material point method to determine the physical density of the control points. The numerical projection function of the aggregation variables is as follows: ; In the formula, is the cubic B-spline kernel function, Con is a constant, is the position of the control point, is the position of the particles within the neighborhood of the control point, represents the material consumption of the h th particle, is the SPH kernel size, is the volume of each control point, is the aggregation variable of the control point, and are the coefficients of two functions, is the interval point of the piecewise function.

4. The topology optimization method for material distribution based on particle flow according to claim 1, wherein The constructing a density distribution function according to the physical density of the control points to represent the material distribution in the design domain, and calculating the density values of the Gauss points specifically includes: Constructing a density distribution function according to the physical density of the control points through the following formula to represent the material distribution in the design domain: ； Calculating the density values of the Gauss points according to the following formula: ； Among them, is the density distribution function of the physical density of the control points, , and are the total number of control points in the corresponding parameter directions respectively, is the NURBS basis function in three dimensions, is the th physical density at the control point, is the position of the control point, is the material consumption of the control point, is the density value of the Gauss point, is the density of the control point that affects the Gauss point.

5. The topology optimization method for material distribution based on particle flow according to claim 1, characterized in that The determining the sensitivities of the objective function and constraint function of the topology optimization to the Gauss points, and according to the sensitivity chain derivative and the density values of the Gauss points, transferring the sensitivities of the objective function and constraint function to the Gauss points to the control points specifically includes: Solving the displacement field of the known displacement vector on the Dirichlet boundary on the isogeometric analysis grid, and determining the sensitivities of the objective function and constraint function of the topology optimization to the Gauss points according to the solved displacement field; According to the sensitivity chain derivative and the density values of the Gauss points through the following formula, using the NURBS basis function to transfer the sensitivities of the objective function and constraint function to the Gauss points to the control points to determine the sensitivities of the objective function and constraint function to the control points: ； Among them, is the density value of the Gauss point, is the physical density at the th control point, is the NURBS basis function in three dimensions.

6. The topology optimization method for material distribution based on particle flow according to claim 5, characterized in that Determining the sensitivity of the objective function to the design variables of the particles through the following formula: ； Determining the sensitivity of the constraint function to the design variables of the particles through the following formula: ； Among them, J is the optimization objective, G is the volume constraint, is the displacement field, is the density distribution function of the physical density of the control points, is the design domain, is the strain tensor, is the elastic tensor matrix, is the position of the particle, is the material consumption of the particle, is the material consumption of the particles within the influence domain, is the aggregation variable of the control points, is the cubic B-spline kernel function, r is the relative distance between the particles within the current influence domain and the control points, where h is the SPH kernel size, is the penalty term, is the physical density of the control points.

7. A topological optimization device for material distribution based on particle flow, characterized in that Comprising: A particle model establishment module, configured to initialize a particle design model according to design parameters, so as to initialize physical information corresponding to design variables of particles in a physical space based on an isogeometric analysis grid; A physical information mapping module, configured to map the physical information of particles to control points of an isogeometric analysis grid through a smooth radial SPH kernel, and perform numerical projection through a control density projection function to obtain the physical density of the control points; A physical density analysis module, configured to represent the material distribution in a design domain by a density distribution function constructed according to the physical density of control points, and calculate the density value of Gaussian points; A sensitivity analysis transfer module, configured to determine the sensitivity of an objective function and a constraint function of topology optimization to Gaussian points, and transfer the sensitivity of the objective function and the constraint function to Gaussian points to control points according to the sensitivity chain derivative and the density value of Gaussian points, and then transfer it to design variables of particles according to the sensitivity chain derivative, to determine the sensitivity of the objective function and the constraint function to the design variables of particles; A topology optimization module, configured to iteratively update the design variables of particles by a moving asymptote method according to the sensitivity of the objective function and the constraint function to the design variables of particles, so as to implement particle flow topology optimization of material distribution.

8. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, and when the computer program is executed by a processor, the method according to any one of claims 1 to 6 is implemented.

9. A computer device, characterized in that, It includes a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the method according to any one of claims 1 to 6 is implemented.

Images