CN115310154A - Continuous fiber 3D printing-oriented topological optimization and fiber path design method - Google Patents

Continuous fiber 3D printing-oriented topological optimization and fiber path design method Download PDF

Info

Publication number
CN115310154A
CN115310154A CN202210809160.5A CN202210809160A CN115310154A CN 115310154 A CN115310154 A CN 115310154A CN 202210809160 A CN202210809160 A CN 202210809160A CN 115310154 A CN115310154 A CN 115310154A
Authority
CN
China
Prior art keywords
fiber
unit
finite element
continuous
rho
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202210809160.5A
Other languages
Chinese (zh)
Inventor
金朋
张丰
昌敏
胡校斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN202210809160.5A priority Critical patent/CN115310154A/en
Publication of CN115310154A publication Critical patent/CN115310154A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B33—ADDITIVE MANUFACTURING TECHNOLOGY
    • B33Y—ADDITIVE MANUFACTURING, i.e. MANUFACTURING OF THREE-DIMENSIONAL [3-D] OBJECTS BY ADDITIVE DEPOSITION, ADDITIVE AGGLOMERATION OR ADDITIVE LAYERING, e.g. BY 3-D PRINTING, STEREOLITHOGRAPHY OR SELECTIVE LASER SINTERING
    • B33Y50/00—Data acquisition or data processing for additive manufacturing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Manufacturing & Machinery (AREA)
  • Materials Engineering (AREA)

Abstract

The invention relates to a topological optimization and fiber path design method for continuous fiber 3D printing, which is based on a continuous fiber 3D printing technology, considers failure constraint and takes structural rigidity as a target; integrating fiber layout, printing spacing and fiber orientation, and performing parallel optimization design on multiple scales; carrying out discrete design on the unit density by using a bidirectional progressive structure optimization method, and carrying out continuous variable optimization on the fiber content in the unit by using a variable density method without punishment on the basis; meanwhile, the main stress direction is used as the continuous fiber orientation, the failure of the composite material structure is considered, and the maximum failure point is approached by using the polymerized P-norm integration structure Tsai-hill rule; taking failure constraint as a supplement of a design target by introducing a Lagrange multiplier; finally, a set of method for planning the continuous fiber path is provided by means of an optimization scheme, and a result finished product is printed.

Description

Continuous fiber 3D printing-oriented topological optimization and fiber path design method
Technical Field
The invention relates to the technical field of continuous fiber 3D printing, in particular to a topological optimization and fiber path design method for continuous fiber 3D printing.
Background
The continuous fiber reinforced composite material has a series of excellent characteristics of high specific strength, large specific modulus, strong designability of material performance and the like, and has wide application in the fields of aviation, ships and automobiles. Due to the limitations of the conventional manufacturing process, the conventional structural optimization problem of the continuous fiber reinforced composite material is to optimize the stacking sequence of the laminated plates by considering the current loading condition under the condition of assuming uniform distribution of fibers and unchanged fiber orientation of each layer, so as to obtain the maximum mechanical property. However, the development of modern manufacturing processes has resulted in advanced composite manufacturing processes, such as continuous fiber-based 3D printing techniques, which are designed based on fiber angles in local regions, which can provide greater design freedom, and which can construct fiber orientation optimization problems for composite structures based on the manufacturing constraints of the designed parts, so as to more effectively distribute the load of the entire structure and maximize the rigidity of the structure.
The topological optimization is one of structural design, and the core idea of the topological optimization is to disperse a design domain into a plurality of units, represent the pseudo density of each unit by data between 0 and 1, and express the existence of the unit in the structure by forming a corresponding interpolation relation between the pseudo density and material properties such as volume, elastic modulus and the like. Although the properties of continuous fiber 3D printed prototypes have not been directly used, advanced manufacturing processes combined with structural design methods often have complementary effects, and isotropic materials create many high-performance, lightweight structures by combining topological optimization with 3D printing processes. However, the study on the topological optimization of the composite material is still in the development stage, and the topological optimization of the composite material needs to be combined with the consideration of the influence of the fiber orientation angle. The existing related research on the topological optimization of the composite material shows that the selection of the fiber angle has important influence on the topological form and the mechanical property, so that the combination of the topological optimization and the fiber orientation design is a promising design method.
In addition, the degree of freedom in the 3d printing process is also considered, the printing distance (which can be equal to the fiber volume fraction of the area) can be actually used as a design variable during printing, and the stress state is improved by theoretical analysis. The research on the printing distance fully utilizes the degree of freedom of the 3d printing composite material topology optimization, demonstrates the help of adding the printing distance as a variable to the optimization result, and develops a new idea for the composite material topology optimization.
The combined optimization study on fiber angle and topological configuration reveals the performance improvement effect of the optimized structure, dani
Figure BDA0003739811710000011
el Peeters et al solve the composite topology optimization problem by using a variable density method, indirectly express fiber angles by using lamination parameters and perform joint optimization, and obtain fiber path distribution through some post-processing, thereby obtaining a better result; haoqing Ding et al uses a normal distribution function to form a material interpolation scheme, thereby realizing continuous fiber angle design and establishing a parallel topology optimization model; xiaolei Yan et al established a structure/material topology and material orientation parallel optimization design method based on a bidirectional progressive topology optimization (BESO) method based on a uniform porous microstructure.
In addition, some researches consider the degree of freedom in the 3D printing process, for example, the printing distance (which can be equal to the fiber volume fraction of the area) is taken as a design variable during printing, and theoretical analysis is carried out to improve the stress state; the fiber is defined as a continuous curve distributed along a main stress motion track by A.V.Malakhov et al, and the influence of the change of the fiber orientation and the fiber distance on the stress concentration coefficient is researched by a finite element method on the basis, and the result shows that the stress concentration coefficient is obviously reduced by the optimization method. Then, the results of a die pressing experiment performed on the test pieces printed by the test pieces show that the strength and the rigidity after optimization are also obviously improved; zhanghao Hou et al design based on stress gradient distribution to make fiber content correspond to stress distribution and fiber direction correspond to maximum principal stress distribution, obtain trajectory lines with unequal fiber spacing, and obtain good results through testing; hang Li et al propose a full-scale optimization scheme, and simultaneously optimize the topology configuration, fiber orientation and fiber spacing of the structure, the scheme is established on a dual-material unit density optimization framework, and the fiber spacing is constrained by constraining the local fiber volume fraction, so that the fiber spacing can realize equidistant change in the whole world, and the full-scale optimization design is realized to a certain extent. The research on the printing distance fully utilizes the degree of freedom of the 3d printing composite material topology optimization, demonstrates the help of adding the printing distance as a variable to the optimization result, and develops a new idea for the composite material topology optimization.
The topology optimization in the prior art mainly has the following defects: (1) the optimization considered scale is incomplete; at present, the optimization design aiming at fiber angles is more, some optimization is carried out simultaneously by considering the fiber angles and topology optimization, the optimization is carried out by considering the fiber spacing, and an optimization framework for the full scale of a continuous fiber printing structure is lacked; (2) Manufacturing problems are not considered in the optimization process, most of the existing optimization methods only obtain some structures theoretically, actual printing and manufacturing are not realized, and manufacturing constraints in the printing process such as fiber continuity and printing feasibility are not considered; (3) the optimization result cannot be directly used for printing; and (4) the optimization does not take the composite material failure constraint into consideration.
Disclosure of Invention
The technical problem to be solved by the invention is to provide a topological optimization and fiber path design method for continuous fiber 3D printing, the method is based on a fiber mixing rate criterion of a composite material, synchronous optimization of fiber content, fiber orientation and topological configuration is carried out under the constraint of aggregation failure, the optimization scale is comprehensive, and a non-equidistant fiber path planning method is further used for generating continuous fiber tracks so as to achieve the result of actual printing and manufacturing.
The technical scheme adopted by the invention is that a topological optimization and fiber path design method for continuous fiber 3D printing, which comprises the following steps:
(1) Establishing a continuous fiber layout optimization model by adopting finite element software, and dispersing a design domain of the continuous fiber layout optimization model into a plurality of finite element units, wherein each finite element unit comprises four nodes;
(2) Initializing design variables of the continuous fiber layout optimization model by using the discrete design domain in the step (1), wherein the design variables comprise unit density rho B Unit fiber orientation angle theta and unit fiber content index rho S Outputting a finite element calculation file by the continuous fiber layout optimization model;
(3) Establishing a continuous fiber material interpolation model according to a composite material mixing ratio criterion, and outputting a material card of the continuous fiber material by the continuous fiber material interpolation model; defining optimization parameters in the matlab, reading a finite element calculation file output by the continuous fiber layout optimization model by the matlab, and obtaining the unit stress, the unit displacement and the unit main stress angle of each finite element unit in the continuous fiber material; each unit stress comprises four node stresses, and each unit displacement comprises four node displacements;
(4) Setting an optimization objective function of the continuous fiber composite material under the constraint of polymerization failure as follows:
Figure BDA0003739811710000031
Figure BDA0003739811710000032
wherein F represents the optimization objective function of the continuous fiber composite material under the constraint of polymerization failure, C represents the overall flexibility of the structure of the continuous fiber composite material, lambda is Lagrange multiplier, and FI max 、FI min Respectively representing the maximum and minimum values of the failure index value range, FI pn To representThe failure indexes are globally aggregated, and U, F and K respectively represent the integral displacement of the continuous fiber composite material structure, the integral load vector of the continuous fiber composite material structure and the integral rigidity matrix of the continuous fiber composite material structure; u. of e And k e Respectively representing the element displacement vector and the element stiffness matrix, p B Representing cell density, p S Represents a unit fiber content index, and theta represents a fiber orientation angle of the unit; v (x) and V * Respectively representing the total structural volume and the target volume; vf and vf * Respectively, the structural fiber content and the target fiber content;
(5) Calculating the element flexibility of each finite element unit according to the node stress and the node displacement of each finite element unit obtained in the step (3), and corresponding element density rho according to the obtained element flexibility of each finite element unit B And a unit fiber content index ρ S Respectively obtaining the cell density rho of each finite element unit by derivation B Sensitivity of and element fiber content index ρ of each finite element S The sensitivity of the method is that finite element analysis is carried out according to the element rigidity matrix and the node displacement of each finite element unit to obtain the stress accompanying variable value of each finite element unit, and the element density rho under the aggregation failure constraint is calculated according to the stress accompanying variable value of each element B Sensitivity and unit fiber content index ρ S The sensitivity of (2);
(6) Polymerizing the continuous fiber composite material optimization objective function under the polymerization failure constraint set in the step (4), wherein the specific process comprises the following steps: for cell density ρ B Carrying out a plurality of iterations on lambda in the optimized objective function F, wherein in each iteration, the unit density rho under the aggregation failure constraint obtained by calculation in the step (5) is subjected to B And updating the cell density rho in the continuous fiber composite optimization objective function F by adopting a BESO method B The index rho of the unit fiber content in each iteration process S Keeping unchanged, updating the orientation angle theta of the unit fiber along with the main stress angle of the unit in each iteration process, exporting a finite element calculation file in each iteration, and exporting the finite element calculation file according to the fact that the finite element calculation file is exportedAnalyzing a stress result by using the limit element calculation file, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches the maximum cycle value;
(7) And (3) polymerizing the continuous fiber composite material optimization objective function set in the step (4) under the polymerization failure constraint again on the basis of the step (6), wherein the specific process is as follows: index p for unit fiber content S Carrying out a plurality of iterations on lambda in the optimized objective function, wherein in each iteration, the unit fiber content index rho under the aggregation failure constraint obtained in the step (5) is subjected to S And updating the unit fiber content index rho in the optimization objective function of the continuous fiber composite material by adopting a SIMP method S Cell density ρ during each iteration B Keeping the orientation angle theta of the unit fiber unchanged, updating along with the main stress angle of the unit in each iteration process, analyzing a stress result according to a finite element calculation file exported in each iteration, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches a maximum cycle value;
(8) Outputting a finite element calculation file after the iteration is finished;
(9) And (5) generating continuous fiber tracks by a non-equidistant fiber path planning method according to the finite element calculation file output in the step (7).
Preferably, in step (8), the specific process of generating continuous fiber tracks by the non-equidistant fiber path planning method according to the finite element calculation file output in step (9) comprises the following steps:
first, the element density rho of each finite element in the finite element calculation file output in the step (8) B All are discrete variables, the value of each discrete variable is 0 or 1, and the unit density rho is taken B Finite element 1 to form a topological configuration;
secondly, distributing sample points equidistantly on the long scale of the topological structure formed in the first step, and dividing the topological structure formed in the first step into a plurality of areas containing different cross section branch numbers according to the cross section size of the sample points; the finite element calculation file output in the step (8) is stored inOf all finite element units of (a) has a unit fiber content index ρ S The fiber content v of the finite element elements in the topological structure formed in step one f (ii) a The unit fiber content index rho S The value of the continuous variable is any value between 0 and 1; the fiber content v of the finite element elements in the topological structure formed in step one f The expression of (a) is:
Figure BDA0003739811710000041
Figure BDA0003739811710000042
thirdly, combining adjacent areas with the same cross-section branch number to obtain a plurality of combined branch areas, and forming a new topological configuration by the combined branch areas again; according to the fiber content v contained in different branches in each combined branch region f Distributing the fiber number to obtain the expression of the fiber number in each merged branch region as follows:
Figure BDA0003739811710000051
wherein R is i Indicating the number of fibers in any one of the branched regions; sum (. Rho.) Si ∩ρ B = 1) represents the cell density ρ in the branch region B The sum of the unit fiber contents of 1;
Figure BDA0003739811710000052
denotes the sum of the fibre contents before branching, R total Represents the total number of fibers before branching;
fourthly, after the distribution of the fiber numbers of different branches in each combined branch area is finished, distributing the fiber points on the sample points, wherein the specific process is as follows: according to the obtained fiber number in each branch region, the unit fiber content index rho before branching S And the fiber content v of the finite element elements in the topological structure obtained in the step two f To distribute the fiber points on the sample points according to the distribution principle that the fiber points are not distributed in the branchThe content v of fibers in finite element elements in the new topological configuration obtained in the third step included between points of the same fiber f Equal; due to the fiber content v of each finite element in the topological configuration f Are unequal, so non-equidistant fiber points are obtained;
fifthly, aiming at one complete fiber and all fiber points on one fiber, filtering the orientation angle theta of the unit fiber finally generated in the step (7) according to the distance weight to obtain the smooth fiber direction; the expression for filtering the final generated unit fiber orientation angle θ is:
Figure BDA0003739811710000053
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003739811710000054
all the sums of distances within the filtering radius are indicated,
Figure BDA0003739811710000055
represents the product of distance and fiber angle;
Figure BDA0003739811710000056
θ i denotes the fiber angle after filtration, θ j Fiber angle within the filtration range, cel the filtration range,
Figure BDA0003739811710000057
representing the distance, i.e. the weight factor, r sen Represents the filtration radius; and finally, according to a Hermite interpolation method, taking the fiber point as a sample point coordinate, taking the fiber direction after fairing as a first derivative, and obtaining a distribution curve of the fiber under the constraint of the minimum printing spacing.
Preferably, in step (3), defining the optimization parameters in matlab includes: design variables, maximum circulation value, target volume, failure index constraints, cell density ρ B Filter radius and unit fibre content index p S The filtration radius of (a);
preferably, in the step (3), the specific process of establishing the continuous fiber material interpolation model according to the composite material mixing ratio criterion and outputting the material card by the continuous fiber material interpolation model comprises the following steps:
(3.1) passing the fiber content value of the continuous composite material in the composite material structure through the unit fiber content index rho S Interpolation is carried out, and v is respectively obtained according to different maximum printing pitches dmax and minimum printing pitches dmin under the determined printing layer thickness fmax And v fmin ，v fmax Represents the maximum value of the fiber content, v, of the continuous composite material within the topological configuration fmin Represents the minimum value of the fiber content of the continuous composite material within the topological configuration; fiber content v for finite element elements of continuous composite material within a topological configuration f Linear interpolation is carried out, and the expression is as follows:
Figure BDA0003739811710000061
in the formula
Figure BDA0003739811710000062
Represents the unit fiber content index of the ith unit;
(3.2) according to the modulus of the composite material structure and the modulus of the fiber material, calculating according to a composite material mixing ratio formula to obtain the rho s The expression of the interpolation model of the interpolated continuous fiber material is as follows:
E1＝v f E f +(1-v f )E m ,
Figure BDA0003739811710000063
Figure BDA0003739811710000064
v12＝v f v f +(1-v f )v m
wherein E is m 、G m 、v m Respectively the elastic modulus, shear modulus and Poisson's ratio of the matrix material constituting the composite material；E m 、G m 、v m Is the elastic modulus, shear modulus and Poisson's ratio, E, of the fibrous material 1 、E 2 、G 12 、v 12 The longitudinal elastic modulus, transverse elastic modulus, in-plane shear modulus and poisson's ratio of the composite material fiber after mixing are respectively expressed.
(3.3) introducing a unit fiber orientation angle theta, and if m = cos theta and n = sin theta, the off-axis rotation matrix of the composite material is as follows:
Figure BDA0003739811710000065
(3.4) introducing a cell density ρ B ，ρ B E {0,1}, indicates the presence or absence of material, will ρ B And (3) bringing the structure model into an off-axis rotation matrix of the composite material to obtain an interpolated continuous fiber material structure model, wherein the expression is as follows:
Figure BDA0003739811710000066
preferably, in step (5), the corresponding cell density ρ is determined according to the cell compliance of each cell B Deriving the cell density rho B The sensitivity of (a) is determined by the adjoint method, and the expression is:
Figure BDA0003739811710000067
wherein the content of the first and second substances,
Figure BDA0003739811710000068
a stiffness matrix representing the ith cell; according to the unit flexibility of each unit and the corresponding unit fiber content rho S Obtaining the unit fiber content index rho by derivation S The sensitivity of (a) is determined by the adjoint method, and the expression is:
Figure BDA0003739811710000069
wherein the content of the first and second substances,
Figure BDA00037398117100000610
Figure BDA00037398117100000611
Figure BDA00037398117100000612
B i 、T i 、D i respectively representing a rigidity matrix, a strain displacement matrix, an off-axis matrix and a constitutive matrix of the ith unit.
Preferably, in step (5), the corresponding stress accompanying variable value is calculated according to the cell stiffness matrix and the node displacement corresponding to each cell, and the stress accompanying variable value is calculated according to the stress accompanying variableValue ofObtaining the unit density rho under the constraint of polymerization failure B And a unit fiber content index ρ S The specific process of the sensitivity of (1) comprises the following steps:
(5.1) Global aggregation FI to obtain failure index PN To p rho Bi The expression of (a) is:
Figure BDA0003739811710000071
Figure BDA0003739811710000072
wherein the content of the first and second substances,
Figure BDA0003739811710000073
Figure BDA0003739811710000074
v is the failure index matrix, σ i X represents tensile strength in the fiber direction, Y represents tensile strength perpendicular to the fiber direction,
Figure BDA0003739811710000075
represents in-plane shear strength;
(5.2) deriving the stress by the failure index of the ith unit, wherein the expression is as follows:
Figure BDA0003739811710000076
(5.3) mixing i ＝D 0 B i u i Substituting the expression in step (4.2) to obtain:
Figure BDA0003739811710000077
wherein L is i Satisfy u i ＝L i U，u i Represents the unit displacement vector, B i Representing a shape function, D i Representing a constitutive matrix;
(5.4) obtaining an expression according to a finite element overall balance equation:
Figure BDA0003739811710000078
(5.5) substituting the expression in the step (5.3) into the expression in the step (5.4) to obtain:
Figure BDA0003739811710000079
(5.6) by introducing an accompanying variable μ into the expression in step (5.5) Bi To obtain the expression:
Figure BDA00037398117100000710
Figure BDA00037398117100000711
wherein K represents a structural overall stiffness matrix;
(5.7) simplifying the sensitivity of the P norm aggregation failure constraint according to the expression of the step (5.6) into the following steps:
Figure BDA00037398117100000712
Figure BDA00037398117100000713
due to the following:
Figure BDA00037398117100000714
finally solving the P norm aggregate failure index constraint pair design variable rho Bi The expression for sensitivity of (a) is:
Figure BDA00037398117100000715
the P norm aggregation failure index pair design variable rho can be obtained by the same method Si The expression for sensitivity of (a) is:
Figure BDA00037398117100000716
wherein
Figure BDA00037398117100000717
A companion node value vector for the ith element versus the unit fiber content index;
(5.8) combining the sensitivity of the aggregation failure constraint to the design variables obtained in the step (5.7) with the sensitivity of the P-norm aggregation failure index to the design variables to finally obtain the final sensitivity of the target function to the unit density, wherein the final sensitivity of the target function to the unit density is
Figure BDA00037398117100000718
Figure BDA00037398117100000719
The final sensitivity of the objective function to the index of the unit fiber content is
Figure BDA00037398117100000720
Figure BDA00037398117100000721
The invention has the beneficial effects that: by adopting the continuous fiber 3D printing-oriented topological optimization process, the topological optimization in the optimization process is based on the fiber mixing rate criterion of the composite material, synchronous optimization of fiber content, fiber orientation and topological configuration is performed under the constraint of polymerization failure, the optimization scale is comprehensively considered, the accuracy is high, and the optimization efficiency is high. On the basis of the topological optimization structure, continuous fiber tracks are further generated by a non-equidistant fiber path planning method, and the method has manufacturability.
Drawings
FIG. 1 is a flow chart of a topology optimization and fiber path design method for continuous fiber 3D printing according to the present invention;
FIG. 2 is an abstract view of a model for optimizing the placement of continuous fibers according to the present invention;
FIG. 3 is an enlarged view of the fiber phase angle at A in FIG. 2;
FIG. 4 is an enlarged view of the non-equidistant fibers at B in FIG. 3;
FIG. 5 is a schematic diagram of a printing path planning process of the cantilever beam according to the present invention;
1. a substrate; 2. fibers; 3. a blank region; 4. the fiber phase angle.
Detailed Description
The invention is further described below with reference to the accompanying drawings in combination with specific embodiments so that those skilled in the art can practice the invention with reference to the description, and the scope of the invention is not limited to the specific embodiments.
The embodiment of the invention provides a topology optimization and fiber path design method for continuous fiber 3D printing, which comprises the following steps of:
(1) Establishing a continuous fiber layout optimization model by adopting finite element software, and dispersing a design domain of the continuous fiber layout optimization model into a plurality of finite element units, wherein each finite element unit comprises four nodes;
(2) Initializing design variables of the continuous fiber layout optimization model by using the discrete design domain in the step (1), wherein the design variables comprise unit density rho B Unit fiber orientation angle theta and unit fiber content index rho S Outputting a finite element calculation file by the continuous fiber layout optimization model;
(3) Establishing a continuous fiber material interpolation model according to a composite material mixing ratio criterion, and outputting a material card of the continuous fiber material by the continuous fiber material interpolation model; defining optimization parameters in the matlab, reading a finite element calculation file output by the continuous fiber layout optimization model by the matlab, and obtaining the unit stress, the unit displacement and the unit main stress angle of each finite element unit in the continuous fiber material; each unit stress comprises four node stresses, and each unit displacement comprises four node displacements;
(4) Setting an optimization objective function of the continuous fiber composite material under the constraint of polymerization failure as follows:
Figure BDA0003739811710000081
Figure BDA0003739811710000091
wherein F represents the continuous fiber composite material optimization objective function under the constraint of polymerization failure, C represents the overall flexibility of the continuous fiber composite material structure, and lambda is Lagrange multiplier, FI max 、FI min Respectively representing the maximum and minimum values of the failure index value range, FI pn Representing the global aggregation of the failure indexes, wherein U, F and K respectively represent the integral displacement of the continuous fiber composite material structure, the integral load vector of the continuous fiber composite material structure and the integral rigidity matrix of the continuous fiber composite material structure; u. of e And k e Respectively representing the element displacement vector and the element stiffness matrix, p B Representing cell density, p S Represents a unit fiber content index, θ represents a fiber orientation angle of the unit; v (x) and V * Respectively representing the total structural volume and the target volume; vf and vf * Respectively, the structural fiber content and the target fiber content;
(5) Calculating the element flexibility of each finite element unit according to the node stress and the node displacement of each finite element unit obtained in the step (3), and corresponding element density rho according to the obtained element flexibility of each finite element unit B And a unit fiber content index ρ S Respectively obtaining the cell density rho of each finite element unit by derivation B Sensitivity of and element fiber content index ρ of each finite element S The finite element analysis is carried out according to the element rigidity matrix and the node displacement of each finite element unit to obtain the stress accompanying variable value of each finite element unit, and the element density rho under the aggregation failure constraint is calculated according to the stress accompanying variable value of each element B Sensitivity and unit fiber content index ρ S The sensitivity of (2);
(6) Polymerizing the continuous fiber composite material optimization objective function under the polymerization failure constraint set in the step (4), wherein the specific process comprises the following steps: for cell density ρ B Carrying out a plurality of iterations on lambda in the optimized objective function F, wherein in each iteration, the unit density rho under the aggregation failure constraint obtained by calculation in the step (5) is subjected to B And updating the cell density rho in the continuous fiber composite optimization objective function F by adopting a BESO method B The index rho of the unit fiber content in each iteration process S Keeping the orientation angle theta of the unit fiber unchanged, updating the orientation angle theta of the unit fiber along with the main stress angle of the unit fiber in each iteration process, exporting a finite element calculation file in each iteration process, analyzing a stress result according to the exported finite element calculation file, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches the maximum cycle value;
(7) And (5) on the basis of the step (6), polymerizing the continuous fiber composite material optimization objective function under the polymerization failure constraint set in the step (4) again, wherein the specific process comprises the following steps: for unit fiber content index ρ S Performing a plurality of iterations on the lambda in the optimized objective function, wherein the lambda in the step (6) has no relation with the lambda in the iteration, and the unit fiber content index rho under the polymerization failure constraint obtained in the step (5) is subjected to each iteration S And updating the unit fiber content index rho in the optimization objective function of the continuous fiber composite material by adopting a SIMP method S Cell density ρ in each iteration B Keeping the orientation angle theta of the unit fiber unchanged, updating along with the main stress angle of the unit in each iteration process, analyzing a stress result according to a finite element calculation file exported in each iteration, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches a maximum cycle value;
(8) Outputting a finite element calculation file after the iteration is finished;
(9) And (5) generating continuous fiber tracks by a non-equidistant fiber path planning method according to the finite element calculation file output in the step (7).
By adopting the topology optimization and fiber path design method facing continuous fiber 3D printing, the topology optimization in the method is based on the fiber mixing rate criterion of the composite material, synchronous optimization of fiber content, fiber orientation and topology configuration is performed under the constraint of polymerization failure, the optimization scale is considered comprehensively, the accuracy is high, on the basis of the topology optimization structure, continuous fiber tracks are generated by a non-equidistant fiber path planning method, and the manufacturability is achieved.
In the steps (1) to (8), a topology optimization framework for continuous fiber 3D printing is described, the topology optimization framework is proposed on the basis of expanding to an orthotropic material after combining a bidirectional gradual optimization method (BESO) and a variable density method (SIMP) in isotropic topology optimization, and an abstract schematic diagram of the design framework is shown in the figure; the design domain is dispersed by using a four-node finite element grid, a plurality of design variables are simultaneously optimized by a parallel optimization frame under the stress constraint, the discrete unit density obtained by a BESO method is used as material distribution, the continuous fiber content distribution is obtained by a non-punished SIMP method, and the main stress direction is used as the fiber direction; then, a material interpolation model with multiple design variables is introduced to express the relation between the objective function and the variables, aggregation constraint is added through a Lagrange multiplier on the basis that the objective function is the maximum in rigidity to prevent failure, and an updating scheme for determining the variables through sensitivity derivation of the design variables is provided. Searching an optimal structure in a given design domain according to the optimization model; in order to further approach the manufacturing requirement, the printing path of the obtained optimized structure is planned and tested.
In the process of topological optimization for continuous fiber 3D printing, fiber density in local areas is not uniform, so that structural fracture and material failure caused by high stress values or stress concentration in local fiber sparse areas are avoided, and test results are seriously influenced. Therefore, the method has important significance in researching the problem of the composite material failure constraint. For the stress constraint, the structure is required to beAll the positions are satisfied, the change of local parameters is sensitive, a large number of constraint conditions need to be satisfied by adopting a discretization numerical optimization method, and the difficulty and the calculation cost of optimization solution are obviously brought by excessive constraint conditions. To reduce the constraint conditions, these numerous local constraints can be written into a global overall constraint form, whose expression is:
Figure BDA0003739811710000101
Figure BDA0003739811710000102
wherein p represents the stress concentration; FI i According to the Tsai-Hill criterion
Figure BDA0003739811710000103
To obtain
Figure BDA0003739811710000111
It should be noted that a larger p can increase the weight of the maximum stress, reduce the chance of local stress overrun, but can also cause iterative oscillation more easily; sigma i Is the stress vector of the cell. Von-mises stress is generally used as a criterion in the isotropic field, and the Tsai-hill criterion is widely applied to strength calculation in anisotropic materials, so that the Tsai-hill criterion is used as a criterion in the composite material failure constraint.
In order to make the obtained topological optimization result have manufacturability, the method of the invention researches the path planning problem, and the current optimization usually involves the planning and processing of the path in the later period rarely; the invention introduces a fiber filling mode suitable for the topological optimization method to solve the fiber discontinuity after optimization and ensure the manufacturability.
In step (9), the specific process of generating continuous fiber tracks by a non-equidistant fiber path planning method according to the finite element calculation file output in step (8) comprises the following steps:
first, the element density rho of each finite element in the finite element calculation file output in the step (8) B All are discrete variables, the value of each discrete variable is 0 or 1, and the unit density rho is taken B A finite element of 1 to form a topological configuration;
secondly, distributing sample points equidistantly on the long scale of the topological structure formed in the first step, and dividing the topological structure formed in the first step into a plurality of areas containing different cross section branch numbers according to the cross section size of the sample points; the element fiber content index rho of all finite element units in the finite element calculation file output in the step (8) S The fiber content v of the finite element elements in the topological structure formed in step one f (ii) a The unit fiber content index rho S The value of the continuous variable is any value between 0 and 1; the fiber content v of the finite element elements in the topological structure formed in step one f The expression of (c) is:
Figure BDA0003739811710000112
Figure BDA0003739811710000113
thirdly, combining adjacent areas with the same number of cross-section branches to obtain a plurality of combined branch areas, and forming a new topological configuration by the combined branch areas again, wherein the number of cross-section branches included in the combined branch areas is equal; according to the fiber content v contained in different branches in each combined branch region f And distributing the fiber number to obtain the expression of the fiber number in each combined branch region as follows:
Figure BDA0003739811710000114
wherein R is i Indicating the number of fibers in any one of the branched regions; sum (. Rho.) Si ∩ρ B = 1) represents the cell density ρ in the branch region B The sum of the unit fiber contents of 1;
Figure BDA0003739811710000115
denotes the sum of the fibre contents before branching, R total Indicates to branchTotal number of fibers before;
fourthly, after the distribution of the fiber numbers of different branches in each combined branch area is finished, distributing the fiber points on the sample points, wherein the specific process is as follows: according to the obtained fiber number in each branch region, the unit fiber content index rho before branching S And the fiber content v of the finite element unit in the topological structure obtained in the step two f The fiber points on the sample points are distributed according to the fiber content v in the finite element unit in the new topological structure obtained in the step three among different fiber points in the branch f Equal; due to the fiber content v of each finite element in the topological configuration f Are unequal, so non-equidistant fiber points are obtained;
fifthly, aiming at one complete fiber and all fiber points on one fiber, filtering the orientation angle theta of the unit fiber finally generated in the step (7) according to the distance weight to obtain the smooth fiber direction; the expression for filtering the final generated unit fiber orientation angle θ is:
Figure BDA0003739811710000121
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003739811710000122
all the sums of distances within the filtering radius are indicated,
Figure BDA0003739811710000123
represents the product of distance and fiber angle;
Figure BDA0003739811710000124
θ i denotes the fiber angle after filtration, θ j Fiber angle within the filtration range, cel the filtration range,
Figure BDA0003739811710000125
representing the distance, i.e. the weight factor, r sen Represents the filtration radius; finally, according to a Hermite interpolation method, the fiber point is used as a sample point coordinateAnd taking the fiber direction after fairing as a first derivative, and obtaining a distribution curve of the fibers under the constraint of the minimum printing space.
In step (3), defining the optimization parameters in matlab includes: design variables, maximum circulation value, target volume, failure index constraints, cell density ρ B Filter radius and unit fiber content index ρ S The filtration radius of (a);
in the step (3), a continuous fiber material interpolation model is established according to a composite material mixing ratio criterion, and a specific process of outputting a material card by the continuous fiber material interpolation model comprises the following steps:
(3.1) passing the fiber content value of the continuous composite material in the composite material structure through the unit fiber content index rho S Interpolation is carried out, and v is respectively obtained according to different maximum printing pitches dmax and minimum printing pitches dmin under the determined printing layer thickness fmax And v fmin ，v fmax Represents the maximum value of the fiber content, v, of the continuous composite material within the topological configuration fmin Represents a minimum value of fiber content of the continuous composite material within the topological configuration; fiber content v for finite element elements of continuous composite material within a topological configuration f Linear interpolation is performed, and the expression is:
Figure BDA0003739811710000126
in the formula
Figure BDA0003739811710000127
Represents the unit fiber content index of the ith unit;
(3.2) according to the modulus of the composite material structure and the modulus of the fiber material, calculating according to a composite material mixing ratio formula to obtain the rho s The interpolation model of the interpolated continuous fiber material has the expression:
E1＝v f E f +(1-v f )E m ,
Figure BDA0003739811710000128
Figure BDA0003739811710000131
v12＝v f v f +(1-v f )v m
wherein, E m 、G m 、v m Respectively the elastic modulus, the shear modulus and the Poisson ratio of a base material forming the composite material; e m 、G m 、v m Is the elastic modulus, shear modulus and Poisson's ratio, E, of the fibrous material 1 、E 2 、G 12 、v 12 The longitudinal elastic modulus, transverse elastic modulus, in-plane shear modulus and poisson's ratio of the composite material fiber after mixing are respectively expressed.
(3.3) introducing a unit fiber orientation angle theta, and if m = cos theta and n = sin theta, the off-axis rotation matrix of the composite material is as follows:
Figure BDA0003739811710000132
(3.4) introducing a cell density ρ B ，ρ B E {0,1}, the presence or absence of material, will ρ B And (3) bringing the structure model into an off-axis rotation matrix of the composite material to obtain an interpolated continuous fiber material structure model, wherein the expression is as follows:
Figure BDA0003739811710000133
when ρ B When the stress strain matrix is 0, the stress strain matrix is empty, and the empty cell does not contribute to rigidity at the empty position; when rho B At 1, the stiffness matrix of the material depends on the values of other design variables. Rho s The fiber fraction is interpolated without punishment, and the rigidity, rho, of the material is influenced s Is actually penalized by p B Realization, for material free region, p s Does not cause the rigidity of the material to change, and because the stress-strain matrix is empty, ρ is s To p rho B Cell sensitivity of =0 approaches 0, and penalty is realizedAnd (6) carrying out the process. Theta affects the coordinate transformation matrix T and thus the cell stiffness matrix of the composite. The design method with completely free angles can cause the problems of too many variables, slow calculation speed and incapability of obtaining a calculation result, so that the topological optimization result is considered to be generally represented as a truss structure, and the stress condition of the branch rod piece is reasonable in the state of main stress. The fiber orientation theta is thus set in the calculation after each iteration as the principal stress direction of the last calculation.
In step (5), the corresponding cell density ρ is measured according to the cell compliance B Deriving the cell density rho B The sensitivity of (2) is determined by the adjoint method, and its expression is:
Figure BDA0003739811710000134
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003739811710000135
a stiffness matrix representing the ith cell; according to the unit flexibility to the corresponding unit fiber content rho S Obtaining the unit fiber content rho by derivation S The sensitivity of (a) is determined by the adjoint method, and the expression is:
Figure BDA0003739811710000136
wherein the content of the first and second substances,
Figure BDA0003739811710000137
Figure BDA0003739811710000138
Figure BDA0003739811710000139
B i 、T i 、D i respectively representing a rigidity matrix, a strain displacement matrix, an off-axis matrix and a constitutive matrix of the ith unit.
In the step (5), the corresponding stress accompanying variable value is obtained by calculation according to the unit stiffness matrix and the node displacement corresponding to each unit, and the stress accompanying variable value is obtained according to the stress accompanying variableValue ofTo obtain poly(s)For cell density ρ under resultant failure constraint B And unit fiber content ρ S The specific process of the sensitivity of (1) comprises the following steps:
(5.1) obtaining Global aggregation FI of failure index PN For rho Bi The expression of (a) is:
Figure BDA0003739811710000141
Figure BDA0003739811710000142
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003739811710000143
v is the failure index matrix, σ i X represents tensile strength in the fiber direction, Y represents tensile strength perpendicular to the fiber direction,
Figure BDA0003739811710000144
represents in-plane shear strength;
(5.2) deriving the stress by the failure index of the ith unit, wherein the expression is as follows:
Figure BDA0003739811710000145
(5.3) mixing i ＝D 0 B i u i Substituting the expression in step (4.2) to obtain:
Figure BDA0003739811710000146
wherein L is i Satisfy u i ＝L i U，u i Representing the unit displacement vector, B i Representing a shape function, D i Representing a constitutive matrix;
(5.4) obtaining an expression according to a finite element integral balance equation:
Figure BDA0003739811710000147
(5.5) substituting the expression in the step (5.3) into the expression in the step (5.4) to obtain:
Figure BDA0003739811710000148
(5.6) by introducing the accompanying variable μ into the expression in step (5.5) Bi And obtaining an expression:
Figure BDA0003739811710000149
Figure BDA00037398117100001410
wherein K represents a structural global stiffness matrix;
(5.7) finally simplifying the sensitivity of the P-norm polymerization stress according to the expression of the step (5.6) as follows:
Figure BDA00037398117100001411
Figure BDA00037398117100001412
because:
Figure BDA00037398117100001413
finally solving a P norm aggregation failure constraint pair design variable rho Bi The expression for sensitivity of (a) is:
Figure BDA00037398117100001414
the P norm aggregation failure constraint pair design variable rho can be solved in the same way Si The expression for sensitivity of (a) is:
Figure BDA00037398117100001415
wherein
Figure BDA00037398117100001416
A unit fiber content index for the ith element is accompanied by a node value vector;
(5.8) combining the sensitivity of the aggregation failure constraint obtained in (5.7) to the design variable with the sensitivity of the P-norm aggregation failure index to the design variable to finally obtain the final sensitivity of the target function to the unit density as
Figure BDA00037398117100001417
Figure BDA0003739811710000151
The final sensitivity of the objective function to the index of the unit fiber content is
Figure BDA0003739811710000152
Figure BDA0003739811710000153
The generation of continuous fiber trajectories by non-equidistant fiber path planning methods is the p in the topology optimization results S And rho B On the basis of dividing the result into a plurality of branch regions, creating a continuous optical fiber path within the range of the optimized result in a non-equidistant manner, the fiber path can be automatically generated according to the number of input fibers and sample points in each structure, and the fiber path arrangement of the cantilever beam is taken as an explanation below.
Taking the optimization result of the cantilever beam as an example, as shown in fig. 3, the path planning is designed. To optimize p in the result S And rho B On the basis, the continuous path is placed along the largest dimension of the structure (in this example, the x-direction); based on boundary conditions, starting from a final topological area, dividing the topological structure into a plurality of branch areas according to a sample point perpendicular line on the maximum dimension, and dividing a fiber point along the cross section of each sample point on the topological structure according to rho S Are not equidistantly distributed. Finally, the points are fairing along the direction of the largest dimension.
According to the method, firstly, sample points are distributed equidistantly on a long scale, and an optimization result is divided into a plurality of sections. And dividing the optimized macroscopic geometric characteristics into a plurality of areas according to the cross section of the sample point according to the unit density variable result fed back by the cross section. And converting the structural form change of the fiber into the change of the branch number, and branching and combining the fiber according to the change of the branch number.
When branching, the design variable ρ S Conversion to fiber content v f Total number of fibers before branchingWith the fibre content v of the individual branches f The result of (2) is to distribute the bifurcations, i.e. the fiber numbers, by using the fiber content as a weighting factor; the number of fibers per branch is:
Figure BDA0003739811710000154
Figure BDA0003739811710000155
wherein R is i Is the number of the branched fibers, R total The total number of fibers before branching. And when the fork is combined, the total number of the fibers after the fork is combined is obtained according to the sum of the branched fibers to be combined. By this method it is ensured that no cross-over between the fibres occurs.
After the total number of fibers to be distributed for each bifurcation has been calculated, the cross-section is determined on a cross-section by cross-section basis according to the sample points of the cross-sections within the bifurcation as described above
Figure BDA0003739811710000156
And rho S Calculating v f The fiber spacing within each bifurcation is not equally distributed, resulting in non-equidistant fiber sample points.
Then a filtering radius is given to the periphery of the distributed fiber points, and the fiber direction of the fiber points is taken as the average value of the fiber directions of all the units with the unit density of 1 in the filtering radius, namely:
Figure BDA0003739811710000157
for a fiber, the fiber is formed by connecting a series of points in the long dimension, in order to ensure the smoothness of the result, all the points on the fiber are the same as the sensitivity filtration, the fiber orientation is filtered by adopting the size of the measured distance, and the expression of filtering the optimized fiber angle theta is as follows:
Figure BDA0003739811710000158
Figure BDA0003739811710000159
the purpose of this filtration is to counteract some fiber pick-upThe effect on the more complex areas of variation, resulting in a smooth curve, as shown in fig. 4.
And finally, according to a Hermite interpolation method, taking the fiber point as a sample point coordinate, taking the fiber direction after fairing as a first derivative, and obtaining a distribution curve of the fiber under the constraint of the minimum printing spacing.

Claims (6)

1. A topology optimization and fiber path design method for continuous fiber 3D printing is characterized in that: the method comprises the following steps:
(1) Establishing a continuous fiber layout optimization model by adopting finite element software, and dispersing a design domain of the continuous fiber layout optimization model into a plurality of finite element units, wherein each finite element unit comprises four nodes;
(2) Initializing design variables of the continuous fiber layout optimization model by using the discrete design domain in the step (1), wherein the design variables comprise unit density rho B Unit fiber orientation angle theta and unit fiber content index rho S Outputting a finite element calculation file by the continuous fiber layout optimization model;
(3) Establishing a continuous fiber material interpolation model according to a composite material mixing ratio criterion, and outputting a material card of the continuous fiber material by the continuous fiber material interpolation model; defining optimization parameters in the matlab, reading a finite element calculation file output by the continuous fiber layout optimization model by the matlab, and obtaining the unit stress, the unit displacement and the unit main stress angle of each finite element unit in the continuous fiber material; each unit stress comprises four node stresses, and each unit displacement comprises four node displacements;
(4) Setting an optimization objective function of the continuous fiber composite material under the constraint of polymerization failure as follows:
Figure FDA0003739811700000011
Figure FDA0003739811700000012
wherein F represents the continuous fiber composite material optimization objective function under the constraint of polymerization failure, C represents the overall flexibility of the continuous fiber composite material structure, and lambda is Lagrange multiplier, FI max 、FI min Respectively representing the maximum and minimum values of the failure index value range, FI pn Representing the global aggregation of the failure indexes, wherein U, F and K respectively represent the integral displacement of the continuous fiber composite material structure, the integral load vector of the continuous fiber composite material structure and the integral rigidity matrix of the continuous fiber composite material structure; u. u e And k e Respectively representing the element displacement vector and the element stiffness matrix, p B Representing cell density, p S Represents a unit fiber content index, and theta represents a fiber orientation angle of the unit; v (x) and V * Respectively representing the total structural volume and the target volume; vf and vf * Respectively, the structural fiber content and the target fiber content;
(5) Calculating the element flexibility of each finite element unit according to the node stress and the node displacement of each finite element unit obtained in the step (3), and corresponding element density rho according to the obtained element flexibility of each finite element unit B And a unit fiber content index ρ S Respectively obtaining the cell density rho of each finite element unit by derivation B Sensitivity of and element fiber content index ρ of each finite element S The sensitivity of the method is that finite element analysis is carried out according to the element rigidity matrix and the node displacement of each finite element unit to obtain the stress accompanying variable value of each finite element unit, and the element density rho under the aggregation failure constraint is calculated according to the stress accompanying variable value of each element B Sensitivity and unit fibre content index p S The sensitivity of (2);
(6) Polymerizing the continuous fiber composite material optimization objective function under the polymerization failure constraint set in the step (4), wherein the specific process is as follows: for cell density ρ B Carrying out a plurality of iterations on lambda in the optimized objective function F, wherein in each iteration, the unit density rho under the aggregation failure constraint obtained by calculation in the step (5) is subjected to B The sensitivity of (a) is filtered out,and updating the unit density rho in the continuous fiber composite material optimization objective function F by adopting a BESO method B The index rho of the unit fiber content in each iteration process S Keeping the orientation angle theta of the unit fiber unchanged, updating the orientation angle theta of the unit fiber along with the main stress angle of the unit fiber in each iteration process, exporting a finite element calculation file in each iteration, analyzing a stress result according to the exported finite element calculation file, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches the maximum cycle value;
(7) And (3) polymerizing the continuous fiber composite material optimization objective function set in the step (4) under the polymerization failure constraint again on the basis of the step (6), wherein the specific process is as follows: for unit fiber content index ρ S Carrying out a plurality of iterations on lambda in the optimized objective function, wherein in each iteration, the unit fiber content index rho under the polymerization failure constraint obtained in the step (5) is subjected to S And updating the unit fiber content index rho in the optimization objective function of the continuous fiber composite material by adopting a SIMP (simple in-process map) method S Cell density ρ during each iteration B Keeping the orientation angle theta of the unit fiber unchanged, updating along with the main stress angle of the unit in each iteration process, analyzing a stress result according to a finite element calculation file exported in each iteration, and stopping iteration when the stress result meets the failure index constraint or the iteration reaches the maximum cycle value;
(8) Outputting a finite element calculation file after the iteration is finished;
(9) And (5) generating continuous fiber tracks by a non-equidistant fiber path planning method according to the finite element calculation file output in the step (7).
2. The continuous fiber 3D printing-oriented topology optimization and fiber path design method according to claim 1, wherein: in step (9), according to the finite element calculation file output in step (8), the specific process of generating continuous fiber tracks by the non-equidistant fiber path planning method comprises the following steps:
first, output in step (8)Element density ρ for each finite element in the finite element calculation file B All are discrete variables, the value of each discrete variable is 0 or 1, and the unit density rho is taken B A finite element of 1 to form a topological configuration;
secondly, distributing sample points equidistantly on the long scale of the topological structure formed in the first step, and dividing the topological structure formed in the first step into a plurality of areas containing different cross section branch numbers according to the cross section size of the sample points; the element fiber content index rho of all finite element units in the finite element calculation file output in the step (8) S The fiber content v of the finite element elements in the topological structure formed in step one f (ii) a The unit fiber content index rho S The value of the continuous variable is any value between 0 and 1; fiber content v of finite element elements in the topological structure formed in step one f The expression of (a) is:
Figure FDA0003739811700000031
Figure FDA0003739811700000032
thirdly, combining adjacent areas with the same number of branches on the cross section to obtain a plurality of combined branch areas, and forming a new topological configuration by the combined branch areas again; according to the fiber content v contained in different branches in each combined branch region f And distributing the fiber number to obtain the expression of the fiber number in each combined branch region as follows:
Figure FDA0003739811700000033
wherein R is i Represents the number of fibers in any branched region; sum (rho) Si ∩ρ B = 1) represents the cell density ρ in the branch region B The sum of the unit fiber contents of 1;
Figure FDA0003739811700000034
indicating the content of the fiber before branchingSum of amounts, R total Represents the total number of fibers before branching;
fourthly, after the fiber numbers of different branches in each combined branch area are distributed, distributing the fiber points on the sample points, wherein the specific process is as follows: according to the obtained fiber number in each branch region, the unit fiber content index rho before branching S And the fiber content v of the finite element unit in the topological structure obtained in the step two f The fiber points on the sample points are distributed according to the fiber content v in the finite element unit in the new topological structure obtained in the step three among different fiber points in the branch f Equal; due to the fiber content v of each finite element in the topological configuration f Are unequal, so non-equidistant fiber points are obtained;
fifthly, aiming at one complete fiber and all fiber points on one fiber, filtering the orientation angle theta of the unit fiber finally generated in the step (7) according to the distance weight to obtain the smooth fiber direction; the expression for filtering the final generated unit fiber orientation angle θ is:
Figure FDA0003739811700000035
wherein the content of the first and second substances,
Figure FDA0003739811700000036
all the sums of distances within the filtering radius are indicated,
Figure FDA0003739811700000037
represents the product of distance and fiber angle;
Figure FDA0003739811700000038
θ i denotes the fiber angle after filtration, θ j Fiber angle within the filtration range, cel the filtration range,
Figure FDA0003739811700000039
representing distance, i.e. weight factor，r sen Represents the filtration radius; and finally, according to a Hermite interpolation method, taking the fiber point as a sample point coordinate, taking the fiber direction after fairing as a first derivative, and obtaining a distribution curve of the fiber under the constraint of the minimum printing spacing.
3. The continuous fiber 3D printing-oriented topology optimization and fiber path design method according to claim 1, wherein: in step (3), defining the optimization parameters in matlab includes: design variables, maximum cycle value, target volume, failure index constraints, cell density ρ B Filter radius and unit fiber content index ρ S The filtration radius of (a);
4. the continuous fiber 3D printing-oriented topology optimization and fiber path design method according to claim 2, wherein: in the step (3), the specific process of establishing the continuous fiber material interpolation model according to the composite material mixing ratio criterion and outputting the material card by the continuous fiber material interpolation model comprises the following steps:
(3.1) passing the fiber content value of the continuous composite material in the composite material structure through the unit fiber content index rho S Interpolation is carried out, and v is respectively obtained according to different maximum printing pitches dmax and minimum printing pitches dmin under the determined printing layer thickness fmax And v fmin ，v fmax Represents the maximum value of the fiber content, v, of the continuous composite material within the topological configuration fmin Represents a minimum value of fiber content of the continuous composite material within the topological configuration; fiber content v for finite element elements of continuous composite material in topological configuration f Linear interpolation is performed, and the expression is:
Figure FDA0003739811700000041
in the formula
Figure FDA0003739811700000042
Represents the unit fiber content index of the ith unit;
(3.2) according to the composite materialThe modulus of the material structure and the modulus of the fiber material are calculated according to a composite material mixing ratio formula to obtain a product with a rho s The interpolation model of the interpolated continuous fiber material has the expression:
E1＝v f E f +(1-v f )E m ，
Figure FDA0003739811700000043
Figure FDA0003739811700000044
ν12＝v f ν f +(1-v f )ν m
wherein, E m 、G m 、v m Respectively the elastic modulus, the shear modulus and the Poisson ratio of a base material forming the composite material; e m 、G m 、v m Is the elastic modulus, shear modulus and Poisson's ratio, E, of the fibrous material 1 、E 2 、G 12 、v 12 The longitudinal elastic modulus, transverse elastic modulus, in-plane shear modulus and poisson's ratio of the composite material fiber after mixing are respectively expressed.
(3.3) introducing a unit fiber orientation angle theta, and if m = cos theta and n = sin theta, the off-axis rotation matrix of the composite material is as follows:
Figure FDA0003739811700000045
(3.4) introducing a cell density ρ B ，ρ B E {0,1}, the presence or absence of material, will ρ B And (3) bringing the structure model into an off-axis rotation matrix of the composite material to obtain an interpolated continuous fiber material structure model, wherein the expression is as follows:
Figure FDA0003739811700000046
5. the continuous fiber 3D printing-oriented topology optimization and fiber path design method according to claim 4, wherein: in step (5), the corresponding cell density rho is matched according to the cell compliance of each cell B Deriving the cell density rho B The sensitivity of (a) is determined by the adjoint method, and the expression is:
Figure FDA0003739811700000051
wherein the content of the first and second substances,
Figure FDA0003739811700000052
a stiffness matrix representing the ith cell; according to the unit flexibility of each unit, corresponding unit fiber content rho S Obtaining the unit fiber content index rho by derivation S The sensitivity of (a) is determined by the adjoint method, and the expression is:
Figure FDA0003739811700000053
wherein the content of the first and second substances,
Figure FDA0003739811700000054
Figure FDA0003739811700000055
Figure FDA0003739811700000056
B i 、T i 、D i respectively representing a rigidity matrix, a strain displacement matrix, an off-axis matrix and a constitutive matrix of the ith unit.
6. The continuous fiber 3D printing-oriented topology optimization and fiber path design method according to claim 5, wherein: in the step (5), corresponding stress accompanying variable values are obtained through calculation according to the unit stiffness matrix and the node displacement corresponding to each unit, and the unit density rho under the aggregation failure constraint is obtained according to the stress accompanying variable values B And sheetElementary fiber content index rho S The specific process of the sensitivity of (1) comprises the following steps:
(5.1) obtaining Global aggregation FI of failure index PN For rho Bi The expression of (a) is:
Figure FDA0003739811700000057
Figure FDA0003739811700000058
wherein the content of the first and second substances,
Figure FDA0003739811700000059
v is the failure index matrix, σ i X represents tensile strength in the fiber direction, Y represents tensile strength perpendicular to the fiber direction,
Figure FDA00037398117000000510
represents in-plane shear strength;
(5.2) deriving the stress by the failure index of the ith unit, wherein the expression is as follows:
Figure FDA00037398117000000511
(5.3) mixing i ＝D 0 B i u i Substituting the expression in step (4.2) to obtain:
Figure FDA00037398117000000512
wherein L is i Satisfy u i ＝L i U，u i Representing the unit displacement vector, B i Representing a shape function, D i Representing an constitutive matrix;
(5.4) obtaining an expression according to a finite element integral balance equation:
Figure FDA00037398117000000513
(5.5) substituting the expression in the step (5.3) into the expression in the step (5.4) to obtain:
Figure FDA00037398117000000514
(5.6) by introducing an accompanying variable μ into the expression in step (5.5) Bi To obtain the expression:
Figure FDA00037398117000000515
Figure FDA00037398117000000516
wherein K represents a structural overall stiffness matrix;
(5.7) simplifying the sensitivity of the P norm aggregation failure constraint according to the expression of the step (5.6) into the following steps:
Figure FDA00037398117000000517
Figure FDA00037398117000000518
due to the following:
Figure FDA00037398117000000519
finally solving the P norm aggregate failure index constraint pair design variable rho Bi The expression for sensitivity of (a) is:
Figure FDA0003739811700000061
the P norm aggregation failure index pair design variable rho can be obtained by the same method Si The expression for sensitivity of (a) is:
Figure FDA0003739811700000062
wherein, therein
Figure FDA0003739811700000063
A companion node value vector for the ith element versus the unit fiber content index;
(5.8) combining the sensitivity of the aggregated failure constraint obtained in (5.7) to the design variable with the sensitivity of the P-norm aggregated failure index to the design variable to obtain the final sensitivity of the objective function to the cell density as
Figure FDA0003739811700000064
Figure FDA0003739811700000065
The final sensitivity of the objective function to the index of the unit fiber content is
Figure FDA0003739811700000066
Figure FDA0003739811700000067
CN202210809160.5A 2022-07-11 2022-07-11 Continuous fiber 3D printing-oriented topological optimization and fiber path design method Pending CN115310154A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210809160.5A CN115310154A (en) 2022-07-11 2022-07-11 Continuous fiber 3D printing-oriented topological optimization and fiber path design method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210809160.5A CN115310154A (en) 2022-07-11 2022-07-11 Continuous fiber 3D printing-oriented topological optimization and fiber path design method

Publications (1)

Publication Number Publication Date
CN115310154A true CN115310154A (en) 2022-11-08

Family

ID=83856762

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210809160.5A Pending CN115310154A (en) 2022-07-11 2022-07-11 Continuous fiber 3D printing-oriented topological optimization and fiber path design method

Country Status (1)

Country Link
CN (1) CN115310154A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115577599A (en) * 2022-11-17 2023-01-06 东南大学 Motor topology optimization method based on component method

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115577599A (en) * 2022-11-17 2023-01-06 东南大学 Motor topology optimization method based on component method

Similar Documents

Publication Publication Date Title
Liu et al. Two-step homogenization of textile composites using mechanics of structure genome
CN107220461A (en) A kind of variation rigidity composite panel shell structure effectively optimizing method
CN107451307B (en) Method for multi-scale calculation of equivalent stiffness matrix of complex composite structure
Reddy On refined theories of composite laminates
CN106066913A (en) Complex composite material structure equivalent material performance multi-dimension computational methods
Esposito et al. Topology optimization-guided stiffening of composites realized through automated fiber placement
CN107451309B (en) Method for multi-scale calculation of equivalent thermal expansion coefficient of complex composite material structure
CN109241562B (en) Method for measuring elastic property of microstructure material based on multi-scale finite element method
CN107451308A (en) A kind of complex composite material structure effective thermal expansion coefficient multiscale simulation method
CN105808829A (en) CPU+GPU heterogeneous parallel computing based natural frequency characteristic analysis method for turbomachinery blade
CN111723457B (en) Level set method for optimization design of fiber curve laying variable-stiffness structure
Bayoumy et al. A continuum based three-dimensional modeling of wind turbine blades
Patni et al. On the accuracy of localised 3D stress fields in tow-steered laminated composite structures
CN115310154A (en) Continuous fiber 3D printing-oriented topological optimization and fiber path design method
CN116011301B (en) Finite element method for geometric state space such as B spline
Chakravarty On the modeling of composite beam cross-sections
Wang et al. Voronoi polygonal hybrid finite elements with boundary integrals for plane isotropic elastic problems
Littell et al. Effect of microscopic damage events on static and ballistic impact strength of triaxial braid composites
Esmaeeli et al. Simultaneous optimization of elastic constants of laminated composites using artificial bee colony algorithm
Li et al. Evaluation of piezoelectric and mechanical properties of the piezoelectric composites with local damages
Drumheller et al. Wave propagation in elastic laminates using a second order microstructure theory
Zhong et al. Isogeometric vibration and material optimization of rotating in-plane functionally graded thin-shell blades with variable thickness
CN103455713A (en) Novel method for designing mechanical modulus of planar weaving type composite material
CN108595724A (en) Composite material revolving meber design method
CN114428983A (en) Design method of bistable curvilinear fiber laminated plate

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination