Topological optimization method and device for large-scale fiber reinforced composite material structure

Abstract

The invention relates to the technical field of structure optimization, in particular to a topological optimization method and device for a large-scale fiber reinforced composite material structure, and the method comprises the steps: obtaining material attribute parameters and the like; calculating a unit stiffness matrix corresponding to each group of sampling angles according to a preset sampling interval and attribute parameters of a preset fiber reinforced composite material, calculating a template stiffness matrix of each grid level, and calculating the current displacement of each node of the densest grid level by adopting a geometric multi-grid preprocessing conjugate gradient method; and calculating a target function according to the current design variable and the current displacement of each node of the densest grid level, so as to update the design variable until a preset number of iterations is reached, and outputting a new optimization design. Therefore, the problems that due to the fact that the memory of a GPU global memory is limited, the element stiffness matrix cannot be stored in the topological optimization problem of the large-scale fiber reinforced composite material structure, and the calculation efficiency of the dynamic calculation element stiffness matrix in the finite element solving process is low are solved.

Description

Technical Field

[0001] The present invention relates to the technical field of structural optimization, and in particular to a topology optimization method and device for a large-scale fiber-reinforced composite material structure. Background Art

[0002] In recent years, the rapid development of 3D printing technology has made it possible to fabricate three-dimensional fiber-reinforced composite structures with continuously varying fiber orientations. This technology leverages the advantages of fiber-reinforced composites, such as high specific stiffness, high specific strength, and high designability, while also posing new challenges to the structural design of fiber-reinforced composites. Topology optimization, a multidisciplinary design method that integrates mathematics, computational mechanics, and computer science, offers new insights into the design of novel fiber-reinforced composite structures.

[0003] With the increasing demands of industrial production and computer computing power, topology optimization methods that are only suitable for solving small-scale two-dimensional problems are no longer able to meet the actual needs of production. The scale of new topology optimization problems is gradually increasing to tens of millions or even hundreds of millions of units, which undoubtedly poses a huge challenge to computer storage and computing resources. GPUs, as high-performance processors, are widely used in large-scale parallel computing problems and have achieved remarkable results in solving large-scale topology optimization problems of isotropic materials.

[0004] However, due to the limitations of GPU storage capacity and the speed of communication with computer memory, traditional large-scale topology optimization methods for isotropic materials cannot be directly applied to the topology optimization design of large-scale fiber-reinforced composite materials. The main reason is that in the optimization problem of isotropic material topology, such as metal, each unit has the same unit stiffness matrix, but the constitutive relationship matrix of each unit of the fiber-reinforced composite material is affected by the angle design variable. As the angle variable changes during the iteration process, the unit stiffness matrix will be different. Storing all the unit stiffness matrices consumes a huge amount of GPU video memory, which greatly limits the computational scale of the fiber-reinforced composite material topology optimization problem. If the unit stiffness matrix is ​​not stored but is dynamically calculated when needed, there will be a large number of repeated cumulative calculations, which will greatly reduce the performance of the optimization algorithm and greatly waste time costs. Summary of the Invention

[0005] The present invention provides a topology optimization method and device for large-scale fiber-reinforced composite materials, so as to solve the problems of being unable to store the stiffness matrices of each unit in the topology optimization problem of large-scale fiber-reinforced composite materials due to the limited global memory of GPU, and the low efficiency of dynamically calculating the stiffness matrix of each unit during the finite element solution process.

[0006] The first aspect of the present invention provides a topology optimization method for a large-scale fiber-reinforced composite material structure, comprising the following steps: obtaining property parameters of a preset fiber-reinforced composite material; initializing design variables of a pre-built topology optimization model, and defining constraints of the topology optimization model, wherein the design variables include density design variables and angle design variables; determining angle variable sampling points according to a preset sampling spacing, so as to calculate the unit stiffness matrix corresponding to each group of sampling angles according to the property parameters of the preset fiber-reinforced composite material and the angle variable sampling points, and storing the unit stiffness matrix corresponding to each group of sampling angles in a preset GPU global memory; performing an optimization iteration process based on the topology optimization model, filtering the density design variables and the angle design variables to obtain filtered density design variables and filtered angle design variables; based on an inter-unit parallel method, calculating the unit stiffness matrix corresponding to each group of sampling angles according to the preset fiber-reinforced composite material and the angle variable sampling points, and storing the unit stiffness matrix corresponding to each group of sampling angles in a preset GPU global memory. The template stiffness matrix of each node of each grid level is calculated based on the corresponding unit stiffness matrix, the template stiffness matrix solution format and the template stiffness matrix transfer relationship between each grid level, and the template stiffness matrices of all other levels except the densest grid level are stored in the preset GPU global memory; the smoothed residual after preprocessing of the densest grid level is calculated according to the template stiffness matrix of each node of each grid level, and the smoothed residual after preprocessing of the densest grid level is processed based on the preprocessing conjugate gradient method to obtain the current displacement of each node of the densest grid level; the objective function is solved according to the current design variables and the current displacement of each node of the densest grid level, and the sensitivity analysis of the objective function value and the constraint conditions is performed to obtain the sensitivity value; the current design variables are updated according to the objective function value and the sensitivity value to obtain a new optimized design, until the preset number of iterations is reached, and the new optimized design is output.

[0007] Optionally, determining the angle variable sampling points according to the preset sampling interval, so as to calculate the unit stiffness matrix corresponding to each group of sampling angles according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, includes:

[0008] Determine the angle variable sampling points within the value range of each angle design variable by using the preset sampling interval;

[0009] Constructing a constitutive relationship matrix of the orthotropic material according to the preset property parameters of the fiber-reinforced composite material;

[0010] The element stiffness matrix corresponding to each group of sampling angles is calculated according to the constitutive relationship matrix of the orthotropic material and the angle variable sampling points.

[0011] Optionally, the single stiffness matrix corresponding to each set of sampling angles is:

[0012]

[0013] Among them, K s (α s ,θ s ) is the element stiffness matrix corresponding to each set of sampling angles, D s (α s ,θ s ) is the material constitutive relationship matrix corresponding to the sampling angle, B is the shape function matrix composed of 8 sub-matrices, Ω e is the integration area represented by each cell.

[0014] Optionally, filtering the density design variable and the angle design variable to obtain filtered density design variable and filtered angle design variable includes:

[0015] The density design variable is filtered using a preset filtering method to obtain the filtered density design variable

[0016] The angle design variable α in the spherical coordinate system is based on trigonometric functions e and θ e Convert to a vector in Cartesian coordinate system

[0017] According to the preset fiber direction filter radius, the vector in the Cartesian coordinate system is respectively Filter to get the filtered direction vector

[0018] The filtered direction vector Perform normalization and transform the normalized direction vector according to the inverse trigonometric function The filtered angle design variables are obtained.

[0019] Optionally, the template stiffness matrix transfer relationship between each grid level is:

[0020]

[0021] in, is the template stiffness matrix of each node on a coarser grid level, i s is the sum index corresponding to the x direction, j s is the sum index corresponding to the y direction, k s is the corresponding sum index in the z direction, ω j =(2-|j1|)(2-|j2|)(2-|j3|) / 8,ω k=(2-|k1|)(2-|k2|)(2-|k3|) / 8, which are the weights of template stiffness matrix transfer between different grid levels. is the template stiffness matrix of each node on the denser mesh level.

[0022] Optionally, the calculating the smoothed residual after preprocessing of the densest grid level according to the template stiffness matrix of each node of each grid level, and processing the smoothed residual after preprocessing of the densest grid level based on the preprocessed conjugate gradient method to obtain the current displacement of each node of the densest grid level includes:

[0023] Based on the inter-element parallel method and the "V" cycle algorithm of geometric multigrid, the initial residual of the densest grid level is calculated according to the template stiffness matrix of the densest grid level, the preset initial displacement solution and the preset load of each node;

[0024] The initial residual of the densest grid level is transferred to the coarsest grid level by a constraint method to obtain the constrained residual of the coarsest grid level, and the constrained displacement of each node of the coarsest grid level is solved by a direct method according to the global stiffness matrix of the coarsest grid level and the constrained residual of the coarsest grid level;

[0025] Transferring the restricted displacement of each node of the coarsest grid level to the densest grid level by a continuation method to obtain the continuation displacement of each node of the densest grid level, and calculating the continuation residual of the densest grid level according to the template stiffness matrix of the densest grid level and the continuation displacement of each node of the densest grid level;

[0026] Processing the extended residual of the densest grid level by a post-smoothing method to obtain a smoothed residual after preprocessing of the densest grid level;

[0027] The pre-processed conjugate gradient method is adopted to obtain the current displacement of each node of the densest grid level according to the smoothed residual of the densest grid level.

[0028] Optionally, updating the current design variables according to the objective function value and the sensitivity value to obtain a new optimized design includes:

[0029] Based on the optimization criterion method, the density design variable in the current design variable is updated according to the objective function value and the sensitivity value to obtain a new density design variable;

[0030] Based on the moving asymptote method, the angle design variables in the current design variables are updated according to the objective function value and the sensitivity value to obtain new angle design variables.

[0031] A second embodiment of the present invention provides a topology optimization device for a large-scale fiber-reinforced composite material structure, comprising:

[0032] An acquisition module, used to obtain property parameters of a preset fiber-reinforced composite material;

[0033] An initialization module, configured to initialize design variables of a pre-built topology optimization model and define constraints of the topology optimization model, wherein the design variables include density design variables and angle design variables;

[0034] a first calculation module, configured to determine angle variable sampling points according to a preset sampling interval, calculate a unit stiffness matrix corresponding to each set of sampling angles according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, and store the unit stiffness matrix corresponding to each set of sampling angles in a preset GPU global memory;

[0035] a filtering module, configured to perform an optimization iteration process based on the topology optimization model, filter the density design variables and the angle design variables, and obtain filtered density design variables and filtered angle design variables;

[0036] A second calculation module is configured to calculate the template stiffness matrix of each node of each grid level based on the unit stiffness matrix corresponding to each group of sampling angles, the template stiffness matrix solution format, and the template stiffness matrix transfer relationship between each grid level in an inter-unit parallel manner, and store the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory;

[0037] a third solving module, configured to calculate a smoothed residual after preprocessing of the densest grid level according to the template stiffness matrix of each node of each grid level, and process the smoothed residual after preprocessing of the densest grid level based on a preprocessed conjugate gradient method to obtain a current displacement of each node of the densest grid level;

[0038] A solution and analysis module, configured to solve an objective function based on current design variables and current displacements of nodes at the densest grid level, and perform sensitivity analysis on the objective function value and the constraint conditions to obtain a sensitivity value;

[0039] The iterative optimization module is used to update the current design variables according to the objective function value and the sensitivity value to obtain a new optimized design until a preset number of iterations is reached, and output the new optimized design.

[0040] A third aspect of the present invention provides an electronic device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the topology optimization method for large-scale fiber-reinforced composite structures as described in the above embodiments.

[0041] A fourth aspect of the present invention provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the above-mentioned topology optimization method for large-scale fiber-reinforced composite material structures.

[0042] The economic device for the topology optimization method of large-scale fiber-reinforced composite materials proposed in an embodiment of the present invention adopts a method of pre-calculating and storing the unit stiffness matrix corresponding to the sampling angle to solve the problem of insufficient GPU video memory. At the same time, a method of rounding each unit angle to the sampling point to read the corresponding unit cell stiffness matrix is ​​adopted to improve the computational performance of the optimization algorithm.

[0043] Additional aspects and advantages of the present invention will be set forth in part in the description which follows and, in part, will be obvious from the description which follows, or may be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments in conjunction with the accompanying drawings, in which:

[0045] Figure 1 A flowchart of a topology optimization method for a large-scale fiber-reinforced composite material structure provided by an embodiment of the present invention;

[0046] Figure 2 A topological schematic diagram of a topological optimization method for a large-scale fiber-reinforced composite material structure provided by an embodiment of the present invention;

[0047] Figure 3 A schematic diagram of a method for defining angle variables provided in an embodiment of the present invention;

[0048] Figure 4 A flow chart of the "V" cycle algorithm for the preprocessing process of the geometric multigrid method provided in an embodiment of the present invention;

[0049] Figure 5 A schematic diagram of the working condition of the topology optimization problem provided by an embodiment of the present invention;

[0050] Figure 6 Iteration history change diagram provided by an embodiment of the present invention;

[0051] Figure 7 A schematic diagram of the optimization results provided by an embodiment of the present invention;

[0052] Figure 8 A schematic diagram showing a time comparison between an embodiment of the present invention and a traditional Gaussian integration method;

[0053] Figure 9 A schematic block diagram of a topology optimization device for a large-scale fiber-reinforced composite material structure provided by an embodiment of the present invention;

[0054] Figure 10 This is a schematic structural diagram of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0055] The following describes embodiments of the present invention in detail, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and are not to be construed as limiting the present invention.

[0056] The following describes a method and apparatus for topology optimization of a large-scale fiber-reinforced composite material structure according to an embodiment of the present invention with reference to the accompanying drawings.

[0057] Specifically, Figure 1 A schematic flow chart of a topology optimization method for a large-scale fiber-reinforced composite material structure provided by an embodiment of the present invention.

[0058] like Figure 1 As shown, the topology optimization method of the large-scale fiber reinforced composite structure includes the following steps:

[0059] In step S101, property parameters of a preset fiber-reinforced composite material are obtained.

[0060] Specifically, the properties of the fiber-reinforced composite material (transversely isotropic material), the size of the design domain (which may contain non-design domains), the voxelized model resolution, the geometric multigrid level, the structural topology optimization parameters, including but not limited to the filter radius of the density and angle variables, the sampling spacing of the angle variables, the volume fraction, etc. are input into the host memory and copied to the GPU global memory.

[0061] In step S102 , the design variables of the pre-built topology optimization model are initialized, and the constraints of the topology optimization model are defined, wherein the design variables include density design variables and angle design variables.

[0062] Specifically, the density and angle design variables of each unit in the pre-built topology optimization model, as well as the variables required for numerical processing (such as design variable filtering, sensitivity filtering, normalization, etc.) are declared and initialized in the GPU global memory, and their boundary conditions and loads are defined.

[0063] In step S103, the angle variable sampling points are determined according to the preset sampling interval, so as to calculate the unit stiffness matrix corresponding to each group of sampling angles according to the preset property parameters of the fiber reinforced composite material and the angle variable sampling points, and the unit stiffness matrix corresponding to each group of sampling angles is stored in the preset GPU global memory.

[0064] Specifically, if Figure 2 As shown, the sampling points of the angle variables are determined within the range of each angle design variable using a preset sampling interval; the constitutive relationship matrix of the orthotropic material is constructed based on the preset property parameters of the fiber-reinforced composite material; the element stiffness matrix corresponding to each set of sampling angles is calculated based on the constitutive relationship matrix of the orthotropic material and the sampling points of the angle variables, and the element stiffness matrix corresponding to each set of sampling angles is stored in the preset GPU global memory. The preset property parameters of the fiber-reinforced composite material include the longitudinal elastic modulus, transverse elastic modulus, Poisson's ratio, and in-plane shear modulus of the fiber-reinforced composite material.

[0065] Among them, the single stiffness matrix corresponding to each set of sampling angles is:

[0066]

[0067] Among them, K s (α s ,θ s ) is the element stiffness matrix corresponding to each set of sampling angles, D s (α s ,θ s ) is the material constitutive relationship matrix corresponding to the sampling angle, B is the shape function matrix composed of 8 sub-matrices, Ω e is the integration area represented by each cell.

[0068] In step S104, an optimization iteration process is performed based on the topology optimization model, and the density design variables and the angle design variables are filtered to obtain filtered density design variables and filtered angle design variables.

[0069] Specifically, the iterative process of topology optimization begins at this time, and in order to avoid numerical problems such as checkerboard phenomena, the density design variables and angle design variables are filtered.

[0070] Among them, density filtering is expressed as follows:

[0071] Use the preset filtering method to filter the density design variables and obtain the filtered density design variables Among them, the preset filtering method is:

[0072]

[0073] in, is the density design variable after filtering, ω(x k ) is the weight function, ω(x k )=max(R-||x k -x e ||,0),Φ e is the set of all units whose center point to the center point of the e-th unit is not greater than the filter radius, V k is the volume of the kth unit whose center point is not greater than the filter radius from the center point of the eth unit, ρ k is the density design variable of the kth unit whose center point is not greater than the center point of the eth unit, N e is the number of units;

[0074] Design variable α based on the angle in the spherical coordinate system according to trigonometric functions e and θ e Convert to a vector in Cartesian coordinate system Among them, the trigonometric functions are: and

[0075] According to the preset fiber direction filter radius, the vector in the Cartesian coordinate system is Filter to get the filtered direction vector

[0076] The filtered direction vector Perform normalization and transform the normalized direction vector according to the inverse trigonometric function The filtered angle design variables are obtained, where the inverse trigonometric function is: and

[0077] In step S105, based on the inter-unit parallel method, the template stiffness matrix of each node of each grid level is calculated according to the unit stiffness matrix corresponding to each group of sampling angles, the template stiffness matrix solution format and the template stiffness matrix transfer relationship between each grid level, and the template stiffness matrices of all other levels except the densest grid level are stored in the preset GPU global memory.

[0078] Specifically, the calculation process of the template stiffness matrix of each level required by the assembly-free geometric multi-grid preprocessing method is as follows: based on the parallel method between nodes, the calculation process of the template matrix corresponding to the 27 adjacent nodes of each node is assigned to a separate GPU thread to complete. The 27 nodes belong to the 8 cells adjacent to the central node. When the stiffness matrix of a certain adjacent cell is used to calculate its contribution to the template matrix of a certain node, the filtered angle variable in the cell will be rounded to the sampling point angle value closest to it. Through the index corresponding to the sampling point angle, the GPU global memory is accessed to obtain the cell stiffness matrix corresponding to the pre-calculated sampling point angle.

[0079] Among them, the template stiffness matrix transfer relationship between each grid level is:

[0080]

[0081] in, is the template stiffness matrix of each node on a coarser grid level, i s is the sum index corresponding to the x direction, j s is the sum index corresponding to the y direction, k s is the corresponding sum index in the z direction, ω j =(2-|j1|)(2-|j2|)(2-|j3|) / 8,ω k =(2-|k1|)(2-|k2|)(2-|k3|) / 8, which are the weights of template stiffness matrix transfer between different grid levels. is the template stiffness matrix of each node on the denser mesh level.

[0082] In step S106, the smoothed residual after preprocessing of the densest grid level is calculated according to the template stiffness matrix of each node of each grid level, and the smoothed residual after preprocessing of the densest grid level is processed based on the preprocessing conjugate gradient method to obtain the current displacement of each node of the densest grid level.

[0083] Specifically, if Figure 4 As shown, the smoothed residual after the densest grid level preprocessing is calculated according to the template stiffness matrix of each node at each grid level, and the smoothed residual after the densest grid level preprocessing is processed based on the preprocessed conjugate gradient method to obtain the current displacement of each node at the densest grid level, including:

[0084] Based on the inter-element parallel method and the "V" cycle algorithm of geometric multigrid, the initial residual of the densest mesh level is calculated according to the template stiffness matrix of the densest mesh level, the preset initial displacement solution and the preset load of each node;

[0085] The constraint method is used to transfer the initial residual of the densest grid level to the coarsest grid level to obtain the constrained residual of the coarsest grid level. The direct method is then used to solve the constrained displacement of each node of the coarsest grid level based on the global stiffness matrix of the coarsest grid level and the constrained residual of the coarsest grid level.

[0086] The restricted displacement of each node of the coarsest mesh level is transferred to the densest mesh level through the continuation method to obtain the continuation displacement of each node of the densest mesh level, and the continuation residual of the densest mesh level is calculated based on the template stiffness matrix of the densest mesh level and the continuation displacement of each node of the densest mesh level;

[0087] The post-smoothing method is used to process the extended residual of the densest grid level to obtain the smoothed residual after the densest grid level preprocessing;

[0088] The preconditioned conjugate gradient method is used to improve the convergence speed. The current displacement of each node at the densest grid level is obtained according to the smoothed residual of the densest grid level.

[0089] The smoothing process adopts the Jacobi method or the Gauss-Seidel method.

[0090] In step S107 , the objective function is solved according to the current design variables and the current displacement of each node at the densest grid level, and a sensitivity analysis is performed on the objective function value and the constraint conditions to obtain a sensitivity value.

[0091] Specifically, based on the current design variables and the current displacement of each node at the densest grid level, the objective function value of the optimized structure, i.e., the structural strain energy, is solved; and a sensitivity analysis of the objective function and each constraint condition is performed to obtain the sensitivity value of each design variable.

[0092] In step S108, the current design variables are updated according to the objective function value and the sensitivity value to obtain a new optimized design, until a preset number of iterations is reached, and the new optimized design is output.

[0093] Specifically, based on the optimization criterion method, the density design variable is updated according to the objective function value and the sensitivity value to obtain a new density design variable. Based on the moving asymptote method, the angle design variable is updated according to the objective function value and the sensitivity value to obtain the angle design variable, and it is judged whether the current number of optimization iterations reaches the preset threshold. If so, the iteration is terminated and the new optimized design is output. Otherwise, the iteration continues and steps S104-S108 are executed.

[0094] like Figure 5As shown below, the three-dimensional cantilever beam structure obtained after topology optimization is used as an example to illustrate the finite element analysis of large-scale fiber reinforced composite materials. In this structure, the area on the left end face of the cube initial design domain close to 1 / 4 of the width on both sides is a fixed constraint, with a height of 1.0m, a length of 2.0m, and a width of 1.0m. The number of units after the design domain is discretized into a finite element model is N e It can be expressed as: N e =128×256×128=4.1943×10 6 The longitudinal elastic modulus E2 of fiber reinforced composite materials is 132×10 9 pa, transverse elastic modulus E1 is 9.6×10 9 pa, Poisson's ratio v 21 is 0.258, v 13 is 0.35, the in-plane shear modulus G 21 5.8×10 9 The structure is subjected to a concentrated load of 20 kN in the vertical direction at the midpoint of the lower end of the right end face. The optimization goal is to minimize the structural strain energy and impose a volume constraint with a volume fraction of 10%.

[0095] Step 1: Establish the mathematical model of the topology optimization problem as follows:

[0096] find: ρ, α, θ

[0097] min:c(ρ, α, θ)=F T U

[0098] st:K(ρ,α,θ)U=F

[0099] :G k (ρ,α,θ)≤0(k=1,…,N c )

[0100] :0≤ρe≤1(e=1,…,N e )

[0101] :-π≤α e ≤π(e=1,…,N e )

[0102] :-π / 2≤θ e ≤π / 2(e=1,…,N e )

[0103] Where c is the structural strain energy, K is the overall stiffness matrix, U is the node displacement vector, and F is the overall structural load vector. e is the number of units, ρ e is the density design variable, α e and θ eThey are the two angle design variables of each unit, which are defined as follows Figure 3 As shown. Volume constraint function G k It can be expressed as follows:

[0104]

[0105] Among them, v e is the volume of each unit, and ξ is the volume fraction.

[0106] Step 2: Initialize design variables The sampling spacing parameter gap = 1°, and the density filter radius and angle filter radius are both set to 1.5 times the cell size.

[0107] Step 3: According to the sampling interval, determine the sampling points of the angle variables within the value range of each angle variable. In this embodiment, 1° is selected as the angle sampling interval, and the sampling angle α is s is an integer between [-180°, 180°], the number of sampling points is 361, and the sampling angle θ s is an integer between [-90°, 90°], and the number of sampling points is 181. The element stiffness matrix corresponding to each set of sampling angles is calculated as follows:

[0108]

[0109] Among them, B=[B0B1B2B3B4B5B6B7] is a shape function matrix composed of 8 sub-matrices.

[0110] Among them, the submatrix B ei It can be expressed as:

[0111]

[0112] Where K s (a s ,θ s ) is the element stiffness matrix corresponding to each set of sampling angles, D s (α s ,θ s ) is the material constitutive relationship matrix corresponding to the sampling angle, B is the shape function matrix composed of 8 sub-matrices, Ω e is the integration area represented by each cell. is the material constitutive relationship matrix corresponding to the sampling angle, which is calculated as follows:

[0113] Where C is the constitutive relation matrix of the orthotropic material, which is expressed by the above material property parameters as follows:

[0114]

[0115] The rotation matrix of the angle variable can be expressed as:

[0116]

[0117] Among them, c α =cos(α s ), s α =sin(α s ), c θ =cos(θ s ), s θ = sin(θ s ).

[0118] Step 4: To avoid numerical problems such as checkerboard phenomenon, the density design variables and angle design variables are filtered. Density filtering is expressed as follows:

[0119]

[0120] in, is the density design variable after filtering, ω(x k ) is the weight function, ω(x k )=max(R-||x k -x e ||, 0), Φ e is the set of all units whose center point to the center point of the e-th unit is not greater than the filter radius, V k is the volume of the kth unit whose center point is not greater than the filter radius from the center point of the eth unit, ρ k is the density design variable of the kth unit whose center point is not greater than the center point of the eth unit, N e is the number of units.

[0121] The operation process of Angular variable filtering is as follows:

[0122] (1) According to the following trigonometric transformation: and The angle design variable α in the spherical coordinate system e and θ e Convert to a vector in Cartesian coordinate system

[0123] (2) According to the preset fiber direction filter radius, Filtered fiber direction vector

[0124] (3) Normalize the filtered direction vector to obtain

[0125] (4) And according to the inverse trigonometric function relationship: and Calculate the filtered angle variable.

[0126] Step five: Calculate the template matrices of each level required for the assembly-free geometric multigrid preprocessing method. The calculation process adopts a node-to-node parallel approach. The calculation process of the template matrices corresponding to the 27 adjacent nodes of each node is assigned to a separate GPU thread to complete. The 27 nodes belong to the 8 cells adjacent to the central node. When the stiffness matrix of a certain adjacent cell is used to calculate its contribution to the template matrix of a certain node, the filtered angle variable in the cell will be rounded to the sampling point angle value closest to it. Through the index corresponding to the sampling point angle, the GPU global memory is accessed to obtain the cell stiffness matrix corresponding to the pre-calculated sampling point angle as an approximate solution to the actual cell stiffness matrix.

[0127] Among them, due to the use of SIMP method, the approximate solution of the element stiffness matrix can be expressed as:

[0128]

[0129] Where, is the actual element stiffness matrix corresponding to the actual design variables, is the density penalty function Usually take ∈=10 -6 and p = 3, K s (α s ,θ s ) is the sampling unit stiffness matrix corresponding to the sampling angle obtained by mapping the actual angle variable, α s and θ s are the actual angle variables after filtering and The mapped sampling angle.

[0130] Step 6: Use Figure 4 The "V-loop" algorithm shown in the figure performs geometric multigrid preprocessing, and its process includes:

[0131] (1) Calculate the preprocessing matrices at each level required for the assembly-free geometric multigrid method and store them in the GPU global memory;

[0132] (2) Smoothing the densest level;

[0133] (3) Calculate the residual of the static equilibrium equation of the densest level grid;

[0134] (4) The residual is transferred to the next densest level as the right-hand side of the static equilibrium equation through the restriction method, and the above process is repeated until the coarsest grid level;

[0135] (5) At the coarsest level, assemble the global stiffness matrix and solve the static equilibrium equations using a direct method based on Cholesky decomposition;

[0136] (6) Correct the displacement value of the sub-coarse level by the continuation method, and repeat the steps until the displacement solution of the densest level is obtained;

[0137] (7) Smoothing the densest level.

[0138] The calculation process of process (5) in step 6 is completed by CPU calculation, and the displacement data of each node is copied to the GPU global memory. All other processes are completed by GPU parallel calculation.

[0139] Step 7: Use the preconditioned conjugate gradient method to solve the static equilibrium equation and calculate the displacement of each node. The convergence tolerance of the conjugate gradient method iterative solver is set to 1×10 -5 .

[0140] Step 8: Based on the stiffness matrix of each unit calculated in step 3 and the displacement of each node calculated in step 7, the objective function of the optimized structure, i.e., the structural strain energy, is solved; and a sensitivity analysis of the objective function and each constraint condition is performed;

[0141] Step 9: Optimize and update the design variables: update the density design variables using the optimization criterion method, and update the angle design variables using the moving asymptote method to obtain a new optimized design;

[0142] Step 10, such as Figure 6 and 7 As shown, determine whether the iteration converges: determine whether the number of optimization iterations reaches 150. If it reaches 150, the convergence condition is met, the iteration is terminated, and the optimization result is output; if the convergence condition is not met, jump to step 4) to continue the iteration.

[0143] The following is an optimization design of a cantilever beam structure using the method of the embodiment of the present invention and the traditional large-scale topology optimization method of online calculation of element stiffness matrix. The embodiment of the present invention is run on a desktop computer equipped with Intel(R) Core(R) i7-13700KF@5.40GHz and NVIDIA RTX 4070GPU (12GB video memory and 5888 CUDA cores). The topology optimization results are shown in Figure 2. Figure 8 The running time statistics of the optimization algorithm are shown in Table 1.

[0144] Table 1 Comparison of results between the method of the present invention and the traditional method

[0145] Method of the present invention Traditional methods Number of grids 128×256×128 128×256×128 Runtime 5176.5s 76917.3s Initial design strain energy 15246.2 15246.2 Final strain energy (150 iterations) 3.36 3.34

[0146] Results demonstrate that the embodiments of the present invention can simultaneously optimize the density and fiber orientation of large-scale fiber-reinforced composite materials, achieving a clear topological structure and a continuously varying fiber orientation design. The objective function values ​​of the optimization results are nearly identical to those of traditional methods. Furthermore, the proposed optimization algorithm significantly improves computational efficiency compared to traditional methods, significantly reducing the time cost of solving large-scale fiber-reinforced composite topology optimization problems. While ensuring computational accuracy, it achieves a reasonable balance between the storage space cost of the anisotropic material unit stiffness matrix and the computational time cost, demonstrating its significant practical significance.

[0147] In summary, the GPU-accelerated topology optimization method for large-scale fiber-reinforced composite structures proposed in accordance with an embodiment of the present invention significantly reduces the consumption of GPU global memory by storing the fiber-reinforced composite unit stiffness matrix. Furthermore, the storage consumption of the unit stiffness matrix is ​​affected only by the number of angular variable sampling points and does not increase with increasing model resolution. This significantly reduces the time cost of iterative optimization. Within the range of angular variable values, when a sufficient number of equally spaced sampling points are selected, the objective function value optimized by this method is nearly identical to the exact value, meeting very high precision requirements.

[0148] Next, a topology optimization device for a large-scale fiber-reinforced composite material structure according to an embodiment of the present invention will be described with reference to the accompanying drawings.

[0149] Figure 9 4 is a block diagram of a topology optimization device for a large-scale fiber-reinforced composite material structure according to an embodiment of the present invention.

[0150] like Figure 9 As shown, the topology optimization device 90 for a large-scale fiber-reinforced composite material structure includes: an acquisition module 901, an initialization module 902, a first calculation module 903, a filtering module 904, a second calculation module 905, a third solution module 906, a solution and analysis module 907 and an iterative optimization module 908.

[0151] The acquisition module 901 is configured to obtain the preset property parameters of the fiber-reinforced composite material. The initialization module 902 is configured to initialize the design variables of a pre-built topology optimization model and define the constraints of the topology optimization model, wherein the design variables include density design variables and angle design variables. The first calculation module 903 is configured to determine the angle variable sampling points based on a preset sampling interval, calculate the element stiffness matrix corresponding to each set of sampling angles based on the preset fiber-reinforced composite material property parameters and the angle variable sampling points, and store the element stiffness matrix corresponding to each set of sampling angles in the preset GPU global memory. The filtering module 904 is configured to perform an optimization iteration based on the topology optimization model, filter the density design variables and the angle design variables, and obtain filtered density design variables and filtered angle design variables. The second calculation module 905 is configured to calculate the template stiffness matrix for each node at each grid level in an inter-unit parallel manner based on the element stiffness matrix corresponding to each set of sampling angles, the template stiffness matrix solution format, and the template stiffness matrix transfer relationship between each grid level. The template stiffness matrices for all levels except the densest grid level are stored in the preset GPU global memory. The third solution module 906 is used to calculate the smoothed residual after preprocessing of the densest mesh level based on the template stiffness matrix of each node at each mesh level, and process the smoothed residual after preprocessing of the densest mesh level based on the preprocessed conjugate gradient method to obtain the current displacement of each node at the densest mesh level. The solution and analysis module 907 is used to solve the objective function based on the current design variables and the current displacement of each node at the densest mesh level, and perform sensitivity analysis on the objective function value and the constraint conditions to obtain the sensitivity value. The iterative optimization module 908 is used to update the current design variables based on the objective function value and the sensitivity value to obtain a new optimized design, until the preset number of iterations is reached, and then output the new optimized design.

[0152] It should be noted that the aforementioned explanation of the embodiment of the topology optimization method for a large-scale fiber-reinforced composite material structure is also applicable to the topology optimization device for a large-scale fiber-reinforced composite material structure of this embodiment, and will not be repeated here.

[0153] The GPU-accelerated topology optimization device for large-scale fiber-reinforced composite structures, proposed in accordance with an embodiment of the present invention, significantly reduces the GPU global memory consumption required to store the fiber-reinforced composite unit stiffness matrix. Furthermore, the unit stiffness matrix storage consumption is affected only by the number of sampling points for the angular variables and does not increase with increasing model resolution. This significantly reduces the time cost of iterative optimization. Within the range of angular variables, when a sufficient number of equally spaced sampling points are selected, the objective function value obtained by this method is nearly identical to the exact value, meeting very high precision requirements.

[0154] Figure 10This is a schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. The electronic device may include:

[0155] A memory 1001 , a processor 1002 , and a computer program stored in the memory 1001 and executable on the processor 1002 .

[0156] When the processor 1002 executes the program, the topology optimization method for a large-scale fiber-reinforced composite material structure provided in the above embodiment is implemented.

[0157] Furthermore, the electronic device further includes:

[0158] The communication interface 1003 is used for communication between the memory 1001 and the processor 1002 .

[0159] The memory 1001 is used to store computer programs that can be run on the processor 1002 .

[0160] The memory 1001 may include a high-speed RAM memory, and may also include a non-volatile memory (non-volatile memory), such as at least one disk memory.

[0161] If the memory 1001, processor 1002, and communication interface 1003 are implemented independently, the communication interface 1003, memory 1001, and processor 1002 can be connected to each other via a bus and communicate with each other. The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus. The bus can be divided into an address bus, a data bus, a control bus, etc. For ease of representation, Figure 10 Only one thick line is used in the diagram, but this does not mean that there is only one bus or one type of bus.

[0162] Optionally, in a specific implementation, if the memory 1001, the processor 1002 and the communication interface 1003 are integrated on a chip, the memory 1001, the processor 1002 and the communication interface 1003 can communicate with each other through an internal interface.

[0163] The processor 1002 may be a central processing unit (CPU), an application specific integrated circuit (ASIC), or one or more integrated circuits configured to implement the embodiments of the present invention.

[0164] An embodiment of the present invention further provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-mentioned topology optimization method for large-scale fiber-reinforced composite material structures.

[0165] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic expressions of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or N embodiments or examples in a suitable manner. In addition, those skilled in the art can combine and combine different embodiments or examples described in this specification and the features of different embodiments or examples without contradiction.

[0166] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be understood to indicate or imply relative importance or implicitly specify the number of technical features indicated. Thus, a feature specified as "first" or "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "N" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0167] Any process or method description in a flowchart or otherwise described herein may be understood to represent a module, segment or portion of code comprising one or N executable instructions for implementing a custom logical function or step of a process, and the scope of the preferred embodiments of the present invention includes alternative implementations in which functions may be performed out of the order shown or discussed, including performing functions in a substantially simultaneous manner or in the reverse order depending on the functions involved, which should be understood by those skilled in the art to which the embodiments of the present invention pertain.

[0168] The logic and / or steps represented in the flowcharts or otherwise described herein, for example, can be considered as a sequenced list of executable instructions for implementing the logical functions, and can be embodied in any computer-readable medium for use by, or in conjunction with, an instruction execution system, apparatus, or device (e.g., a computer-based system, a system including a processor, or other system that can fetch and execute instructions from an instruction execution system, apparatus, or device). For purposes of this specification, a "computer-readable medium" can be any device that can contain, store, communicate, propagate, or transport a program for use by, or in conjunction with, an instruction execution system, apparatus, or device. More specific examples (a non-exhaustive list) of computer-readable media include the following: an electrical connection with one or N wires (electronic devices), a portable computer disk cartridge (magnetic device), random access memory (RAM), read-only memory (ROM), erasable and programmable read-only memory (EPROM or flash memory), fiber optic devices, and a portable compact disc read-only memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program is printed, since the program can be obtained electronically by optically scanning the paper or other medium and then editing, interpreting or processing it in other suitable ways as necessary, and then storing it in a computer memory.

[0169] It should be understood that various parts of the present invention can be implemented using hardware, software, firmware, or a combination thereof. In the above embodiment, the N steps or methods can be implemented using software or firmware stored in a memory and executed by a suitable instruction execution system. If implemented using hardware, as in another embodiment, any one of the following technologies known in the art or a combination thereof can be used: a discrete logic circuit having logic gate circuits for implementing logic functions on data signals, an application-specific integrated circuit having suitable combinational logic gate circuits, a programmable gate array (PGA), a field programmable gate array (FPGA), etc.

[0170] Those skilled in the art will understand that all or part of the steps in the method of the above embodiment can be completed by instructing related hardware through a program, and the program can be stored in a computer-readable storage medium. When the program is executed, it includes one or a combination of the steps of the method embodiment.

[0171] In addition, the functional units in the various embodiments of the present invention may be integrated into a single processing module, or each unit may exist physically separately, or two or more units may be integrated into a single module. The aforementioned integrated modules may be implemented in the form of hardware or in the form of software functional modules. If the integrated modules are implemented in the form of software functional modules and sold or used as independent products, they may also be stored in a computer-readable storage medium.

[0172] The storage medium mentioned above may be a read-only memory, a magnetic disk, or an optical disk, etc. Although the embodiments of the present invention have been shown and described above, it should be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Persons skilled in the art may make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A topology optimization method for large-scale fiber-reinforced composite structures, characterized in that: The following steps are involved: Obtaining property parameters of preset fiber reinforced composite materials; Initializing design variables of a pre-built topology optimization model and defining constraints of the topology optimization model, wherein the design variables include density design variables and angle design variables; Determining angle variable sampling points according to a preset sampling interval, calculating a unit stiffness matrix corresponding to each set of sampling angles according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, and storing the unit stiffness matrix corresponding to each set of sampling angles in a preset GPU global memory; Performing an optimization iteration process based on the topology optimization model, filtering the density design variables and the angle design variables to obtain filtered density design variables and filtered angle design variables; Based on an inter-unit parallel approach, the template stiffness matrix of each node of each grid level is calculated according to the unit stiffness matrix corresponding to each group of sampling angles, the template stiffness matrix solution format, and the template stiffness matrix transfer relationship between each grid level, and the template stiffness matrices of all other levels except the densest grid level are stored in the preset GPU global memory; the smoothed residual after preprocessing of the densest grid level is calculated according to the template stiffness matrix of each node of each grid level, and the smoothed residual after preprocessing of the densest grid level is processed based on the preprocessed conjugate gradient method to obtain the current displacement of each node of the densest grid level; Solving the objective function according to the current design variables and the current displacement of each node of the densest grid level, and performing sensitivity analysis on the objective function value and the constraint condition to obtain a sensitivity value; The current design variables are updated according to the objective function value and the sensitivity value to obtain a new optimized design, until a preset number of iterations is reached, and the new optimized design is output.

2. The topology optimization method for large-scale fiber-reinforced composite material structures according to claim 1, characterized in that: The step of determining the angle variable sampling points according to the preset sampling interval, and calculating the element stiffness matrix corresponding to each group of sampling angles according to the preset property parameters of the fiber reinforced composite material and the angle variable sampling points, includes: Determine the angle variable sampling points within the value range of each angle design variable by using the preset sampling interval; Constructing a constitutive relationship matrix of the orthotropic material according to the preset property parameters of the fiber-reinforced composite material; The element stiffness matrix corresponding to each group of sampling angles is calculated according to the constitutive relationship matrix of the orthotropic material and the angle variable sampling points.

3. The topology optimization method for large-scale fiber-reinforced composite material structures according to claim 2, characterized in that: The single stiffness matrix corresponding to each set of sampling angles is: Among them, K s (α s ,θ s ) is the element stiffness matrix corresponding to each set of sampling angles, D s (α s ,θ s ) is the material constitutive relationship matrix corresponding to the sampling angle, B is the shape function matrix composed of 8 sub-matrices, Ω e is the integration area represented by each cell.

4. The topology optimization method of large-scale fiber-reinforced composite material structure according to claim 1, characterized in that: The filtering of the density design variable and the angle design variable to obtain filtered density design variable and filtered angle design variable includes: The density design variable is filtered using a preset filtering method to obtain the filtered density design variable The angle design variable α in the spherical coordinate system is based on trigonometric functions e and θ e Convert to a vector in Cartesian coordinate system According to the preset fiber direction filter radius, the vector in the Cartesian coordinate system is respectively Filter to get the filtered direction vector The filtered direction vector Perform normalization and transform the normalized direction vector according to the inverse trigonometric function The filtered angle design variables are obtained.

5. The topology optimization method of large-scale fiber-reinforced composite material structure according to claim 1, characterized in that: The template stiffness matrix transfer relationship between each grid level is: in, is the template stiffness matrix of each node on a coarser grid level, i s is the sum index corresponding to the x direction, j s is the sum index corresponding to the y direction, k s is the corresponding sum index in the z direction, ω j =(2-|j1|)(2-|j2|)(2-|j3|) / 8,ω k =(2-|k1|)(2-|k2|)(2-|k3|) / 8, which are the weights of template stiffness matrix transfer between different grid levels. is the template stiffness matrix of each node on the denser mesh level.

6. The topology optimization method of large-scale fiber-reinforced composite material structure according to claim 1, characterized in that: The calculation of the smoothed residual after the densest grid level preprocessing according to the template stiffness matrix of each node of each grid level, and processing the smoothed residual after the densest grid level preprocessing based on the preprocessing conjugate gradient method to obtain the current displacement of each node of the densest grid level includes: Based on the inter-element parallelism and the "V" cycle algorithm of geometric multigrid, the initial residual of the densest grid level is calculated according to the template stiffness matrix of the densest grid level, the preset initial displacement solution and the preset load of each node; The initial residual of the densest grid level is transferred to the coarsest grid level by a constraint method to obtain the constrained residual of the coarsest grid level, and the constrained displacement of each node of the coarsest grid level is solved by a direct method according to the global stiffness matrix of the coarsest grid level and the constrained residual of the coarsest grid level; Transferring the restricted displacement of each node of the coarsest grid level to the densest grid level by a continuation method to obtain the continuation displacement of each node of the densest grid level, and calculating the continuation residual of the densest grid level according to the template stiffness matrix of the densest grid level and the continuation displacement of each node of the densest grid level; Processing the extended residual of the densest grid level by a post-smoothing method to obtain a smoothed residual after preprocessing of the densest grid level; The pre-processed conjugate gradient method is adopted to obtain the current displacement of each node of the densest grid level according to the smoothed residual of the densest grid level.

7. The topology optimization method of large-scale fiber reinforced composite material structure according to claim 1, characterized in that: The updating of the current design variables according to the objective function value and the sensitivity value to obtain a new optimized design includes: Based on the optimization criterion method, the density design variable in the current design variable is updated according to the objective function value and the sensitivity value to obtain a new density design variable; Based on the moving asymptote method, the angle design variables in the current design variables are updated according to the objective function value and the sensitivity value to obtain new angle design variables.

8. A topology optimization device for large-scale fiber-reinforced composite materials, characterized in that: include: An acquisition module, used to obtain property parameters of a preset fiber-reinforced composite material; An initialization module, configured to initialize design variables of a pre-built topology optimization model and define constraints of the topology optimization model, wherein the design variables include density design variables and angle design variables; a first calculation module, configured to determine angle variable sampling points according to a preset sampling interval, calculate a unit stiffness matrix corresponding to each set of sampling angles according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, and store the unit stiffness matrix corresponding to each set of sampling angles in a preset GPU global memory; a filtering module, configured to perform an optimization iteration process based on the topology optimization model, filter the density design variables and the angle design variables, and obtain filtered density design variables and filtered angle design variables; A second calculation module is configured to calculate the template stiffness matrix of each node of each grid level based on the unit stiffness matrix corresponding to each group of sampling angles, the template stiffness matrix solution format, and the template stiffness matrix transfer relationship between each grid level in an inter-unit parallel manner, and store the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory; a third solving module, configured to calculate a smoothed residual after preprocessing of the densest grid level according to the template stiffness matrix of each node of each grid level, and process the smoothed residual after preprocessing of the densest grid level based on a preprocessed conjugate gradient method to obtain a current displacement of each node of the densest grid level; A solution and analysis module, configured to solve an objective function based on current design variables and current displacements of nodes at the densest grid level, and perform sensitivity analysis on the objective function value and the constraint conditions to obtain a sensitivity value; The iterative optimization module is used to update the current design variables according to the objective function value and the sensitivity value to obtain a new optimized design until a preset number of iterations is reached, and output the new optimized design.

9. An electronic device, characterized in that: include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the topology optimization method for a large-scale fiber-reinforced composite material structure according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: The program is executed by a processor to implement the topology optimization method for a large-scale fiber-reinforced composite material structure according to any one of claims 1 to 7.

Images