CN113312821B - B-spline density method-based topological optimization design method for three-dimensional self-supporting structure - Google Patents
B-spline density method-based topological optimization design method for three-dimensional self-supporting structure Download PDFInfo
- Publication number
- CN113312821B CN113312821B CN202110620489.2A CN202110620489A CN113312821B CN 113312821 B CN113312821 B CN 113312821B CN 202110620489 A CN202110620489 A CN 202110620489A CN 113312821 B CN113312821 B CN 113312821B
- Authority
- CN
- China
- Prior art keywords
- density
- region
- points
- point
- function
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/10—Additive manufacturing, e.g. 3D printing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Complex Calculations (AREA)
Abstract
The invention relates to a three-dimensional self-supporting structure topology optimization design method based on a B-spline density method. By arranging uniform detection points in the structural domain, the relationship between density values at the detection points is controlled, and the elimination of the reverse taper-like region is realized. By constraining the direction of the density field gradient at the design domain boundary, the special V-shaped cross-sectional features present at the design domain boundary are eliminated. The calculation of the inclination angle of the structural boundary in the invention is more convenient and accurate. The suspension angle constraints introduced by the invention are irrelevant to the finite element grids, so that finite element analysis can be carried out by adopting irregularly-shaped units aiming at a complex-shaped design domain in engineering, and the optimal design of the self-supporting structure under any forming direction is realized.
Description
Technical Field
The invention relates to a three-dimensional structure topology optimization design method, in particular to a three-dimensional self-supporting structure topology optimization design method based on a B spline density method.
Background
Additive manufacturing technology is a manufacturing technology that has emerged in recent years to shape structural components by the build-up of materials layer by layer. Compared with the traditional manufacturing technology, the additive manufacturing technology has higher manufacturing freedom degree, and parts with complex shapes can be molded, so that the design space of the structure is enlarged, and the combination of the additive manufacturing technology and the topological optimization design method becomes a hot spot. However, there are still several limitations in additive manufacturing technology, and the dangling portion is one of them. In a conventional additive manufacturing technology, if a suspended portion with a small inclination angle exists in a structure, an auxiliary support structure needs to be added below the suspended portion to prevent collapse in the manufacturing process. The use of support structures results in increased manufacturing costs and time, and may result in damage to the surface of the structure when later removed. Therefore, the self-supporting constraint is introduced into the structural topology optimization design stage, so that the designed structure can be manufactured without adding a support, the practicability of a design result can be effectively improved, and the method has very important engineering application value.
The document "Zhang K, chemical g.three-dimensional high resolution timing dimensioning constraints [ J ] adaptive Manufacturing engineering, 2020,35 101224" discloses a method for the topologically optimised design of three-dimensional self-supporting structures. The method proposed by the literature is based on a density method, firstly adopts an 18-element strategy, and estimates the gradient direction of a density field according to local density distribution around elements in a finite element model, thereby calculating the suspension angle and constructing suspension angle constraint. And secondly, introducing square constraint in the horizontal direction to enable the density of each unit to be larger than the average value of the densities of the units in the four support sets in the horizontal direction, thereby eliminating a suspended area in the structure. By combining the two constraints, the topological optimization design of the three-dimensional self-supporting structure is realized.
Although the method proposed by the literature can obtain a self-supporting structure, the 18-unit strategy adopted by the method obtains the gradient direction of the density field through linear interpolation, and the calculation mode is complex. In addition, the constraints proposed by the method are based on the regular grids and the regular shape design domain, and when the unit shapes and the arrangement modes are irregular or any forming direction is considered, certain limitations exist, so that the method is difficult to apply to engineering practice.
Disclosure of Invention
The invention solves the technical problems that: in order to overcome the defect that the existing topological optimization design method of the three-dimensional self-supporting structure is poor in practicability, the invention provides a topological optimization design method of the three-dimensional self-supporting structure based on a B spline density method. According to the method, three relevant constraints are constructed on the basis of a B-spline parameterized density field so as to realize the optimal design of a self-supporting structure. Firstly, obtaining the gradient of a structure density field by derivation of a B spline function, identifying a boundary region of the structure according to a gradient mode and the curvature radius of a space curve, and controlling the surface inclination angle of the structure in the boundary region to be larger than a specified critical angle of equipment. Secondly, according to different forming directions, a cuboid region containing a structural domain is constructed, a certain number of uniformly distributed detection points are arranged in the region, and the relationship between density values of the detection points is controlled to eliminate the reverse-tapered-like region. And finally, eliminating the special V-shaped section characteristics at the boundary of the design domain by restricting the gradient direction of the density field at the boundary of the design domain.
The technical scheme of the invention is as follows: a three-dimensional self-supporting structure topological optimization design method based on a B spline density method comprises the following steps:
step 1: structuring the structural density field: building the structural design Domain omega of the white support Structure d Defining a level set function phi representing the region d And constructing a density field rho representing the structural design domain d (ii) a Designing Domain Ω according to Structure d The shape of (2), the structural design domain omega d A rectangular parallelepiped region Ω included therein; constructing a three-dimensional B-spline function field inside a cuboid region omega as a B-spline density field rho of the region; and (4) performing Heaviside projection on the B spline density field rho to obtain a projected B spline density field
Wherein β is the Heaviside projectionη is the projection limit of the Heaviside projection; projecting the B-spline density field
And a density field ρ representing a structural design domain d Performing Boolean intersection operation to obtain the final density field of the structure
Step 2: identifying the boundary region of the structure:
and step 3: constructing inclination angle constraint: for the boundary region omega obtained in step 2 b Arranging a plurality of integration points therein and calculating the gradient of the B-spline density field rho at each integration point
And further calculating cosine values cos alpha of the boundary inclination angles alpha at the integration points:
wherein b is a unit vector with a direction pointing in the direction of additive manufacturing molding,
a gradient mode representing the B-spline density field; according to the inclination angle alpha and the critical angle alpha specified by the equipment 0 Difference of cosine values of (a) in a boundary region omega of the structure b And (3) internally integrating to obtain an inclination angle constraint value q:
where L (δ) is a piecewise function with respect to a variable δ, L (δ) =0 when δ < 0; when the value of delta is greater than 0,L(δ)＝δ 2 (ii) a And establishing an inclination angle constraint function by using the inclination angle constraint value q:
wherein epsilon q Is an extremely small number introduced for relaxing the tilt angle constraint, and the tilt angle constraint is applied to ensure that the tilt angle of any point at the structure boundary is larger than the specified critical angle of the equipment;
and 4, step 4: the elimination of the reverse taper-like region in the optimization result is realized, and the method comprises the following substeps:
step 4.1: constructing a new rectangular parallelepiped region Ω including the rectangular parallelepiped region Ω according to the molding direction of additive manufacturing t ：Ω t One side of the rectangular area is parallel to the additive manufacturing forming direction, and the other two sides are not limited and can contain the rectangular area omega;
and 4.2: in a rectangular region omega t A plurality of detection points are uniformly arranged along the directions of the three edges, and the detection points are positioned in a structural design domain omega d The detection point at one side of the bottom edge towards the additive manufacturing and forming direction is an effective detection point; final density field according to structure
Calculating the density value of the effective detection point, wherein the density value is positioned in the structural design field omega d The density value of the outer effective detection point is 0 to obtain the density value vector of the effective detection point
Step 4.3: in the effective detection points, n is remained after excluding the lowest detection point in the additive manufacturing molding direction t The detection points are detection points needing to participate in calculation, and for the ith detection point needing to participate in calculation, the density value of the point is constrained
Should be less than n in the region of the pyramid below it c The maximum of the density values at the points is formulated as:
wherein n is c The point refers to n in the layer below the ith detection point needing to participate in calculation c The points comprise a plurality of points in a square area taking a point right below the ith detection point needing to participate in calculation as a central point, and the n points are formed together c The point of the light source is a point,
respectively representing the density values of nc points in a pyramid region below the ith detection point, the first item i in the subscript representing the corresponding serial number of the detection point at the top of the pyramid region, and the second items 1-n c Respectively representing n inside the pyramid region c And (3) in order to facilitate optimization by using a gradient algorithm, approximating a max function in the above formula by using a KS function, and calculating an integral constraint value:
where w is a parameter in the KS function,. Epsilon 1 Is a very small number introduced to balance the approximation error of the KS function, and h (δ) is a function of the variable δ in the form:
wherein the magnitude of the parameter mu determines the severity of the change of the function at delta =0,
step 4.4: with this constraint value g, an overall constraint function can be established:
wherein epsilon 2 Is a very small number introduced to relax the tilt constraint; by applying the constraint, the elimination of the reverse-tapered-like region in the optimization result is realized;
and 5: for the domain located in the structural design Ω d Is of d A point with an upper coordinate of (x, y, z), which is perpendicular to the structure boundary and points to the interior of the structure, can be represented as
Let the plane perpendicular to the additive manufacturing forming direction and passing through the point (x, y, z) be denoted as m 3 Then the level set function phi d Gradient of (2)
In plane m 3 The projection vector d on can be expressed as:
and ensure that:
according to the projection vector d of a plurality of integral points on the boundary of the structural design domain and the gradient of the B-spline density field
At structural design domain boundary Γ d And integrating to obtain an integral constraint value:
wherein the L (δ) function is the same as in step 3; using the constraint value, an overall constraint function can be established:
by applying the constraint, the elimination of the special V-shaped section characteristics at the boundary of the structural design domain in the optimization result is realized;
and 6: in the structural design domain omega d Dispersing the materials into finite element units, and defining load and boundary conditions simultaneously;
and 7: the topology optimization problem is defined as:
where P denotes the value of the control point of the B-spline function field constructed in the rectangular region omega, n x 、n y 、n z Respectively representing the number of control points in three directions of a cuboid region omega; p is i The middle subscript i represents the serial number of the control point; c represents the degree of structural compliance;
respectively representing the final density fields of three different structures, and respectively representing the total rigidity matrix, the total load vector and the displacement vector of the structure by K, F and U; v 0 Denotes the total volume of the structure, V d And
respectively representing projection limits of eta d Time structure final density field
Current volume fraction and volume fraction upper limit;
and 8: and (6) solving the topological optimization problem in the step (6) by adopting an optimization algorithm to obtain an optimal result.
The further technical scheme of the invention is as follows: the three density fields in step 7 are obtained by selecting three different projection limits η when performing the Heaviside projection in step 1 e 、η i 、η d And obtaining the product, wherein the values of the product are between 0 and 1.
Effects of the invention
The invention has the technical effects that: compared with the design method in the background art, the density field gradient is directly obtained through derivation calculation, so that the calculation of the inclination angle of the structure boundary is more convenient and accurate. According to the method, a B-spline function is used for constructing a structure density field, the gradient of the density field is obtained through derivation, a structure boundary area is further identified, and a boundary inclination angle is obtained in the boundary area according to the gradient direction of the density field, so that the size of a suspension angle is controlled. Meanwhile, the elimination of the reverse taper-like region is realized by arranging uniform detection points in the structural domain and controlling the relation between density values at the detection points. Finally, by constraining the direction of the density field gradient at the design domain boundary, the special V-shaped cross-sectional features appearing at the design domain boundary are eliminated. Compared with the design method of the background art, the gradient of the density field is directly obtained through derivation calculation, so that the calculation of the inclination angle of the structure boundary is more convenient and accurate. In addition, the suspension angle constraints introduced by the invention are irrelevant to the finite element grids, so that the finite element analysis can be carried out by adopting units with irregular shapes aiming at the design domain with complicated shapes in the engineering, the optimal design of the self-supporting structure in any forming direction is realized, and the invention has stronger engineering practicability. Without considering self-support, the structural compliance obtained with the same initial structural optimization was 43.4188, thus achieving structural self-support at the expense of 12.28% compliance.
Drawings
FIG. 1 is a diagram illustrating the geometry and boundary conditions of a model according to an embodiment of the present invention.
FIG. 2 is a diagram illustrating distribution of detection points within a domain according to an embodiment of the present invention.
FIG. 3 is a schematic diagram of the detection points and the pyramid area below the detection points in an embodiment of the present invention.
FIG. 4 is a diagram of the design results of an embodiment of the present invention.
Detailed Description
In the description of the present invention, it is to be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", and the like, indicate orientations or positional relationships based on those shown in the drawings, merely for convenience of description and simplification of the description, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed in a particular orientation, and be operated, and thus, are not to be construed as limiting the present invention.
Referring to fig. 1-4, the technical solution adopted by the present invention to solve the technical problem is as follows: a three-dimensional self-supporting structure topological optimization design method based on a B spline density method is characterized by comprising the following steps:
step 1: a structural density field is constructed. Building a structural design Domain omega of a self-supporting Structure d Writing a level set function phi representing the region d And constructing a density field rho representing the structural design domain d . Design of Domain Ω according to Structure d Is configured to design the structure into a domain omega d A rectangular parallelepiped region Ω included therein. And constructing a three-dimensional B-spline function field inside the cuboid region omega as a B-spline density field rho of the region. Subjecting the B-spline density field rho to Heaviside projection to obtain a projected B-spline density field
Where β is the parameter of the Heaviside projection and η is the projection limit of the Heaviside projection. The B-spline density field after projection
And a density field ρ representing a structural design domain d Performing Boolean intersection operation to obtain the final density field of the structure
Step 2: a boundary region of the structure is identified. For a point with coordinates (x, y, z) in the structural design domain, its final density field
Gradient of (2)
And gradient mode
Can be obtained by derivation of B-spline functions:
wherein x, y, z are the coordinates in a cartesian coordinate system, respectively. A lower limit of a gradient mode is preset
When the gradient norm at point (x, y, z) is greater than the lower limit, it is temporarily identified as a boundary. Structural boundaries omega identified from gradient modes b Can be expressed as:
in order to eliminate the small corner region pair above the hole in the structure or at the intersection of two rodsThe influence of the tilt angle constraint to be applied requires the above-mentioned structural boundary Ω b The region with the smaller inner radius of curvature is excluded. Also for a point with coordinates (x, y, z) in the design domain of the structure, the iso-surface m1 of the final density field of the structure at that point and the gradient at that point are calculated
Plane m formed by additive manufacturing and forming direction 2 M can be obtained by a curvature radius calculation formula of the space curve 1 And m 2 The radius of curvature r (x, y, z) at the point (x, y, z). A lower limit r of the radius of curvature is specified min Then the boundary region omega after the small fillet is removed b Can be re-expressed as:
and step 3: and constructing an inclination angle constraint. For the boundary region omega obtained in step 2 b Arranging a plurality of integration points therein and calculating gradients of the B-spline density field rho at each integration point
And further calculating cosine values cos alpha of the boundary inclination angles alpha at the integration points:
wherein b is a unit vector with a direction pointing in the direction of additive manufacturing molding,
the gradient mode of the B-spline density field is represented. According to the inclination angle alpha and the critical angle alpha specified by the equipment 0 Difference of cosine values of (a) in a boundary region omega of the structure b And (3) internally integrating to obtain an inclination angle constraint value q:
where L (δ) is a piecewise function with respect to a variable δ, L (δ) =0 when δ < 0; when δ > 0, L (δ) = δ 2 . And establishing an inclination angle constraint function by using the inclination angle constraint value q:
wherein epsilon q Is a very small number introduced to relax the tilt constraint. By applying the tilt angle constraint, the tilt angle of any point at the structure boundary is made larger than the critical angle specified by the device.
And 4, step 4: constructing a new rectangular parallelepiped region Ω including the rectangular parallelepiped region Ω according to the molding direction of additive manufacturing t ，Ω t One side of the rectangular area is parallel to the additive manufacturing and forming direction, and the other two sides are not limited and can contain the rectangular area omega. In a rectangular region omega t A plurality of detection points are uniformly arranged along the directions of three edges and positioned in a structural design domain omega d The detection point at one side of the bottom edge towards the additive manufacturing and forming direction is an effective detection point. Final density field according to structure
Calculating the density value of the effective detection point, wherein the density value is positioned in the structural design domain omega d The density value of the outer effective detection point is 0 to obtain the density value vector of the effective detection point
In the effective detection points, n is remained after excluding the lowest detection point in the additive manufacturing molding direction t The detection points are detection points which need to participate in calculation.
For the ith detection point needing to participate in calculation, the density value at the point is constrained
Should be smaller than the pyramid area below itN within a domain c The maximum of the density values at the points is formulated as:
the lower region here means a lower region defined in the additive manufacturing forming direction. n is a radical of an alkyl radical c The point refers to n in the layer below the ith detection point needing to participate in calculation c The points comprise a plurality of points in a square area taking a point right below the ith detection point needing to participate in calculation as a central point, and the n points are formed together c And (4) points.
Respectively represents n in a pyramid region below the ith detection point c Density value of points, the first item i in subscript represents the corresponding serial number of the detection point at the top of pyramid region, and the second items 1-n c Respectively representing n inside the pyramid region c And detecting points. To facilitate optimization using the gradient algorithm, the KS function is used to approximate the max function in the above equation and the overall constraint value is calculated:
where w is a parameter in the KS function,. Epsilon 1 Is a very small number introduced to balance the KS function approximation error, and h (δ) is a function of the variable δ, of the form:
wherein the magnitude of the parameter μ determines the severity of the change of the function at δ =0. With this constraint value g, an overall constraint function can be established:
wherein epsilon 2 Is a very small number introduced to relax the tilt constraint. By applying this constraint, elimination of the reverse-tapered-like region in the optimization result is achieved.
And 5: for a domain located in the structural design domain Ω d Is of d A point with an upper coordinate (x, y, z), which is perpendicular to the structure boundary and points in the direction of the structure interior, can be represented as
I.e. the level set function phi representing the structural design domain d Of the gradient of (c). Let the plane perpendicular to the additive manufacturing forming direction and passing through the point (x, y, z) be denoted as m 3 Then the level set function phi d Gradient of (2)
In plane m 3 The projection vector d on can be expressed as:
in order to avoid the appearance of special V-shaped section characteristics at the boundary of the structural design domain, the gradient of the B-spline density field rho at the boundary of the design domain needs to be ensured
The dot product with the projection vector d is equal to or greater than 0, that is:
according to the projection vector d of a plurality of integral points on the boundary of the structural design domain and the gradient of the B-spline density field
At structural design domain boundary Γ d And integrating to obtain an integral constraint value:
where the L (δ) function is the same as in step 3. Using the constraint values, an overall constraint function can be established:
by applying the constraint, elimination of special V-shaped section characteristics at the boundary of the structural design domain in the optimization result is realized.
Step 6: in the structural design domain omega d Discretized into finite element elements, defining load and boundary conditions.
And 7: the topology optimization problem is defined as:
where P denotes the value of the control point of the B-spline function field constructed in the rectangular region omega, n x 、n y 、n z The number of control points in three directions of the rectangular region Ω is respectively indicated. P i The middle subscript i indicates the number of the control points, and all control points have values between 0 and 1. C represents the degree of structural compliance.
Three different structural final density fields are respectively represented, which are obtained by selecting three different projection limits eta when performing the Heaviside projection in step 1 e 、η i 、η d So as to obtain the compound with the characteristics of,the values are all between 0 and 1. K. F and U represent the overall stiffness matrix, the overall load vector and the displacement vector of the structure, respectively. V 0 Denotes the total volume of the structure, V d And
respectively representing projection limits of eta d Time structure final density field
And a current volume fraction of (d) and a volume fraction upper limit.
And 8: and (6) solving the topological optimization problem in the step (6) by adopting an optimization algorithm to obtain an optimal result.
Reference is made to fig. 1-4. The invention is illustrated by taking a three-dimensional simple beam structure under fixed load as an example. Due to the symmetry of the simply supported beam, only the right half part of the simply supported beam is taken for optimal design. Its Young's modulus is 1 and Poisson's ratio is 0.3. The molding direction defining the structure is represented by a vector b = (0, 1), and the critical molding angle of the printing apparatus is α 0 =45 °. The material distribution in the simply supported beam is selected and designed to enable the flexibility of the simply supported beam to be minimum, and the volume fraction of the total material consumption is maximum 20%. The method comprises the following specific steps:
1. a structural density field is constructed. First, write out the structural design domain Ω d Level set function of (phi) d Density field rho of structural design domain d . Constructing the rectangular region Ω will design the region Ω d And (4) carrying out interpolation by adopting the B spline basis function of the uniform node for 3 times to obtain a B spline density field rho in the cuboid region omega. Performing Heaviside projection on the B-spline density field, wherein the projection parameter beta is 1, the projection limit eta is 0.5, and obtaining the projected B-spline density field
Projecting the B-spline density field
And the design domain density field ρ d Performing Boolean intersection operation to obtain the final density field of the structure
2. Structural boundaries are identified. The lower limit of the gradient mode predefined in this example is
Lower limit of curvature radius r min =2, thus the boundary region Ω b Is shown as
3. And constructing an inclination angle constraint. In the boundary region omega b Selecting integral points, calculating the gradient at the integral points by derivation and further calculating the cosine value of the inclination angle, wherein the calculation formula is
Based on the difference between the tilt angle cosine and the cosine of the equipment-specified critical angle, the method is performed in a region omega b And (3) internally integrating to construct a dip angle constraint function:
the allowable region Ω in this embodiment b The angle of inclination of the inner part is subject to a tolerance of 0.5 deg., so that q Is calculated as
4. Constructing a new rectangular region omega including the rectangular region omega according to the molding direction of additive manufacturing t . The new cuboid region omega in this example t The same as the rectangular parallelepiped region Ω. In a rectangular region omega t A plurality of detection points are uniformly arranged along the directions of the three edges,the total number of detection points in this example is 60 × 24 × 60, and all the detection points are valid detection points. Final density field according to structure
Calculating to obtain the density values at the effective detection points
Wherein is located in the structural design domain Ω d The density value of the outer effective detection point is 0. The number of detection points which are left to participate in calculation after the lowest detection point in the effective detection points is removed is 60 multiplied by 24 multiplied by 59. The overall constraints are:
in the formula n t ＝84960，w＝20，n c ＝25，ε 1 ＝0.15，ε 2 ＝5×10 -6 . The value of the parameter μ in the h (δ) function is μ =40.
5. At structural design domain boundary Γ d And selecting an integral point and calculating a projection vector d of the integral point. According to the projection vector d of a plurality of integral points on the boundary of the structural design domain and the gradient of the B-spline density field
At structural design domain boundary Γ d And integrating, constructing an integral constraint:
6. design domain omega of structure by adopting patch grids d Discretizing and defining relevant boundary conditions.
7. Defining a topology optimization problem:
where P represents the value of all control points, with an initial value of 0.5 and upper and lower limits of 1 and 0, respectively. n is x ＝30，n y ＝12，n z =30. Three different structural final density fields
By selecting the projection limit η at the Heaviside projection, respectively e ＝0.7、η i ＝0.5、η d And = 0.3. The projection parameter β of the Heaviside projection takes 1 at 100 steps before optimization and then multiplies by a factor of 1.7 every 100 steps until it reaches its maximum value 25.
8. Finite element analysis and optimization solution. And carrying out finite element analysis on the established model, and respectively carrying out sensitivity analysis on the target function and the constraint function. And (4) selecting a global convergence moving asymptote optimization algorithm (GCMMA) to carry out optimization design to obtain an optimization result.
The method obtains an optimized design result after iteration of the step 450 of the embodiment. The volume fraction of the initial structure is 20%, the flexibility is 842.0703, and the volume fraction of the structure in the optimized result is 20%, the flexibility is 48.7511, and the flexibility is reduced by 94.21%. Without considering self-support, the structural compliance achieved with the same initial structural optimization was 43.4188, thus achieving a structural self-support sacrifice of 12.28% compliance. Compared with the design method in the background art, the density field gradient is directly obtained through derivation calculation, so that the calculation of the inclination angle of the structure boundary is more convenient and accurate. In addition, the suspension angle constraints introduced by the invention are irrelevant to the finite element grids, so that the finite element analysis can be carried out by adopting units with irregular shapes aiming at the design domain with complicated shapes in the engineering, the optimal design of the self-supporting structure in any forming direction is realized, and the invention has stronger engineering practicability.
Claims (2)
1. A three-dimensional self-supporting structure topological optimization design method based on a B spline density method is characterized by comprising the following steps:
step 1: structuring the structural density field: building a structural design domain omega of a self-supporting structure d Defining a level set function phi representing the region d And constructing a density field rho representing the structural design domain d (ii) a Designing Domain Ω according to Structure d The structural design domain omega is constructed d A rectangular parallelepiped region Ω included therein; constructing a three-dimensional B-spline function field inside a cuboid region omega as a B-spline density field rho of the region; subjecting the B-spline density field rho to Heaviside projection to obtain a projected B-spline density field
Wherein β is a parameter of the Heaviside projection and η is a projection limit of the Heaviside projection; projecting the B-spline density field
And a density field ρ representing a structural design domain d Performing Boolean intersection operation to obtain the final density field of the structure
Step 2: identifying the boundary region of the structure:
and 3, step 3: constructing inclination angle constraint: for the edge obtained in step 2Boundary region omega b Arranging a plurality of integration points therein and calculating the gradient of the B-spline density field rho at each integration point
And further calculating cosine values cos alpha of the boundary inclination angles alpha at the integration points:
wherein b is a unit vector with a direction pointing in the direction of additive manufacturing molding,
a gradient mode representing the B-spline density field; according to the inclination angle alpha and the specified critical angle alpha of the equipment 0 The difference between cosine values of (a) in a boundary region omega of the structure b And (3) internally integrating to obtain an inclination angle constraint value q:
where L (δ) is a piecewise function with respect to a variable δ, when δ < 0, L (δ) =0; when δ > 0, L (δ) = δ 2 (ii) a And establishing an inclination angle constraint function by using the inclination angle constraint value q:
wherein epsilon q The minimum number is introduced for relaxing the inclination angle constraint, and the inclination angle constraint is applied to ensure that the inclination angle of any point at the structure boundary is larger than the specified critical angle of the equipment;
and 4, step 4: the elimination of the reverse taper-like region in the optimization result is realized, and the method comprises the following substeps:
step 4.1: constructing a new rectangular region including a rectangular region omega according to a molding direction of additive manufacturingBody region omega t ：Ω t One side of the rectangular area is parallel to the additive manufacturing forming direction, and the other two sides are not limited and can contain the rectangular area omega;
step 4.2: in a rectangular region omega t A plurality of detection points are uniformly arranged along the directions of the three edges, and the detection points are positioned in a structural design domain omega d The detection points on one side of the bottom edge towards the additive manufacturing forming direction are effective detection points; final density field according to structure
Calculating the density value of the effective detection point, wherein the density value is positioned in the structural design domain omega d The density value of the outer effective detection point is 0 to obtain the density value vector of the effective detection point
Step 4.3: in the effective detection points, n is remained after excluding the lowest detection point in the additive manufacturing molding direction t The detection points are detection points needing to participate in calculation, and for the ith detection point needing to participate in calculation, the density value of the point is constrained
Should be less than n in the region of the pyramid below it c The maximum of the density values at the points is formulated as:
wherein n is c The point refers to n in the next layer of the ith detection point needing to participate in calculation c The points comprise a plurality of points in a square area taking the point right below the ith detection point needing to participate in calculation as a central point, and the n points are formed together c The point of the light beam is the point,
respectively representN in pyramid area below ith detection point c Density value of points, the first item i in subscript represents the corresponding serial number of the detection point at the top of pyramid region, and the second items 1-n c Respectively representing n inside the pyramid region c And (3) in order to facilitate optimization by using a gradient algorithm, approximating a max function in the above formula by using a KS function, and calculating an integral constraint value:
where w is a parameter in the KS function,. Epsilon 1 Is a very small number introduced to balance the KS function approximation error, and h (δ) is a function of the variable δ, of the form:
wherein the magnitude of the parameter mu determines the severity of the change of the function at delta =0,
step 4.4: with this constraint value g, an overall constraint function can be established:
wherein epsilon 2 Is a very small number introduced to relax the tilt constraint; by applying the constraint, the elimination of the reverse-tapered-like region in the optimization result is realized;
and 5: for the domain located in the structural design Ω d Is of d A point with an upper coordinate (x, y, z), which is perpendicular to the structure boundary and points in the direction of the structure interior, can be represented as
Let the plane perpendicular to the additive manufacturing forming direction and passing through the point (x, y, z) be denoted as m 3 The gradient of the level set function phid
In plane m 3 The projection vector d on can be expressed as:
and ensure that:
according to the projection vector d of a plurality of integral points on the boundary of the structural design domain and the gradient of the B-spline density field
At structural design domain boundary Γ d And integrating to obtain an integral constraint value:
wherein the L (δ) function is the same as in step 3; using the constraint values, an overall constraint function can be established:
by applying the constraint, the elimination of the special V-shaped section characteristics at the boundary of the structural design domain in the optimization result is realized;
and 6: dispersing the structure design domain omega d into finite element units, and defining load and boundary conditions at the same time;
and 7: the topology optimization problem is defined as:
where P denotes the value of the control point of the B-spline function field constructed in the rectangular region omega, n x 、n y 、n z Respectively representing the number of control points in three directions of a cuboid region omega; p i The middle subscript i represents the serial number of the control point; c represents the degree of structural compliance;
respectively representing the final density fields of three different structures, and respectively representing the total rigidity matrix, the total load vector and the displacement vector of the structure by K, F and U; v 0 Denotes the total volume of the structure, V d And
respectively representing projection limits of eta d Time structure final density field
Current volume fraction and volume fraction upper limit;
and 8: and (6) solving the topological optimization problem in the step (6) by adopting an optimization algorithm to obtain an optimal result.
2. The method for topologically optimally designing a three-dimensional self-supporting structure based on the B-spline density method as claimed in claim 1, wherein the three density fields in the step 7 are obtained by selecting three different projection limits η when performing the Heaviside projection in the step 1 e 、η i 、η d And obtaining the product, wherein the values of the product are between 0 and 1.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110620489.2A CN113312821B (en) | 2021-06-03 | 2021-06-03 | B-spline density method-based topological optimization design method for three-dimensional self-supporting structure |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110620489.2A CN113312821B (en) | 2021-06-03 | 2021-06-03 | B-spline density method-based topological optimization design method for three-dimensional self-supporting structure |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113312821A CN113312821A (en) | 2021-08-27 |
| CN113312821B true CN113312821B (en) | 2023-01-24 |
Family
ID=77377242
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110620489.2A Active CN113312821B (en) | 2021-06-03 | 2021-06-03 | B-spline density method-based topological optimization design method for three-dimensional self-supporting structure |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113312821B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116522562B (en) * | 2023-05-22 | 2024-01-30 | 上海祉元社企业管理合伙企业(有限合伙) | Supporting structure construction method, system, electronic equipment and storage medium |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013207785A (en) * | 2012-03-29 | 2013-10-07 | Toshiba Corp | Video processing apparatus and video processing method |
| WO2014074954A2 (en) * | 2012-11-08 | 2014-05-15 | Suman Das | Systems and methods for fabricating three-dimensional objects |
| CN110096756A (en) * | 2019-04-08 | 2019-08-06 | 河海大学 | A kind of free form surface structural form creation method considering load uncertainties |
| CN111319268A (en) * | 2020-02-20 | 2020-06-23 | 西北工业大学 | Self-supporting structure optimization design method considering additive manufacturing printing direction |
| CN111814383A (en) * | 2020-07-25 | 2020-10-23 | 西北工业大学 | B-spline density method-based self-supporting structure topology optimization design method |
| TW202116072A (en) * | 2019-10-01 | 2021-04-16 | 美商杜拜研究特許公司 | Tensor-product b-spline predictor |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9561622B2 (en) * | 2008-05-05 | 2017-02-07 | Georgia Tech Research Corporation | Systems and methods for fabricating three-dimensional objects |
-
2021
- 2021-06-03 CN CN202110620489.2A patent/CN113312821B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013207785A (en) * | 2012-03-29 | 2013-10-07 | Toshiba Corp | Video processing apparatus and video processing method |
| WO2014074954A2 (en) * | 2012-11-08 | 2014-05-15 | Suman Das | Systems and methods for fabricating three-dimensional objects |
| CN110096756A (en) * | 2019-04-08 | 2019-08-06 | 河海大学 | A kind of free form surface structural form creation method considering load uncertainties |
| TW202116072A (en) * | 2019-10-01 | 2021-04-16 | 美商杜拜研究特許公司 | Tensor-product b-spline predictor |
| CN111319268A (en) * | 2020-02-20 | 2020-06-23 | 西北工业大学 | Self-supporting structure optimization design method considering additive manufacturing printing direction |
| CN111814383A (en) * | 2020-07-25 | 2020-10-23 | 西北工业大学 | B-spline density method-based self-supporting structure topology optimization design method |
Non-Patent Citations (5)
| Title |
|---|
| Shapes modeling of 3D objects based on a hybrid representation using extended B-spline surface model;M. Maeda et al;《IEEE》;20021210;全文 * |
| 叶轮快速成型及熔模铸造过程仿真研究;张皓;《中国优秀硕士学位论文全文数据库（电子期刊）（工程科技I辑）》;20200315;全文 * |
| 基于固定网格和拓扑导数的结构拓扑优化自适应泡泡法;蔡守宇 等;《力学学报》;20190430;第51卷(第4期);1235-1244 * |
| 特征驱动的结构拓扑优化方法;张卫红 等;《中国科学：技术科学》;20191031;第49卷(第10期);1177-1185 * |
| 面向增材制造的拓扑优化技术发展现状与未来;朱继宏 等;《航空制造技术》;20201031(第10期);24-38 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113312821A (en) | 2021-08-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109002598B (en) | Self-supporting microstructure topology optimization method considering overhanging angle and minimum size constraint | |
| Xu et al. | Curved layer based process planning for multi-axis volume printing of freeform parts | |
| CN111319268B (en) | Self-supporting structure optimization design method considering additive manufacturing printing direction | |
| Thrimurthulu et al. | Optimum part deposition orientation in fused deposition modeling | |
| CN110517193B (en) | Seabed sonar point cloud data processing method | |
| CN111444640B (en) | Structural topology optimization method considering inclination constraint of additive manufacturing | |
| CN110795873B (en) | Cross-scale topology optimization method considering size control | |
| Sabiston et al. | 3D topology optimization for cost and time minimization in additive manufacturing | |
| CN107391824B (en) | Topological optimization design method of self-supporting structure in additive manufacturing | |
| CN110543654B (en) | Method for determining distributed process parameters of laser shot blasting forming complex curved surface | |
| CN107341316B (en) | Structural shape-topology combined optimization method under design related pressure load effect | |
| CN111475976B (en) | Robust topology optimization method for particle reinforced material member considering mixing uncertainty | |
| CN113312821B (en) | B-spline density method-based topological optimization design method for three-dimensional self-supporting structure | |
| CN108629140B (en) | Porous composite material structure design optimization method based on geodesic distance | |
| CN110598834A (en) | Binocular vision detection system structure optimization method | |
| CN115358009A (en) | Method for generating NURBS curve surface fairing in shell molding of vehicle | |
| CN111814383B (en) | B-spline density method-based self-supporting structure topology optimization design method | |
| CN112942837B (en) | Cantilever structure concrete 3D printing method and system | |
| Kimmel et al. | Finding the shortest paths on surfaces by fast global approximation and precise local refinement | |
| CN109299499B (en) | Multi-step structure optimization design method considering correction factors and aircraft | |
| CN112818488B (en) | Geometric-dimensional collaborative optimization design method for structural reinforcing rib distribution | |
| CN116422905A (en) | Laser scanning method, device and storage medium for powder bed fusion process | |
| CN107194998B (en) | A kind of method of multi-layer three-dimension point cloud single-layered | |
| CN115081113A (en) | Grid reinforcing rib structure optimization method based on multi-level set function parameterization | |
| CN108920785A (en) | A kind of monomeric complex model regularization modeling method and system |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |