Fiber reinforced composite material structure topological optimization method based on Taylor expansion

Abstract

The invention relates to the technical field of structure optimization, in particular to a fiber reinforced composite material structure topological optimization method based on Taylor expansion, which comprises the following steps: acquiring material attribute parameters and the like; calculating an element stiffness matrix corresponding to each group of sampling angles and a partial derivative matrix of the element stiffness matrix to the sampling angles according to the material attribute parameters and a preset sampling interval; respectively filtering the design variables; a geometric multi-grid preprocessing conjugate gradient method iteration solver is adopted to calculate the current displacement of each node of the densest grid level; 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, each 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 method for topological optimization of fiber-reinforced composite material structures based on Taylor expansion. Background Art

[0002] In recent years, fiber-reinforced composites have become an important component of high-tech aerospace products such as aircraft, launch vehicles, and satellites. With the rapid development of 3D printing technology, it has become possible to process composite components with complex structures and continuous spatial changes in fiber orientation.

[0003] For large aerospace components, the characteristic dimensions of different dimensions in the overall structure differ by several orders of magnitude, resulting in the need for high-resolution meshing of the overall structure during finite element analysis to ensure the accuracy of the analysis results. This usually results in the scale of the finite element model reaching tens of millions or even billions of units, which far exceeds the computing power and applicability of current commercial software, and also poses a huge challenge to the storage and computing resources of ordinary computers.

[0004] As a high-performance processor, GPUs are widely used in large-scale parallel computing and have achieved remarkable results in solving large-scale finite element analysis problems for isotropic materials. However, for fiber-reinforced composite materials, the continuous spatial variation of the fiber orientation results in different element stiffness matrices in each finite element model. Storing each element stiffness matrix consumes a huge amount of GPU memory. If the element stiffness matrix is ​​not stored but calculated dynamically, a large number of cumulative calculations will be repeated, resulting in a significant reduction in the performance of the optimization algorithm. These problems make it difficult to solve high-resolution fiber-reinforced composite finite element problems using GPU-based parallel acceleration methods. Summary of the Invention

[0005] The present invention provides a fiber-reinforced composite material structure topology optimization method based on Taylor expansion, so as to solve the problems of being unable to store the stiffness matrix of each unit in the topology optimization problem of large-scale fiber-reinforced composite material structure due to the limited memory of GPU global memory, and the low efficiency of dynamic calculation of unit stiffness matrix during finite element solution.

[0006] A first embodiment of the present invention provides a method for topological optimization of a fiber-reinforced composite material structure based on Taylor expansion, comprising the following steps:

[0007] Obtaining property parameters of preset fiber reinforced composite materials;

[0008] 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;

[0009] Determining angle variable sampling points according to a preset sampling interval, calculating a unit stiffness matrix corresponding to each group of sampling angles and a partial derivative matrix of each unit stiffness matrix with respect to the 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 group of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angles in a preset GPU global memory;

[0010] 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;

[0011] Based on an inter-unit parallel approach, the unit stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angle are processed according to a Taylor expansion formula of a preset order to solve an approximate value of the unit stiffness matrix;

[0012] Calculating the template stiffness matrix of each node of each grid level according to the approximate value of the unit stiffness matrix, the template stiffness matrix solution format and the template stiffness matrix transfer relationship between each grid level, and storing the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory;

[0013] 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;

[0014] 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;

[0015] 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.

[0016] Optionally, 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 and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, includes:

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

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

[0019] The element stiffness matrix corresponding to each group of sampling angles and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle are calculated according to the constitutive relationship matrix of the orthotropic material and the angle variable sampling points.

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

[0021]

[0022] 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.

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

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

[0025] 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

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

[0027] 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.

[0028] Optionally, the SIMP method is used to solve the approximate value of the unit stiffness matrix, wherein the approximate value of the unit stiffness matrix is:

[0029]

[0030]

[0031] in, 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 e and θ e The mapped sampling angle.

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

[0033]

[0034] in, is the template stiffness matrix of each node on a coarser grid level, i s is the summation 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.

[0035] 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:

[0036] 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 each node of each grid level, the preset initial displacement solution and the preset load of each node;

[0037] 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;

[0038] 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;

[0039] 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;

[0040] 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.

[0041] A second embodiment of the present invention provides a fiber-reinforced composite material structure topology optimization device based on Taylor expansion, comprising:

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

[0043] 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;

[0044] 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 group of sampling angles and a partial derivative matrix of each unit stiffness matrix with respect to the 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 group of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angles in a preset GPU global memory;

[0045] 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;

[0046] A second calculation module is configured to process the unit stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angles according to a Taylor expansion formula of a preset order in an inter-unit parallel manner, and solve an approximate value of the unit stiffness matrix;

[0047] a third solving module, configured to calculate the template stiffness matrix of each node of each grid level according to the approximate value of the unit stiffness matrix, the template stiffness matrix solving format, and the template stiffness matrix transfer relationship between each grid level, and store the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory;

[0048] a fourth calculation 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;

[0049] 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;

[0050] 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.

[0051] 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 fiber-reinforced composite material structure topology optimization method based on Taylor expansion as described in the above embodiment.

[0052] 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 Taylor expansion-based fiber-reinforced composite material structure topology optimization method.

[0053] The fiber-reinforced composite material structure topology optimization method based on Taylor expansion 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 and the partial derivative matrix of the unit stiffness matrix with respect to different angle variables, so as to solve the problems of the inability to store the unit stiffness matrix in the topology optimization problem of large-scale fiber-reinforced composite material structure due to the limited memory of the GPU global memory, and the low efficiency of the dynamic calculation of the unit stiffness matrix during the finite element solution process; at the same time, the method of Taylor expansion of each unit stiffness matrix at the sampling point after the angle variable is rounded is adopted to calculate the approximate unit stiffness matrix, thereby improving the computational efficiency of the finite element analysis and topology optimization algorithm; the Taylor expansion order can also be flexibly adjusted according to the need for the fitting accuracy of the unit stiffness matrix to balance the contradiction between the storage consumption of the GPU global memory and the computing performance. Usually, a first-order Taylor expansion can meet the accuracy requirements of the computational accuracy of general engineering problems. After optimization, a clear topological structure and continuously changing fiber direction can be obtained. This method can obtain an objective function value that is almost consistent with the accurate solution of the unit stiffness matrix.

[0054] 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

[0055] 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:

[0056] Figure 1 A flowchart of a method for topology optimization of a fiber-reinforced composite material structure based on Taylor expansion provided in an embodiment of the present invention;

[0057] Figure 2 A topological schematic diagram of a fiber-reinforced composite material structure topology optimization method based on Taylor expansion provided by an embodiment of the present invention;

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

[0059] 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;

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

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

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

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

[0064] Figure 9 A schematic block diagram of a fiber-reinforced composite material structure topology optimization device based on Taylor expansion provided in an embodiment of the present invention;

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

[0066] 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.

[0067] The following describes a fiber-reinforced composite material structure topology optimization method based on Taylor expansion according to an embodiment of the present invention with reference to the accompanying drawings.

[0068] Specifically, Figure 1 A schematic flow chart of a method for topology optimization of a fiber-reinforced composite material structure based on Taylor expansion provided in an embodiment of the present invention.

[0069] like Figure 1 As shown in FIG, the fiber reinforced composite material structure topology optimization method based on Taylor expansion includes the following steps:

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

[0071] Specifically, the parameterization of the preset fiber-reinforced composite material (transversely isotropic material), the size of the finite element model to be solved, the size of the design domain (which may include non-design domains), the resolution of the Cartesian grid model, the geometric multigrid level, and the structural topology optimization parameters, including but not limited to the filter radius of the density and angle variables, the sampling point spacing of the angle variables, and other parameters are input into the host memory and copied to the GPU global memory.

[0072] 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.

[0073] Specifically, the density and angle design variables of each unit are declared and initialized in the GPU global memory, the boundary conditions and loads are input, the field variables for storing the unit stiffness matrix corresponding to each sampling point angle and the matrix of each order partial derivatives of the unit stiffness matrix with respect to each angle design variable are declared, the field variables for storing the block template matrix and preprocessing matrix corresponding to all nodes in each grid level are declared and initialized, and the variables required for the numerical processing process such as (design variable filtering, sensitivity filtering, regularization, etc.), and the field variables required for defining boundary conditions and loads, etc.

[0074] 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 and the partial derivative matrix of each unit stiffness matrix with respect to the 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 and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angles are stored in the preset GPU global memory.

[0075] Furthermore, in one embodiment of the present invention, the angle variable sampling points are determined according to a preset sampling interval, so as to calculate the element stiffness matrix corresponding to each group of sampling angles and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle according to the preset property parameters of the fiber-reinforced composite material and the angle variable sampling points, including:

[0076] Determine the angle variable sampling points within the value range of each angle design variable by presetting the sampling interval;

[0077] Constructing the constitutive relation matrix of orthotropic materials based on the preset property parameters of fiber-reinforced composite materials;

[0078] According to the constitutive relationship matrix of the orthotropic material and the angle variable sampling points, the element stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle are calculated.

[0079] Specifically, if Figure 2 As shown in the figure, a certain number of angle sampling points are generated within the range of each angle based on the preset angle sampling spacing. The preset fiber-reinforced composite material property parameters are substituted into the constitutive relationship matrix of the orthotropic material to calculate the material constitutive relationship matrix corresponding to each angle sampling point. The material parameters involved in the calculation process include the longitudinal elastic modulus, transverse elastic modulus, Poisson's ratio, and in-plane shear modulus of the fiber-reinforced composite material. The element stiffness matrix corresponding to each angle sampling point and the partial derivative matrix of the element stiffness matrix with respect to each variable are then calculated and stored in the GPU global memory.

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

[0081]

[0082] Where 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.

[0083] The partial derivative matrix of the element stiffness matrix corresponding to each set of sampling point angles with respect to the sampling angle is as follows:

[0084]

[0085] Where 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.

[0086] The partial derivative matrix of the material's elastic coefficient matrix with respect to the sampling angle can be calculated as follows:

[0087]

[0088] in,

[0089] 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.

[0090] 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.

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

[0092] 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:

[0093]

[0094] 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;

[0095] 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

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

[0097] 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

[0098] In step S105, based on the inter-unit parallel method, the unit stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angle are processed according to the Taylor expansion formula of the preset order to solve the approximate value of the unit stiffness matrix.

[0099] Among them, since the SIMP method is used to solve the approximate value of the element stiffness matrix, its expression is:

[0100]

[0101]

[0102] in, 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 e and θ e The mapped sampling angle.

[0103] In step S106, the template stiffness matrix of each node of each grid level is calculated according to the approximate value of the 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.

[0104] Specifically, the calculation process of the template stiffness matrix of each level required by the assembly-free geometric multigrid 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 adjacent units of the central node. When the stiffness matrix of a certain adjacent unit is used to calculate the contribution of its node template matrix, the filtered angle variable in the unit It will be rounded to the sampling point angle value closest to its value, and the sampling angle (α s ,θ s ) The corresponding unit stiffness matrix and the matrix of partial derivatives of the unit stiffness matrix with respect to the sampling angle are obtained. The order of Taylor expansion is selected according to the required accuracy, and Taylor expansion is performed at this point to obtain the approximate value of the accurate unit stiffness matrix.

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

[0106]

[0107] in, 2h A i x is the template stiffness matrix of each node on a coarser grid level, i s is the summation 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.

[0108] In step S107, 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.

[0109] Specifically, if Figure 4 As shown in Figure 2, the smoothed residual after preprocessing of the densest grid level is calculated based on the template stiffness matrix of each node at 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 at the densest grid level, including:

[0110] Based on the inter-element parallelism 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 each node at each mesh level, the preset initial displacement solution and the preset load of each node;

[0111] 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.

[0112] 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;

[0113] 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;

[0114] 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.

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

[0116] In step S108 , 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.

[0117] 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.

[0118] In step S109, 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.

[0119] 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-S109 are executed.

[0120] like Figure 5 As 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, and the 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%.

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

[0122] find: ρ, α, θ

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

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

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

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

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

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

[0129] 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 θ e They 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:

[0130]

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

[0132] 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.

[0133] 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 single stiffness matrix corresponding to each group of sampling angles is calculated as follows:

[0134]

[0135] Among them, B = [B0 B1 B2 B3 B4 B5 B6 B7] is a shape function matrix composed of 8 sub-matrices.

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

[0137]

[0138] Where 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. is the material constitutive relationship matrix corresponding to the sampling angle, which is calculated as follows:

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

[0140]

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

[0142]

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

[0144] The partial derivative matrix of the element stiffness matrix corresponding to each group of sampling point angles with respect to the sampling angle is calculated as follows:

[0145]

[0146] Among them, the partial derivative matrix of the material's elastic coefficient matrix with respect to the sampling angle can be calculated as follows:

[0147]

[0148]

[0149] in,

[0150]

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

[0152]

[0153] 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.

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

[0155] (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

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

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

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

[0159] Step 5: Calculate the template matrices of each level required by the assembly-free geometric multigrid preprocessing method. The calculation process is in parallel 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 adjacent units of the central node. When the stiffness matrix of a neighboring unit is used to calculate the contribution of its node template matrix, the filtered angle variable in the unit It will be rounded to the sampling point angle value closest to its value, and the sampling angle (α) calculated in step 3 can be directly accessed from the GPU global memory according to the index. s ,θ s) The corresponding element stiffness matrix and the matrix of partial derivatives of the element stiffness matrix with respect to the angle variable are calculated. The order of Taylor expansion is selected according to the required accuracy. Taylor expansion is performed at this point to obtain the approximate value of the accurate element stiffness matrix. The calculation process is as follows:

[0160]

[0161]

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

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

[0164] (2) Smoothing the densest level;

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

[0166] (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;

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

[0168] (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;

[0169] (7) Smoothing the densest level.

[0170] 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.

[0171] 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 .

[0172] 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;

[0173] Step nine, 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.

[0174] Step 10, such as Figure 6 and 7 As shown, determine whether the iteration has converged: determine whether the number of optimization iterations has reached 150. If it has reached 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.

[0175] 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.

[0176] Table 1 Comparison of results between the method according to the embodiment of the present invention and the traditional method

[0177] Method according to an embodiment of the present invention Traditional methods Number of grids 128×256×128 128×256×128 Runtime 1438.6s 4226.4s Initial design strain energy 15246.2 15246.2 Final strain energy (150 iterations) 3.34 3.34

[0178] 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 optimization algorithm proposed in the embodiments of the present invention significantly improves computational efficiency compared to traditional methods, significantly reducing the time cost of solving large-scale fiber-reinforced composite topology optimization problems. Furthermore, it achieves a reasonable balance between the storage space cost of the anisotropic material element stiffness matrix and the computational time cost, which is of great significance in practical applications.

[0179] In summary, the GPU-accelerated fiber-reinforced composite material topology optimization method based on Taylor expansion of the unit stiffness matrix proposed in an embodiment of the present invention greatly reduces the consumption of GPU global memory for the storage of the fiber-reinforced composite material unit stiffness matrix, and makes the storage consumption of the unit stiffness matrix only affected by the number of sampling points of the angle variable and no longer positively correlated with the model resolution; at the same time, the Taylor expansion is used to calculate the approximate unit stiffness matrix, compared with the method of calculating the accurate unit stiffness matrix by Gaussian integral, which greatly reduces the number of items involved in the cumulative calculation in the GPU single thread, saves memory access time and calculation time, and thus greatly improves the computational efficiency of the algorithm; the unit stiffness matrix approximation obtained by the Taylor expansion method has high accuracy, which is similar to the calculation of the accurate unit stiffness matrix, and the conjugate gradient iterative solver with geometric multigrid preconditioning converges to 10 -8 The number of iterations is almost the same, and the calculated node displacement values ​​can meet very high accuracy requirements. The objective function value is almost consistent with the method of solving the unit stiffness matrix with accurate angle.

[0180] Next, a fiber-reinforced composite material structure topology optimization device based on Taylor expansion according to an embodiment of the present invention will be described with reference to the accompanying drawings.

[0181] Figure 9 It is a block diagram of a fiber reinforced composite material structure topology optimization device based on Taylor expansion according to an embodiment of the present invention.

[0182] like Figure 9 As shown, the fiber-reinforced composite material structure topology optimization device 90 based on Taylor expansion 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 fourth calculation module 907, a solution and analysis module 908 and an iterative optimization module 909.

[0183] 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 and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle 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 and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle 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 process the element stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle using a Taylor expansion formula of a preset order in an inter-element parallel manner to obtain an approximate value of the element stiffness matrix. The third solution module 906 is used to calculate the template stiffness matrix of each node of each grid level based on the approximate value of the unit stiffness matrix, the template stiffness matrix solution format and the template stiffness matrix transfer relationship between each grid level, and store the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory. The fourth calculation module 907 is used to calculate the smoothed residual after preprocessing of the densest grid level based on 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 the preprocessed conjugate gradient method to obtain the current displacement of each node of the densest grid level. The solution and analysis module 908 is used to solve the objective function based on the current design variables and the current displacement of each node of the densest grid level, and perform sensitivity analysis on the objective function value and constraint conditions to obtain the sensitivity value. The iterative optimization module 909 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 the new optimized design is output.

[0184] It should be noted that the aforementioned explanation of the embodiment of the fiber-reinforced composite material structure topology optimization method based on Taylor expansion is also applicable to the fiber-reinforced composite material structure topology optimization device based on Taylor expansion in this embodiment, and will not be repeated here.

[0185] According to the GPU-accelerated fiber-reinforced composite material topology optimization device based on Taylor expansion of unit stiffness matrix proposed in an embodiment of the present invention, the consumption of GPU global memory for storage of fiber-reinforced composite material unit stiffness matrix is ​​greatly reduced, and the storage consumption of unit stiffness matrix is ​​only affected by the number of sampling points of angle variables and is no longer positively correlated with the model resolution; at the same time, compared with the method of calculating the accurate unit stiffness matrix by Gaussian integral, the method of calculating the approximate unit stiffness matrix by Taylor expansion greatly reduces the number of items involved in the cumulative calculation in a single thread of GPU, saves memory access time and calculation time, and thus greatly improves the calculation efficiency of the algorithm; the approximate value of the unit stiffness matrix fitted by the Taylor expansion method has high accuracy, which is similar to the calculation of the accurate unit stiffness matrix, and the conjugate gradient iterative solver with geometric multigrid preconditioning converges to 10 -8 The number of iterations is almost the same, and the calculated node displacement values ​​can meet very high accuracy requirements. The objective function value is almost consistent with the method of solving the unit stiffness matrix with accurate angle.

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

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

[0188] When the processor 1002 executes the program, the fiber-reinforced composite material structure topology optimization method based on Taylor expansion provided in the above embodiment is implemented.

[0189] Furthermore, the electronic device further includes:

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

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

[0192] 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.

[0193] 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.

[0194] 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.

[0195] 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.

[0196] An embodiment of the present invention further provides a computer-readable storage medium having a computer program stored thereon. When the program is executed by a processor, the method for optimizing the fiber-reinforced composite material structure topology based on Taylor expansion is implemented.

[0197] 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.

[0198] 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.

[0199] 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.

[0200] 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.

[0201] 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.

[0202] 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.

[0203] 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.

[0204] 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 is understood that the above embodiments are exemplary and are not to 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 method for topological optimization of fiber-reinforced composite materials based on Taylor expansion, 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 group of sampling angles and a partial derivative matrix of each unit stiffness matrix with respect to the 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 group of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the 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 unit stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angle are processed according to a Taylor expansion formula of a preset order to solve an approximate value of the unit stiffness matrix; Calculating the template stiffness matrix of each node of each grid level according to the approximate value of the unit stiffness matrix, the template stiffness matrix solution format and the template stiffness matrix transfer relationship between each grid level, and storing the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory; 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; 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 method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion 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 unit stiffness matrix corresponding to each group of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angle 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 and the partial derivative matrix of each element stiffness matrix with respect to the sampling angle are calculated according to the constitutive relationship matrix of the orthotropic material and the angle variable sampling points.

3. The method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion 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 method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion 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 method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion according to claim 1, characterized in that: The SIMP method is used to solve the approximate value of the element stiffness matrix, where the approximate value of the element stiffness matrix is: in, is the actual element stiffness matrix corresponding to the actual design variables, is the density penalty function 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 e and θ e The mapped sampling angle.

6. The method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion 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 summation 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.

7. The method for topology optimization of fiber-reinforced composite material structure based on Taylor expansion 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 mesh level is calculated according to the template stiffness matrix of each node of each mesh 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.

8. A fiber-reinforced composite material structure topology optimization device based on Taylor expansion, 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 group of sampling angles and a partial derivative matrix of each unit stiffness matrix with respect to the 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 group of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the 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 process the unit stiffness matrix corresponding to each set of sampling angles and the partial derivative matrix of each unit stiffness matrix with respect to the sampling angles according to a Taylor expansion formula of a preset order in an inter-unit parallel manner, and solve an approximate value of the unit stiffness matrix; a third solving module, configured to calculate the template stiffness matrix of each node of each grid level according to the approximate value of the unit stiffness matrix, the template stiffness matrix solving format, and the template stiffness matrix transfer relationship between each grid level, and store the template stiffness matrices of all levels except the densest grid level in the preset GPU global memory; a fourth calculation 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 Taylor expansion-based fiber-reinforced composite material structure topology optimization method 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 fiber-reinforced composite material structure topology optimization method based on Taylor expansion as described in any one of claims 1 to 7.

Images