A topology optimization method, device, medium and equipment based on CutFEM and SIMP
By combining CutFEM and SIMP, a topology optimization method is developed that solves the computational accuracy and efficiency problems of the SIMP method under complex boundary conditions, achieving high-precision and stable topology optimization, which is suitable for the design of complex structural components.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- DALIAN UNIV OF TECH
- Filing Date
- 2025-09-19
- Publication Date
- 2026-06-09
AI Technical Summary
Existing SIMP methods lack computational accuracy and efficiency under complex geometry and boundary conditions. The presence of intermediate density cells leads to blurred boundaries, and the numerical methods are unstable and dependent on mesh quality, making it difficult to handle irregular boundaries or porous structures.
A topology optimization method combining CutFEM and SIMP is adopted, with node density as the design variable. Node horizontal values are generated through linear mapping, and the nodes are divided into solid, blank, and cut elements. Sub-mesh generation is then performed, inner and outer boundaries are identified, and a total stiffness matrix equation is formed. The design variables are then iteratively optimized.
It improves boundary optimization accuracy and computational stability, adapts to complex geometry and boundary conditions, reduces mesh reconstruction time, reduces numerical oscillation problems, enhances manufacturing feasibility, and is suitable for large-scale structural components.
Smart Images

Figure CN121234745B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of engineering structure topology optimization simulation calculation technology, and in particular to a topology optimization method, apparatus, medium and equipment based on CutFEM and SIMP. Background Technology
[0002] Structural topology optimization is an important technique in modern engineering design, aiming to find the optimal material distribution under given design space, loads, and boundary conditions to maximize structural performance (e.g., minimum flexibility, maximum stiffness, minimum mass). Currently, the most widely used method is density field-based topology optimization, with SIMP (Solid Isotropic Microstructures with Penalization) being a typical example. This method introduces material density variables at elements or nodes and applies nonlinear penalties to material stiffness, thereby driving the optimization results towards a "0-1" discretized design.
[0003] Although the SIMP method has demonstrated good stability and ease of implementation in engineering practice, it also has many limitations, including: 1) the existence of intermediate density cells (i.e., gray areas) blurs the boundaries and seriously affects the manufacturability of the optimization results; 2) increasing the nonlinear penalty factor can suppress intermediate density, but it will exacerbate the numerical instability in the optimization process. In order to obtain clear boundaries, complex processing steps such as filters and projection functions are usually required, which increases the computational cost; 3) the SIMP method is highly dependent on the mesh generation quality and is difficult to maintain accuracy in irregular boundaries or porous structures. Summary of the Invention
[0004] Based on this, it is necessary to propose a topology optimization method, apparatus, medium, and device based on CutFEM and SIMP to address the above problems, aiming to solve the problem of insufficient computational accuracy and efficiency of existing SIMP methods under complex geometry and boundary conditions.
[0005] A topology optimization method based on CutFEM and SIMP includes the following process:
[0006] S1: Define the design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set the material parameters;
[0007] S2: Obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value;
[0008] S3: Based on the node horizontal values, the background mesh unit is divided into solid units, blank units, and cutting units. The solid unit is a unit in which the node horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; the blank unit is a unit in which the node horizontal values of all nodes in the background mesh unit are less than and not equal to 0; the remaining cases are cutting units.
[0009] S4: Identify the material-containing area of the cutting unit and perform sub-mesh division on the material-containing area;
[0010] S5: Identify the boundary lines between the material-containing area and the material-free area of the cutting unit as the inner boundary and the outer boundary;
[0011] S6: Assemble the stiffness matrix in the material-containing region of the solid unit and the cutting unit to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid unit and the cutting unit.
[0012] S7: Generate and solve the topology optimization formula, which includes an objective function and constraints, wherein the objective function is a function of the node displacements;
[0013] S8: Calculate the sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, determine whether the change of the design variables before and after the update is less than a preset value. If yes, complete the optimization and output the design variables. If not, proceed to steps S2~S8 for iterative optimization. If the number of iterations exceeds a predetermined value, complete the optimization and output the design variables.
[0014] A topology optimization device based on CutFEM and SIMP includes:
[0015] An initialization module is used to define a design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set material parameters.
[0016] The topology explicit module is used to obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value.
[0017] A sub-mesh division module is used to divide the background mesh unit into solid units, blank units, and cutting units based on the node horizontal values. Solid units are those where the horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; blank units are those where the horizontal values of all nodes in the background mesh unit are less than and not equal to 0; and the remaining units are cutting units. The module also identifies material-containing areas within the cutting units and performs sub-mesh division on these material-containing areas.
[0018] The inner and outer boundary recognition module is used to identify the boundary line between the material-containing area and the material-free area of the cutting unit as the inner boundary and the outer boundary.
[0019] The finite element calculation module is used to assemble the stiffness matrix in the material region of the solid element and the cutting element to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid element and the cutting element.
[0020] A topology optimization module is used to generate and solve topology optimization formulas, which include objective functions and constraints, wherein the objective function is a function of the node displacements.
[0021] The iteration module is used to calculate sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, and determine whether the change of the design variables before and after the update is less than a preset value. If so, the optimization is completed and the design variables are output. If not, the process proceeds to the explicit topology module, the sub-mesh generation module, the inner and outer boundary identification module, the finite element calculation module, the topology optimization module, and the iteration module for iterative optimization. If the number of iterations exceeds a predetermined value, the optimization is completed and the design variables are output.
[0022] A computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the following steps:
[0023] S1: Define the design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set the material parameters;
[0024] S2: Obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value;
[0025] S3: Based on the node horizontal values, the background mesh unit is divided into solid units, blank units, and cutting units. The solid unit is a unit in which the node horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; the blank unit is a unit in which the node horizontal values of all nodes in the background mesh unit are less than and not equal to 0; the remaining cases are cutting units.
[0026] S4: Identify the material-containing area of the cutting unit and perform sub-mesh division on the material-containing area;
[0027] S5: Identify the boundary lines between the material-containing area and the material-free area of the cutting unit as the inner boundary and outer boundary.
[0028] S6: Assemble the stiffness matrix in the material-containing region of the solid unit and the cutting unit to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid unit and the cutting unit.
[0029] S7: Generate and solve the topology optimization formula, which includes an objective function and constraints, wherein the objective function is a function of the node displacements;
[0030] S8: Calculate the sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, determine whether the change of the design variables before and after the update is less than a preset value. If yes, complete the optimization and output the design variables. If not, proceed to steps S2~S8 for iterative optimization. If the number of iterations exceeds a predetermined value, complete the optimization and output the design variables.
[0031] A computer device includes a memory and a processor, the memory storing a computer program that, when executed by the processor, causes the processor to perform the following steps:
[0032] S1: Define the design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set the material parameters;
[0033] S2: Obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value;
[0034] S3: Based on the node horizontal values, the background mesh unit is divided into solid units, blank units, and cutting units. The solid unit is a unit in which the node horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; the blank unit is a unit in which the node horizontal values of all nodes in the background mesh unit are less than and not equal to 0; the remaining cases are cutting units.
[0035] S4: Identify the material-containing area of the cutting unit and perform sub-mesh division on the material-containing area;
[0036] S5: Identify the boundary lines between the material-containing area and the material-free area of the cutting unit as the inner boundary and the outer boundary;
[0037] S6: Assemble the stiffness matrix in the material-containing region of the solid unit and the cutting unit to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid unit and the cutting unit.
[0038] S7: Generate and solve the topology optimization formula, which includes an objective function and constraints, wherein the objective function is a function of the node displacements;
[0039] S8: Calculate the sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, determine whether the change of the design variables before and after the update is less than a preset value. If yes, complete the optimization and output the design variables. If not, proceed to steps S2~S8 for iterative optimization. If the number of iterations exceeds a predetermined value, complete the optimization and output the design variables.
[0040] Implementing the embodiments of the present invention will have the following beneficial effects:
[0041] This invention first uses node density as a design variable, retaining the advantages of the SIMP method. Then, by converting node density into node level values, it supports complex topological evolutions such as interface splitting and merging without explicit mesh reconstruction, facilitating engineering applications. Subsequently, based on the node level values, the background mesh elements are divided into solid elements, blank elements, and cut elements. Sub-meshing is performed on the cut elements located at the boundaries, which improves the optimization accuracy of the boundaries and is more adaptable to the simulation of complex geometries and boundary conditions. Furthermore, during the iteration process, the background mesh elements remain unchanged, and only the cut elements at the boundaries are sub-meshed, avoiding the reconstruction of meshes for other structural parts. On the one hand, this saves mesh reconstruction time, improves computational efficiency, and is more adaptable to the simulation of large-scale structural components. On the other hand, by sub-meshing the material-containing regions of the cut elements, it can better adapt to changes in free boundaries, improving the stability and accuracy of the calculation. Compared with the unstructured meshes of existing technologies, it can effectively reduce the oscillation problem in numerical solutions during free boundary evolution.
[0042] By using the node horizontal values of all nodes in the background mesh to determine whether a cell is a blank cell, a solid cell, or a cut cell, the gray cells at the boundary can be avoided compared to the SIMP method.
[0043] By identifying internal and external boundaries, the optimization results have higher manufacturing feasibility and are more in line with the needs of actual engineering applications.
[0044] This invention improves the existing SIMP topology optimization method with CutFEM, which not only avoids boundary ambiguity and improves computational efficiency and optimization accuracy, but also enhances computational stability and can handle large-scale structural components and complex boundary problems. Attached Figure Description
[0045] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.
[0046] in:
[0047] Figure 1 This is a flowchart of a topology optimization method based on CutFEM and SIMP in one embodiment;
[0048] Figure 2 This is a schematic diagram of the initial structure within the design domain in one embodiment;
[0049] Figure 3 This is a schematic diagram of the boundary unit processing in one embodiment;
[0050] Figure 4 This is a schematic diagram of a cutting boundary unit in one embodiment;
[0051] Figure 5 This is a schematic diagram of background mesh division and sub-mesh division in one embodiment;
[0052] Figure 6 This is a schematic diagram of boundary identification technology in one embodiment;
[0053] Figure 7 This is the initial structure diagram of Example 1;
[0054] Figure 8 This is a graph showing the optimization results of Example 1;
[0055] Figure 9 This is a structural block diagram of a computer device in one embodiment. Detailed Implementation
[0056] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0057] refer to Figure 1 In one embodiment, a topology optimization method based on CutFEM and SIMP is provided. Specifically, it includes the following steps:
[0058] S1: Define the design domain, discretize the design domain into structured background mesh elements, and set the node density of the background mesh elements. As design variables, material parameters such as elastic modulus and Poisson's ratio are set.
[0059] In this step, discretizing the design domain into structured background mesh elements means discretizing the design domain into finite elements through regular mesh partitioning. In one specific embodiment, the background mesh element is a square element, which is the simplest mesh element and remains unchanged in subsequent calculations. Of course, in other embodiments, the background mesh element can also be other mesh element forms, such as quadrilateral or hexahedral elements.
[0060] refer to Figure 2 Set the design domain to a rectangular outer border, set the material Poisson's ratio to 0.3, initialize the design variables using holes, and specifically set the node density. .
[0061] S2: Node density Obtain node level values by performing a linear mapping of the design domain. Based on node level value The topological structure of the structure is obtained.
[0062] In this step, the node density is... Transformed into node level values through linear mapping This approach reduces the dimensionality of optimization variables and lowers computational complexity. The level set method avoids explicit boundary tracking through implicit interface representations (such as zero-level sets), effectively preventing checkerboard patterns and mesh dependency errors. The zero-value set of the level set function can directly represent structural boundaries, supporting complex topological evolution such as interface splitting and merging without explicit mesh reconstruction. Linear mapping ensures the continuity of the density field gradient, avoiding stress concentration problems caused by abrupt changes in element density and improving structural mechanical performance.
[0063] To further improve computational stability and suppress numerical oscillations, in a preferred embodiment, step S2 above includes the following steps:
[0064] S21: Node density The node density is transformed after filtering. In one specific embodiment, the specific formula is as follows:
[0065]
[0066]
[0067] in, For the filtered nodes i density, It is the node density The weighted average filter function, The nodes of the background grid cells Centered on, with radius R The set of neighboring nodes, There are two nodes. and The linear weight function between them represents the neighborhood nodes. With the central node The closer they are, the greater their weight. It is a node of the background mesh unit. density.
[0068] The above filtering function smooths or filters the original node density value by weighted averaging of the nodes in the neighborhood, resulting in smoother boundaries and a more stable optimization process.
[0069] S22: Filter the node density Node density after projection In one specific embodiment, the specific formula is as follows:
[0070]
[0071] in, It is the projection function of the filtered node density. It is a threshold that determines the center location of the density distribution (usually 0.5). It is the steepness of the projection, which controls the smoothness of the projection curve. The larger the value, the steeper the transition.
[0072] The projection transformation function described above is calculated using a normalized combination of the hyperbolic tangent function (tanh), which is applied to the filtered node density. A nonlinear transformation is performed. The hyperbolic tangent function tanh compresses the input value to the interval [-1, 1], suppressing small local features.
[0073] S23: Projected node density Obtain node level values by performing a linear mapping of the design domain. In one specific embodiment, the specific formula is as follows:
[0074]
[0075] in, For a linear mapping function, A linear transformation from the interval [0,1] to [ , [After normalization] The mapping is applied to the actual engineering scope (such as material thickness, porosity, etc.) to ultimately generate a manufacturable gradient structure.
[0076] In this invention, ,in It is the side length of the unit.
[0077] S3: Based on the node horizontal values, the background mesh elements are divided into solid elements, blank elements, and cut elements. Solid elements are elements in which the node horizontal values of all nodes in the background mesh element are greater than or equal to 0; blank elements are elements in which the node horizontal values of all nodes in the background mesh element are less than or equal to 0; and cut elements are the rest.
[0078] The above method determines the cell type (solid, blank, or cut cell) based on the node's horizontal value to avoid generating gray cells.
[0079] S4: Identify the material-containing regions of the cutting unit, and perform sub-mesh division on the material-containing regions, specifically including the following steps:
[0080] S41: Reference Figure 3 Interpolate the horizontal values of adjacent nodes of the cutting unit to obtain the boundary nodes with a horizontal value of 0 on the edge of the cutting unit, and calculate the side length of the material region, for example, the side length. .
[0081] S42: Obtain the coordinates of the boundary nodes based on the node coordinates of the cutting unit. Specifically, the coordinates of the boundary nodes can be obtained by adding the side length of the material-containing region to the node coordinates.
[0082] S43: Connect the two boundary nodes of a cutting element to obtain the boundary line between the material-containing region and the material-free region of the cutting element.
[0083] S44: Based on the node coordinates of the cutting element, the coordinates of the boundary nodes, and the boundary line, the material region is divided into triangular sub-meshes, such as... Figure 4 As shown. Figure 5 A schematic diagram of the background grid and subgrid division is provided.
[0084] For two-dimensional problems, there are 16 types of boundary line cutting elements, including: Type 1, which has a node horizontal value with one node. If the value is greater than 0, there are four types of cutting lines. Type 2 involves nodes with two adjacent nodes whose horizontal values are equal. A value greater than 0 indicates four types of cutting lines. Type 3 involves nodes with two non-adjacent nodes whose horizontal values are... A value greater than 0 indicates four types of cutting lines. Type four involves a node with three adjacent nodes and a horizontal value of 0. If the value is greater than 0, there are four types of cutting lines. Figure 3 The image shows an example of a cutting process where only the top-left node has a positive value. The first value is positive, while the other three are negative. We are based on four nodes. The value is used to determine the distance of the cutting point using the linearization assumption. It can obtain the node coordinates and boundary information (outer normal, boundary node coordinates) of a solid triangle (with material). (Reference) Figure 4 The element is triangulated based on the coordinates of the element nodes and the coordinates of the boundary nodes.
[0085] To more accurately capture the integration domain within the cutting unit, an adaptive quadtree meshing method can preferably be used to adaptively refine the integration domain, thereby improving integration accuracy. Specifically, the adaptive quadtree method recursively subdivides each cutting unit based on its position, shape, and relationship to the boundary. During the subdivision process, it begins with the boundary region of the cutting unit, based on the nodes... The region is divided using linearization assumptions and then further refined based on whether the sub-regions meet the required precision, until the precision requirement is met or the maximum number of subdivisions is reached. This method allows for finer mesh generation near the boundaries while maintaining computational efficiency. By combining quadtree meshing with linearization assumptions for cell cutting, a finer and more accurate integration domain can be generated near the structural boundaries, thereby improving the accuracy and stability of subsequent numerical calculations.
[0086] S5: Identify the boundary lines between the material-containing and material-free regions of the cutting unit as the inner and outer boundaries.
[0087] In a specific embodiment, the following steps are used to identify the inner and outer boundaries:
[0088] S51: Determine whether all element nodes in the material-free region where the boundary line is located are connected to the outside. If any element node is connected to the outside, then the boundary line is determined to be the outer boundary. If none of the element nodes are connected to the outside, then the boundary line is determined to be the inner boundary.
[0089] To determine whether all element nodes in the material-free region where the boundary line is located are connected to the outside world, the following steps can be used:
[0090] S511: Set up an initial queue of external unit nodes connected to the outside world. The queue of external unit nodes includes at least one external unit node with a horizontal value less than and not equal to 0.
[0091] S512: Starting from any initial external unit node, search for the unit nodes connected to it along a search direction. If the node level value of a unit node is less than and not equal to 0, it is added to the queue of external unit nodes, and the search continues to expand along the search direction to the next unit node. If the node level value of a unit node is greater than 0, the expansion search in that direction is stopped.
[0092] S513: Search other search directions according to the method in step S512;
[0093] S514: Following the method in step S512, search in all directions starting from other initial external unit nodes until all unit nodes have been traversed.
[0094] S515: If a unit node belongs to the queue of external unit nodes, then it is determined that the unit node is connected to the outside world.
[0095] Identifying the boundary lines between the material-containing and material-free regions of the cutting unit as the inner and outer boundaries, in another embodiment, the following method may also be used:
[0096] S51: Obtain the external source direction of the cutting unit. The external source direction is the direction from the external blank unit connected to the cutting unit to the cutting unit. The external blank unit is a blank unit that is connected to the outside world.
[0097] In one specific embodiment, obtaining the external source direction of the cutting unit includes the following steps:
[0098] S511: Set an initial external blank cell queue, which includes at least one initial external blank cell;
[0099] S512: Starting from any initial external blank cell, search for the background grid cells connected to it along a search direction. If the background grid cell is a blank cell, it is added to the queue of external blank cells, and the search continues to expand along the search direction to the next background grid cell. If the background grid cell is a solid / cut cell, the search in that direction is stopped, and the search direction is recorded as the external source direction of the cut cell.
[0100] S513: Search other search directions according to the method in step S512;
[0101] S514: Following the method in step S512, search in all directions starting from other initial external blank cells until all background grid cells have been traversed.
[0102] In the above implementation, the boundary elements are identified while the direction of external origin is marked on the boundary elements, and all boundary elements (i.e. cutting elements) and external blank elements can be identified. The remaining blank elements are blank elements located inside the structure.
[0103] In one specific embodiment, the direction of the external source is represented by a direction vector. For example, "left → right" is recorded as .
[0104] S52: Obtain the outer normal direction of the boundary line of the cutting unit. The outer normal direction is the direction from the material area to the non-material area.
[0105] S53: When the outward normal direction of the boundary line of the cutting unit is opposite to the external source direction of the cutting unit, the boundary line of the cutting unit is identified as the outer boundary; otherwise, it is identified as the inner boundary.
[0106] In the above implementation, the opposite direction of the outward normal to the direction of the external source is not an absolute opposite, but means opposite or basically opposite. For example, in the XY four-quadrant coordinate system, those located in adjacent quadrants are considered to have the same direction, otherwise they are considered to have opposite directions.
[0107] In a specific embodiment of the present invention, when The direction of the outer normal is opposite to the direction of the external source, and the boundary is the outer boundary; otherwise, it is the inner boundary.
[0108] In the two implementation methods described above, the search direction can be up, down, left, or right, or it can be the x or y coordinate direction, and so on.
[0109] This invention, through the two embodiments described above, enables the identification of inner and outer boundaries, ensuring accurate identification of the correct outer and inner boundaries in complex structures. This provides accurate boundary information for subsequent geometric modeling, topology optimization, and manufacturing processes. The identification of inner and outer boundaries is as follows: Figure 6 As shown.
[0110] S6: Assemble the stiffness matrix in the material regions of the solid element and the cut element to form the total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid element and the cut element.
[0111] In one specific embodiment, forming the overall stiffness matrix equations includes the following steps:
[0112] S61: Divide the material area of the cutting unit into triangular meshes to form one or more triangular units.
[0113] S62: Obtain the element stiffness matrix of solid elements and triangular elements based on the equivalent integral weak form.
[0114] S63: Assemble the element stiffness matrix to obtain the overall stiffness matrix equation.
[0115] Specifically, the equivalent weak integral form of the linear elasticity problem is:
[0116]
[0117]
[0118]
[0119] in, For stress tensor, For volume forces, For displacement boundary, As the stress boundary, The direction of the outer normal to the boundary. For displacement.
[0120] Introduction The symmetrization term yields the modified weak form:
[0121]
[0122] The symmetry term ensures that the application of boundary conditions does not affect the variational structure of the weak form, while preserving the nature of the physical problem (such as energy conservation).
[0123] When boundary conditions are applied to the mesh nodes, the equivalent weak integral form is:
[0124]
[0125] in, For virtual strain tensor, Representation domain Integrals on; It is a trial function.
[0126] When the boundary When the boundary does not fit the mesh (e.g., the boundary passes through the element), traditional displacement boundary conditions
[0127] This can lead to numerical instability (such as boundary flux oscillations), so a Nitsche weak form is used with a penalty function. Forced displacement boundary conditions avoid the extra degrees of freedom of the Lagrange multiplier method while satisfying the principle of virtual work:
[0128]
[0129] in, The stability coefficient is typically set to a value of [value missing]. , Let be the element size. To improve the condition number of the stiffness matrix, we introduce the Ghost Penalty term. To prevent pathological problems.
[0130]
[0131] in,
[0132]
[0133] in, For the Ghost Penalty term, the integral edge is a 3D surface, and the common edge between 2D elements is a 2D element. The value is usually between 1 and 10. , This represents the displacement on both sides of the cutting unit.
[0134] For a two-dimensional four-node rectangular element, the stress matrix is:
[0135]
[0136] Among them, subscript Indicates the x and y axis directions. For the elasticity matrix, The strain matrix, , Let be Poisson's ratio, and u be the displacement.
[0137] The displacement matrix is:
[0138]
[0139] in, For shape functions, is the nodal displacement vector.
[0140] The strain-displacement matrix is:
[0141]
[0142] The strain-displacement matrix in the n-direction of the outer normal to the boundary is:
[0143]
[0144] stress
[0145]
[0146] Stress in the n-direction outside the boundary normal:
[0147]
[0148] The matrix expression in the equivalent weak integral form:
[0149]
[0150]
[0151]
[0152]
[0153]
[0154]
[0155]
[0156]
[0157]
[0158]
[0159]
[0160]
[0161] Derivation of the overall stiffness matrix equation for:
[0162]
[0163]
[0164]
[0165] in, U is the element stiffness (Nodal load) The row and column indices (degrees of freedom numbers) are assembled in the overall stiffness matrix K (F terms on the right-hand side). The terms on the right-hand side of the overall stiffness equation represent the strain energy term, boundary flux term, symmetry term, Nitsche stability term, and Ghost Penalty term, respectively. It is the direction of the outer normal to the boundary.
[0166] S7: Generate and solve topology optimization formulas, which include objective functions and constraints. The objective function is a function of nodal displacements.
[0167] Specifically, in one embodiment, the topology optimization formula is as follows:
[0168]
[0169]
[0170]
[0171]
[0172]
[0173] in, It means to The objective function to be minimized is There are m inequality constraints, with the last row corresponding to the constraints on the design variables.
[0174] Objective function: , representing structural flexibility, where, Let be the displacement vectors of the nodes of the solid element and the cut element. For the total stiffness matrix.
[0175] constraint: This indicates that the volume fraction does not exceed the set upper limit of volume. ,in, The number of solid units and cut units, This is the upper limit threshold for the area of the material-containing region in both solid and cut elements. This represents the area of the material-containing region within the solid element and the cut element.
[0176] S8: Calculate sensitivity, obtain the derivatives of the objective function and constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, and determine whether the change of the design variables before and after the update is less than the preset value. If so, the optimization is completed and the design variables are output. If not, proceed to steps S2~S8 for iterative optimization. If the number of iterations exceeds the preset value, the optimization is completed and the design variables are output.
[0177] In a specific embodiment of the present invention, the sensitivity analysis employs a strategy combining the finite difference method and the analytical method to improve computational accuracy and efficiency. Element stiffness and element area are compared with nodes. The derivative is calculated using the finite difference method, and the node... The derivatives with respect to the design variables are obtained analytically. This strategy leverages the advantages of both methods, ensuring the accuracy and efficiency of the sensitivity analysis.
[0178] Specifically, calculating sensitivity includes the following steps:
[0179] S81: Calculate the element stiffness and element area with respect to node horizontal values for solid and cut elements using the finite difference method. The derivative of is given by the following formula:
[0180]
[0181]
[0182] in, It is the element stiffness matrix of the i-th solid element or triangular element. It is the cell area of the i-th solid cell or triangular cell. Is it a solid unit or a triangular unit? The node's horizontal value. yes Tiny perturbations.
[0183] The above step S81 uses the finite difference method for calculation, which can effectively capture the element stiffness and element area relative to the nodal horizontal values. Sensitivity to change.
[0184] S82: Apply the chain rule to calculate the derivative of the node level values of solid elements and cut elements with respect to design variables, as shown in the following formula:
[0185]
[0186] S83: Apply the chain rule to calculate the derivatives of element stiffness and element area with respect to design variables for solid elements and cut elements, as shown in the following formulas:
[0187]
[0188]
[0189] Steps S82 and S83 above employ analytical methods for calculation. By directly deriving the derivatives from the known mathematical relationship between the design variables and the node level values ϕ, the multiple calculations required by the finite difference method are avoided, thus improving computational efficiency. This invention combines the finite difference method with analytical methods to solve for sensitivity, thereby improving computational efficiency while ensuring computational accuracy.
[0190] S84: Based on the results of steps S81 to S83 above, the derivatives of the objective function and constraint function with respect to the design variables can be further derived, as shown in the following formulas:
[0191]
[0192]
[0193] The following is a specific example.
[0194] Calculation example 1
[0195] Standard MBB beam optimization
[0196] Design area: rectangle, length 2, width 1.
[0197] Boundary conditions: Apply displacement constraints to the left mesh nodes.
[0198] Material properties: Elastic modulus E=1, Poisson's ratio ν=0.3
[0199] Nodal force: F=1 at the bottom right corner.
[0200] Objective function: Minimize compliance
[0201] Constraint: Volume fraction V = 0.4 (i.e., approximately 40% of the material region).
[0202] Calculation method: The topology of the MBB beam was optimized using the CutFEM method.
[0203] Optimization results: (attached) Figure 7 Appendix Figure 8 As shown
[0204] This invention also discloses a topology optimization device based on CutFEM and SIMP, comprising:
[0205] The initialization module is used to define the design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as the design variable, and set the material parameters.
[0206] The explicit topology module is used to obtain the node level value by linearly mapping the node density to the design domain, and then obtain the topology of the structure based on the node level value.
[0207] The sub-mesh generation module is used to divide the background mesh cells into solid cells, blank cells, and cut cells based on the node horizontal values. Solid cells are cells in which all nodes of the background mesh cell have a node horizontal value greater than or equal to 0; blank cells are cells in which all nodes of the background mesh cell have a node horizontal value less than or equal to 0; and cut cells are the remaining cases. The module identifies the material regions within the cut cells and performs sub-mesh generation on the material regions.
[0208] The inner and outer boundary recognition module is used to identify the boundary lines between the material-containing and non-material-containing areas of the cutting unit as the inner and outer boundaries.
[0209] The finite element calculation module is used to assemble the stiffness matrix in the material region of solid elements and cutting elements to form the total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of solid elements and cutting elements.
[0210] The topology optimization module is used to generate and solve topology optimization formulas, which include objective functions and constraints. The objective function is a function of nodal displacements.
[0211] The iteration module is used to calculate sensitivity, obtain the derivatives of the objective function and constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, and determine whether the change of the design variables before and after the update is less than the preset value. If so, the optimization is completed and the design variables are output. If not, the module switches to the explicit topology module, submesh generation module, inner and outer boundary identification module, finite element calculation module, topology optimization module, and iteration module for iterative optimization. If the number of iterations exceeds the preset value, the optimization is completed and the design variables are output.
[0212] The specific execution methods for each of the above modules can be found in the aforementioned topology optimization method based on CutFEM and SIMP.
[0213] Figure 9 An internal structural diagram of a computer device in one embodiment is shown. This computer device can specifically be a terminal or a server. Figure 9As shown, the computer device includes a processor, memory, and network interface connected via a system bus. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system and may also store a computer program that, when executed by the processor, enables the processor to implement a topology optimization method. The internal memory may also store a computer program that, when executed by the processor, enables the processor to implement the topology optimization method. Those skilled in the art will understand that... Figure 9 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.
[0214] In one embodiment, a computer device is provided, including a memory and a processor, the memory storing a computer program that, when executed by the processor, causes the processor to perform the aforementioned topology optimization method based on CutFEM and SIMP.
[0215] In one embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, causes the processor to perform the aforementioned topology optimization method based on CutFEM and SIMP.
[0216] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and RAMbus dynamic RAM (RDRAM), etc.
[0217] The technical features of the above embodiments can be combined in any way. For the sake of brevity, 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, they should be considered to be within the scope of this specification.
[0218] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of this patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application. Therefore, the protection scope of this patent application should be determined by the appended claims. Please enter the specific implementation details.
Claims
1. A topology optimization method based on CutFEM and SIMP, characterized in that, Includes the following processes: S1: Define the design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set the material parameters; S2: Obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value; S3: Based on the node horizontal values, the background mesh unit is divided into solid units, blank units, and cutting units. The solid unit is a unit in which the node horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; the blank unit is a unit in which the node horizontal values of all nodes in the background mesh unit are less than and not equal to 0; the remaining cases are cutting units. S4: Identify the material-containing area of the cutting unit and perform sub-mesh division on the material-containing area; S5: Identify the boundary lines between the material-containing area and the material-free area of the cutting unit as the inner boundary and the outer boundary; S6: Assemble the stiffness matrix in the material-containing region of the solid unit and the cutting unit to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid unit and the cutting unit. S7: Generate and solve the topology optimization formula, which includes an objective function and constraints, wherein the objective function is a function of the node displacements; S8: Calculate the sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, determine whether the change of the design variables before and after the update is less than a preset value. If yes, complete the optimization and output the design variables. If not, proceed to steps S2~S8 for iterative optimization. If the number of iterations exceeds a predetermined value, complete the optimization and output the design variables.
2. The topology optimization method based on CutFEM and SIMP according to claim 1, characterized in that, Step S2 includes the following steps: S21: The node density is filtered to convert it into a filtered node density; S22: The filtered node density is transformed into a projected node density through projection; S23: Obtain the node level value by linearly mapping the projected node density to the design domain.
3. The topology optimization method based on CutFEM and SIMP according to claim 2, characterized in that, The filtering transformation function is expressed as follows: in, For the filtered nodes i density, The node density is The weighted average filter function, The nodes of the background mesh unit Centered on, with radius R The set of neighboring nodes, There are two nodes. and The linear weighting function between them It is a node of the background mesh unit. density; The function representing the projection transformation is as follows: in, The projected node density. It is the projection function of the filtered node density. It is a threshold. It is the steepness of the projection; The functional representation of the linear mapping is: in, The horizontal value of the node. It is a linear mapping function. and These are the minimum and maximum values of the node's horizontal value, respectively.
4. The topology optimization method based on CutFEM and SIMP according to claim 1, characterized in that, In step S4, the material-containing regions of the cutting unit are identified, and the material-containing regions are divided into sub-meshes. Includes the following processes: S41: Interpolate the horizontal values of the adjacent nodes of the cutting unit to obtain the boundary nodes with a horizontal value of 0 on the edge of the cutting unit. S42: Obtain the coordinates of the boundary node based on the node coordinates of the cutting unit; S43: Connect the two boundary nodes of a cutting unit to obtain the boundary line between the material-containing area and the material-free area of the cutting unit; S44: Divide the material region into sub-mesh based on the node coordinates of the cutting unit, the coordinates of the boundary node, and the boundary line.
5. The topology optimization method based on CutFEM and SIMP according to claim 1, characterized in that, In step S5, the boundary lines between the material-containing area and the material-free area of the cutting unit are identified as the inner and outer boundaries, including the following steps: S51: Determine whether all unit nodes in the material-free area where the boundary line is located are connected to the outside. If any unit node is connected to the outside, then the boundary line is determined to be the outer boundary. If none of the unit nodes are connected to the outside, then the boundary line is determined to be the inner boundary. Alternatively, in step S5, identifying the boundary line between the material-containing area and the material-free area of the cutting unit as the inner and outer boundaries includes the following steps: S51: Obtain the external source direction of the cutting unit, wherein the external source direction is from the external blank unit connected to the cutting unit to the cutting unit, and the external blank unit is a blank unit connected to the outside. S52: Obtain the outer normal direction of the boundary line of the cutting unit, wherein the outer normal direction is the direction from the area with material to the area without material; S53: When the outward normal direction of the boundary line of the cutting unit is opposite to the external source direction of the cutting unit, the boundary line of the cutting unit is identified as the outer boundary; otherwise, it is identified as the inner boundary.
6. The topology optimization method based on CutFEM and SIMP according to claim 1, characterized in that, In step S7, the objective function is: in, Let be the displacement vector of the nodes of the solid unit and the cutting unit. For the total stiffness matrix, For structural flexibility; The constraints are as follows: in, The number of the solid units and the cutting units. The upper limit threshold for the area of the material-containing region in the solid unit and the cutting unit. The area of the material-containing region in the solid unit and the cutting unit.
7. The topology optimization method based on CutFEM and SIMP according to claim 1, characterized in that, In step S8, the calculation of sensitivity includes the following steps: S81: Calculate the derivatives of the element stiffness and element area of the solid element and the cutting element with respect to the horizontal value of the node using the finite difference method; S82: Apply chain rules to calculate the derivative of the node level value of the solid element and the cutting element with respect to the design variable; S83: Apply the chain rule to calculate the derivatives of the element stiffness and element area of the solid element and the cut element with respect to the design variables; S84: Calculate the derivatives of the objective function and the constraint function with respect to the design variables.
8. A topology optimization device based on CutFEM and SIMP, characterized in that, include: An initialization module is used to define a design domain, discretize the design domain into structured background mesh elements, use the node density of the background mesh elements as a design variable, and set material parameters. The topology explicit module is used to obtain the node level value by linearly mapping the node density to the design domain, and obtain the topology of the structure based on the node level value. A sub-mesh division module is used to divide the background mesh unit into solid units, blank units, and cutting units based on the node horizontal values. Solid units are those where the horizontal values of all nodes in the background mesh unit are greater than and not equal to 0; blank units are those where the horizontal values of all nodes in the background mesh unit are less than and not equal to 0; and the remaining units are cutting units. The module also identifies material-containing areas within the cutting units and performs sub-mesh division on these material-containing areas. The inner and outer boundary recognition module is used to identify the boundary line between the material-containing area and the material-free area of the cutting unit as the inner boundary and the outer boundary. The finite element calculation module is used to assemble the stiffness matrix in the material region of the solid element and the cutting element to form a total stiffness matrix equation, solve the total stiffness matrix equation, and obtain the nodal displacements of the solid element and the cutting element. A topology optimization module is used to generate and solve topology optimization formulas, which include objective functions and constraints, wherein the objective function is a function of the node displacements. The iteration module is used to calculate sensitivity, obtain the derivatives of the objective function and the constraints with respect to the design variables, obtain the updated design variables by solving the topology optimization formula, and determine whether the change of the design variables before and after the update is less than a preset value. If so, the optimization is completed and the design variables are output. If not, the process proceeds to the explicit topology module, the sub-mesh generation module, the inner and outer boundary identification module, the finite element calculation module, the topology optimization module, and the iteration module for iterative optimization. If the number of iterations exceeds a predetermined value, the optimization is completed and the design variables are output.
9. A computer-readable storage medium storing a computer program that, when executed by a processor, causes the processor to perform the steps of the topology optimization method as described in any one of claims 1 to 7.
10. A computer device comprising a memory and a processor, the memory storing a computer program that, when executed by the processor, causes the processor to perform the steps of the topology optimization method as claimed in any one of claims 1 to 7.