Topological optimization design method and system for failure safety of fiber composite material under dynamic load

Abstract

The invention belongs to the technical field of structure optimization design, and particularly relates to a fiber composite failure safety topological optimization design method and system under a dynamic load, and the method comprises the steps: defining a design domain and related design parameters based on the stress working condition and geometric characteristics of a to-be-optimized structure; carrying out transient dynamic analysis and constructing an equivalent static load for topological optimization, constructing a multi-local failure working condition by taking material density, angle subinterval selection and a fiber angle as design variables, introducing an equivalent directional stiffness matrix continuously related to the fiber angle, and aggregating a multi-working-condition flexibility target by adopting a normalized KS function; failure safety topological optimization iteration and internal and external circulation convergence judgment are carried out in combination with a moving asymptote method, and collaborative optimization of structural topology and fiber orientation is achieved under the dynamic load condition. The problem that it is difficult to consider structural redundancy, failure safety and fiber orientation collaborative design in fiber reinforced composite topological optimization under the dynamic load is solved, and structural safety redundancy and local damage resistance are improved.

Description

Technical Field

[0001] This application belongs to the field of structural optimization design technology, and specifically relates to a method and system for topology optimization design for failure safety of fiber composites under dynamic loads. Background Technology

[0002] Fiber-reinforced composites are widely used in aerospace, automotive manufacturing, and civil engineering due to their high specific strength, high specific stiffness, excellent fatigue resistance, and corrosion resistance. The mechanical properties of these materials largely depend on the fiber arrangement and orientation design; structural performance can be improved through proper orientation optimization. Therefore, combining fiber orientation design with structural optimization methods has become an important research direction for enhancing the structural performance of composite materials.

[0003] Topology optimization, as an advanced design technique for efficiently utilizing materials and achieving lightweight engineering structures, can automatically find the optimal material distribution within a given design domain. However, classical topology optimization methods often focus on increasing overall stiffness or reducing compliance, which may lead to insufficient structural safety redundancy and reduced tolerance to local damage, defects, or manufacturing deviations.

[0004] For example, patent number CN119066908A proposes a topology optimization method for fiber-reinforced composite materials that considers residual stress. By introducing residual stress constraints into the topology optimization model, it achieves control over the residual stress level and avoids manufacturing defects. However, the optimized structure has relatively low safety redundancy, and its overall safety may still be challenged under dynamic loads or local failure conditions.

[0005] In engineering applications, structures typically require a certain degree of redundancy to ensure continued load-bearing capacity in the event of failure. Patent number CN117408077A proposes an optimization design method for offshore wind turbine jackets that considers failure safety. By optimizing the failure and non-failure regions with different objectives, the structure maintains good stability even after localized damage in marine environments, providing a framework for redundancy design.

[0006] On the other hand, traditional topology optimization is mostly designed for static load conditions, while dynamic conditions such as vibration, impact, and cyclic loads are widely present in actual engineering, causing the structural response to exhibit obvious time-varying characteristics. Topology optimization under dynamic loads not only has high computational complexity but also places higher demands on the evolution of material distribution. Patent No. CN106372347A proposes a dynamic response topology optimization method based on an improved two-way progressive method using an equivalent static load method. By introducing the improved two-way progressive method into the equivalent static load framework, the computational efficiency of the dynamic optimization process is improved.

[0007] In summary, while existing technologies have explored aspects such as topology optimization of composite materials, the influence of residual stress, failure-safe design, and topology optimization under dynamic loads, they still suffer from insufficient structural redundancy, difficulty in balancing dynamic stability, and insufficient coupling in composite material orientation design. Therefore, a new method is urgently needed that can coordinate structural redundancy, fiber orientation optimization, and topology optimization under dynamic loads to improve the overall performance and resistance to localized damage of fiber-reinforced composite structures. Summary of the Invention

[0008] To address the technical problems in existing technologies, such as insufficient redundancy in topology optimization of fiber-reinforced composite structures, difficulty in effectively handling dynamic load conditions, insufficient coupling between fiber orientation optimization and material distribution optimization, and significant performance degradation under local damage conditions, this application provides a failure-safe topology optimization design method and system for fiber composites under dynamic loads. The technical solution of this application can fully consider the local damage and failure conditions that may occur in the actual service process of fiber-reinforced composite structures under dynamic conditions such as vibration, impact, and periodic loads. By introducing a dynamic response characterization method of equivalent static load, a multi-local failure condition construction strategy, and an equivalent directional stiffness matrix that is continuously related to fiber angles, a unified solution framework for fiber orientation optimization, structural topology optimization, and failure-safe design is realized. Through the above optimization strategies, this application can not only improve the load-bearing capacity and stability of fiber-reinforced composite structures under dynamic loads, but also enhance their continuous load-bearing capacity and safety redundancy when local damage occurs, so that the structure will not suffer catastrophic failure after local damage. This enables lightweight design of fiber-reinforced composite structures while meeting performance requirements, and is particularly suitable for aerospace, vehicle equipment, marine equipment and other transportation equipment fields with extremely high requirements for structural safety and dynamic adaptability.

[0009] On the one hand, this application provides a failure-safe topology optimization design method for fiber composites under dynamic loads, the method comprising:

[0010] Step 1: Based on the stress conditions and geometric characteristics of the structure to be optimized, define the design parameters of the design domain for the fiber-reinforced composite structure;

[0011] Step 2: Based on the design domain constructed in Step 1, perform dynamic mechanical analysis on the fiber-reinforced composite structure and obtain the equivalent static load for topology optimization;

[0012] Step 3: Replace the dynamic load of the original structure with the equivalent static load. Take the model with the equivalent static load as the analysis object. Use the material density variable, angle sub-interval selection variable and fiber angle variable in the above model as design variables. Construct multiple local failure conditions and, based on the equivalent directional stiffness matrix and normalized KS objective function that are continuously related to the fiber angle, perform failure-safe topology optimization iteration on the design variables using the moving asymptote method until the inner loop converges.

[0013] Step 4: Determine whether the external loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load and the structural flexibility when constructing the equivalent static load is less than 0.1%. If it has converged, output the failure safety topology optimization result of the fiber-reinforced composite material at this time. Otherwise, perform dynamic analysis again on the updated structure and construct a new equivalent static load before returning to step 3 to continue optimization.

[0014] In a preferred implementation, step 1 further includes:

[0015] Step 1.1: Based on the load-bearing form, main deformation mode and typical stress characteristics of the fiber-reinforced composite structure to be optimized in practical applications, a three-point bending beam is selected as the structure to be optimized.

[0016] Step 1.2: Based on the geometric and force symmetry of the three-point bending beam, the right half of the three-point bending beam is selected as the design domain to construct a semi-three-point bending beam model;

[0017] Step 1.3 determines the design parameters of the design domain to establish the input conditions for topology optimization and fiber orientation optimization.

[0018] In a preferred implementation, step 2 further includes:

[0019] Step 2.1: Based on the geometric dimensions, boundary constraint arrangement, material properties and mesh generation of the design domain determined in Step 1, the design domain is discretized by finite element method, and the mass matrix M, damping matrix C and stiffness matrix K of the structure are generated to construct the transient dynamic control equations.

[0020] The transient dynamic governing equations for fiber-reinforced composite structures are:

[0021]

[0022] In the formula: The structural mass matrix; Here is the structural damping matrix; Here is the stiffness matrix of the structure; It is the displacement vector; It is the velocity vector; It is the acceleration vector; The dynamic load on the structure at time t;

[0023] Step 2.2: Based on the transient dynamic control equations, solve for the displacement vector of the structure at time m under dynamic load. ;

[0024] Step 2.3: Based on the displacement at time m Calculate the equivalent static load at that moment;

[0025] Step 2.4: Equivalent static load at the m-th time step obtained in Step 2.3 The equivalent static load data for n time steps are obtained, and the n equivalent static load vectors are weighted, summed, or averaged to obtain the final equivalent static load for topology optimization. ;

[0026] Final equivalent static load for:

[0027]

[0028] In the formula: This is the final equivalent static load; Let be the equivalent static load vector at time step m; K is the overall structural stiffness matrix. Let be the nodal displacement vector of the structure at time step m; n is the number of time steps used to construct the final equivalent static load.

[0029] In a preferred implementation, step 3 further includes:

[0030] Step 3.1: Establish a failure safety topology optimization model for fiber-reinforced composite materials with material density variable, angle sub-interval selection variable and fiber angle variable as design variables, and use the maximum value of the compliance under each damage condition as the optimization objective;

[0031] Step 3.2: Based on the finite element mesh, the candidate local damage blocks are laid out and translated according to the preset mesh size to generate multiple local failure regions. The finite element elements inside and outside each region are assigned the minimum weakening coefficient and 1 respectively to form the corresponding weakening coefficient field, so as to construct the local failure scenario under multiple damage conditions required in Step 3.1.

[0032] Step 3.3: Divide the original design interval of the fiber direction into several sub-intervals according to parameter n, calculate the directional elastic matrix corresponding to the central angle of each sub-interval, and introduce discrete-continuous weight variables and small deflection angles within the sub-intervals. Weight and rotate the directional elastic matrix to obtain the equivalent directional stiffness matrix that is continuously related to the fiber angle design variable. Substitute the equivalent directional stiffness matrix into the overall stiffness matrix expression to form the stiffness matrix used in Step 3.1.

[0033] Step 3.4: Based on the failure-safe topology optimization model established in Step 3.1 with the maximum flexibility of each local failure condition as the objective, and according to the structural flexibility under each local failure scenario obtained in Step 3.2, the maximum flexibility objective is transformed into a continuously differentiable objective function using the normalized KS aggregation function.

[0034] Step 3.5: Based on the KS objective function constructed in Step 3.4 and the explicit derivative relationship of compliance with design variables, calculate the sensitivity of the objective function to material density variables, angle sub-interval selection variables, and fiber angle variables, and use the moving asymptote method to iteratively update the three types of design variables based on the sensitivity until the inner loop converges.

[0035] In the preferred implementation, further, in step 3.1, the failure-safe topology optimization model of the fiber-reinforced composite material is constructed through the following constraints:

[0036]

[0037] In the formula: The set of variables that the solver needs to solve; The density of the material; Choose variables for the angle sub-intervals, n sub-intervals; For fiber angle variables; To minimize the objective function f in all m damage scenarios and improve the stiffness under the weakest condition; For the first The flexibility under each damage scenario; m is the number of damage scenarios; To find the worst softness in all damage scenarios; In the first Damage scenario, displacement Mechanical equilibrium must be satisfied; This is the equivalent static load; Here is the stiffness matrix for the l-th damage scenario; The displacement vector is obtained under the l-th damage scenario; Material usage constraints; The sum of the densities of all elements; N is the total number of elements; f is the maximum allowable material volume fraction; ; The variable chosen for the angle interval is also in [0,1]. To limit the fiber laying angle Scope; It is the fiber orientation angle; To limit the range of angle change.

[0038] In the preferred implementation, further, in step 3.2, the weakening coefficient of the e-th finite element is:

[0039]

[0040] In the formula: The design domain is given in step 1; For the first There are local failure regions; e is the e-th finite element element within the design domain; For this unit in a partial failure scenario The weakening coefficient below; This is the preset minimum weakening coefficient.

[0041] In the preferred implementation, further, in step 3.4, the maximum compliance objective is transformed into a continuously differentiable objective function using a normalized KS aggregation function:

[0042]

[0043] In the formula: The objective function value after aggregation by the KS function; is the natural logarithm function; m is the total number of local failure scenarios described in step 3.2; This represents the sequence number of the partial failure scenario. =1-m; For the first The structural compliance under equivalent static load in a local failure scenario; e is the base of the natural logarithm, representing the exponential function. ; For regularization parameters; This is the normalization coefficient.

[0044] In the preferred implementation, further, in step 3.5, the sensitivity of the objective function to the material density variable is:

[0045]

[0046] The sensitivity of the angle sub-interval selection variable is:

[0047]

[0048] The sensitivity of the fiber angle variable is:

[0049]

[0050] In the formula: Design variables for the material density of the e-th element for the objective function f. Sensitivity; is the objective function value after aggregation by the KS function; m is the total number of local failure conditions; This is the sequence number of the partial failure condition; For the first The structural flexibility under equivalent static load in a partial failure condition; For softness For material density variables The partial derivatives; Design variables for the volume distribution of the e-th unit. Right now ; This refers to the regularization parameter in the KS function; For the first The weighting factor corresponding to each local failure condition; e is the base of the natural logarithm, representing the exponential function. ; The sign for summing over all partial failure conditions; Choose variables for the objective function f for the i-th angle sub-interval in the e-th unit. Sensitivity; Design variables associated with the i-th angular sub-interval of the e-th unit; For softness Choose variables for angle sub-intervals The partial derivatives; The objective function f is the fiber angle variable in the e-th element. Sensitivity; Design variables for the fiber angle of the e-th finite element; In the first Compliance under partial failure conditions Design variables for unit fiber angle The partial derivatives of .

[0051] In a preferred implementation, further, in step 1.3, the stress condition of the semi-three-point bending beam under dynamic load is set as follows: a fixed constraint is applied to the lower right corner of the design domain, a time-varying dynamic load is applied to the upper left corner of the design domain, and a horizontal degree-of-freedom constraint is applied to the left end face of the design domain; the areas of a preset unit width on both sides of the design domain are set as non-failure regions. : x and y are spatial coordinate variables in the structural coordinate system where the design domain is located, where x is the coordinate in the horizontal direction along the length of the semi-three-point curved beam, and y is the coordinate in the vertical direction along the height of the semi-three-point curved beam; L is the length of the design domain along the beam length; and b is the preset unit width of the non-failure regions set on the left and right sides of the design domain along the beam length.

[0052] On the other hand, this application also provides a topology optimization design system for the failure safety of fiber composites under dynamic loads, the system comprising:

[0053] The design parameter definition module is used to define the design parameters of the design domain of the fiber-reinforced composite structure based on the stress conditions and geometric characteristics of the structure to be optimized.

[0054] The dynamic analysis and equivalent static load generation module is used to perform dynamic mechanical analysis on fiber-reinforced composite structures based on the design domain constructed by the design parameter definition module, and obtain the equivalent static load for topology optimization.

[0055] The failure-safe topology optimization iteration module is used to replace the dynamic load of the original structure with the equivalent static load obtained by the dynamic analysis and equivalent static load generation module. The model with the applied equivalent static load is used as the analysis object. The material density variable, angle sub-interval selection variable and fiber angle variable in the above model are used as design variables to construct multiple local failure conditions. Based on the equivalent directional stiffness matrix and normalized KS objective function that are continuously related to the fiber angle, the failure-safe topology optimization iteration is performed on the design variables in combination with the moving asymptote method until the inner loop converges.

[0056] The outer loop convergence judgment and result output module is used to determine whether the outer loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load obtained by the failure-safe topology optimization iteration module is less than 0.1%. When converged, the failure-safe topology optimization result of the fiber-reinforced composite material is output. Otherwise, the dynamic analysis is re-performed on the updated structure and a new equivalent static load is constructed before returning to the failure-safe topology optimization iteration module to continue optimization.

[0057] The beneficial effects of this application are:

[0058] First, the failure safety topology optimization design method for fiber composites under dynamic loads proposed in this application, compared with the prior art, transforms the complex dynamic response optimization problem into a series of static topology optimization problems by constructing an equivalent static load based on transient dynamic analysis. This reduces the computational load while preserving the time-varying characteristics of the dynamic load. Based on this, material density variables, angle sub-interval selection variables, and fiber angle variables are used as unified design variables. An equivalent directional stiffness matrix continuously related to the fiber angle is introduced, and a normalized KS objective function is constructed. Combined with modeling of multiple local failure conditions and iterative solution using the moving asymptote method, the collaborative optimization design of structural topology, fiber orientation, and failure safety constraints is achieved. This effectively improves the structure's tolerance to local damage, manufacturing defects, and orientation deviations, enhancing safety redundancy. Simultaneously, by achieving convergence and updating of design variables in the inner loop and determining the consistency between the equivalent static load and the structural response based on relative changes in flexibility in the outer loop, the obtained optimization results balance lightweighting and failure safety under dynamic load conditions, improving the overall load-bearing capacity and dynamic stability of fiber-reinforced composite structures.

[0059] Secondly, in the preferred implementation, in step 1.1, this application combines the actual load-bearing form and main deformation mode of the structure to be optimized, and takes a typical and standardized three-point bending beam as the object to be optimized, so that the test conditions are highly consistent with engineering applications, which facilitates the comparison and evaluation of the optimization effect of this method; in step 1.2, by utilizing the geometric and stress symmetry of the three-point bending beam, only the right half is selected as the design domain, which not only avoids the repeated modeling and solving of symmetrical regions, reducing the number of finite element meshes and computational costs, but also naturally satisfies the symmetry requirements of the overall structure by mirroring the optimization results of the design domain; furthermore, by pre-determining the design parameters such as the size, load, and boundary constraints of the design domain in step 1.3, the input conditions for topology optimization and fiber orientation optimization are uniformly constructed, which is conducive to improving the stability and feasibility of subsequent numerical analysis and optimization solutions.

[0060] Third, in the preferred implementation, this application further establishes the transient dynamic control equations of the fiber-reinforced composite material structure based on the mass matrix, damping matrix, and stiffness matrix through steps 2.1-2.4, and obtains the displacement response at each moment through numerical integration, so that the construction of the equivalent static load fully reflects the dynamic effects such as inertia and damping; on the other hand, the equivalent static load sequence at each moment is summed or averaged to obtain the final equivalent static load that can comprehensively represent the entire vibration process, thereby avoiding the randomness and conservatism caused by selecting only the peak moment load, and improving the representativeness and stability of the load condition to the real dynamic environment.

[0061] Fourth, in the preferred implementation, further through steps 3.1-3.5, on the one hand, a failure-safe topology optimization model is constructed with the material density variable, angle sub-interval selection variable, and fiber angle variable as unified design variables, aiming at maximizing the compliance of multiple local failure conditions. By generating multiple local weakening regions and corresponding weakening coefficients on the finite element mesh, a unified description of different failure locations and failure degrees is achieved, ensuring that the optimization results can maintain sufficient load-bearing capacity under various local failure conditions. On the other hand, using the fiber angle directional elastic matrix and the discrete-continuous parameterization method, a finite number of angle sub-intervals are smoothly mapped to an equivalent directional stiffness matrix continuously related to the fiber angle design variable. Then, the multi-condition compliance maximization problem is transformed into a continuously differentiable objective function through the normalized KS aggregation function, and the sensitivity expression for the three types of design variables is derived, which facilitates efficient gradient optimization by coupling with the moving asymptote method. Thus, while ensuring the stability and convergence of the solution, the collaborative optimization design of structural topology, fiber arrangement, and failure safety is achieved.

[0062] Fifth, the failure-safe topology optimization design system for fiber composites under dynamic loads of this application modularly integrates design parameter definition, dynamic analysis and equivalent static load generation, failure-safe topology optimization iteration, and external loop convergence judgment and result output into the same platform. This realizes an automated closed-loop process of dynamic load condition modeling, equivalent static load conversion, fiber orientation and failure-safe constraint collaborative optimization, and multi-level convergence control, which greatly reduces the difficulty of engineering implementation and improves the efficiency and reliability of optimization design. Attached Figure Description

[0063] Figure 1 This is a flowchart of the failure safety topology optimization design method for fiber-reinforced composite materials under dynamic load as described in an embodiment of the present invention;

[0064] Figure 2 This is a schematic diagram of the analysis domain outer loop and design domain inner loop of the failure safety topology optimization method for fiber-reinforced composite materials described in this embodiment of the invention;

[0065] Figure 3 This is a schematic diagram showing the geometric parameters, boundary conditions, and design domain division of the structure in an embodiment of the present invention;

[0066] Figure 4 This is a schematic diagram illustrating a localized structural failure scenario in an embodiment of the present invention.

[0067] Figure 5 This is a schematic diagram showing the distribution of localized structural failure areas in an embodiment of the present invention;

[0068] Figure 6 This is a schematic diagram of the topology optimization design results in an embodiment of the present invention, wherein: Figure 6In the diagram, a1-a3 represent the topology optimization design results of fiber-reinforced composite materials under dynamic load. Figure 6 b1-b3 in the figure represent the failure safety topology optimization design results of fiber-reinforced composite materials under dynamic load when the damage size L = 30.

[0069] Figure 7 In the embodiments of the present invention, respectively in Figure 6 A schematic diagram of the remaining structure obtained after applying a 30² post-damage patch to structures a2 and b2. Detailed Implementation

[0070] To enable those skilled in the art to better understand the technical solutions of this application, the following will provide a more detailed description of this application in conjunction with the accompanying drawings and embodiments.

[0071] The directional terms such as above, below, left, right, front, and back used in this application are based on the positional relationships shown in the attached drawings. Different attached drawings may result in different positional relationships, therefore they should not be interpreted as limitations on the scope of protection.

[0072] In this application, the terms "installation," "connection," "interlocking," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, an integral connection, a mechanical connection, an electrical connection, or a connection that allows communication between components. They can also refer to a direct connection or an indirect connection through an intermediate medium. They can refer to the internal connection of two components or the interaction between two components. For those skilled in the art, the specific meaning of the above terms in this application can be understood according to the specific circumstances.

[0073] This invention discloses a failure-safe topology optimization design method and system for fiber-reinforced composite materials under dynamic loads. This method simulates local damage by removing the stiffness of the composite material in a predefined damage patch region during topology optimization, effectively characterizing potential failure modes during actual service. Simultaneously, this invention employs a discrete-continuous fiber orientation parameterization method, improving the safety of fiber-reinforced composite structures under complex loading environments while enabling continuous control of fiber orientation, and achieving parallel optimization design of material distribution and fiber orientation. Furthermore, this invention transforms the original dynamic load condition into several equivalent static loads based on the equivalent static load method. This not only accurately reflects the structural response characteristics under dynamic conditions but also improves the computational efficiency of the failure-safe topology optimization solution process, making safe redundancy design of complex composite structures under dynamic load conditions possible.

[0074] As per the instruction manual Figure 1-2 The failure-safe topology optimization design method for fiber composites under dynamic loads of the present invention includes:

[0075] Step 1: Based on the stress conditions and geometric characteristics of the structure to be optimized, define the design parameters of the design domain of the fiber-reinforced composite structure.

[0076] The purpose of step 1 is to determine the design domain range, load conditions, boundary conditions, and non-failure regions based on the three-point bending beam, so that the design domain can form a unified, clear, and computable initial structural model that can be used for dynamic response analysis and failure-safe topology optimization, thereby ensuring the effectiveness and feasibility of material distribution design and fiber orientation design in the subsequent optimization process.

[0077] It should be noted that the design domain refers to the spatial region where material distribution optimization is permitted, typically represented as a computational region with defined geometry, boundary conditions, and load conditions. The load case refers to the various external loads a structure experiences during use or operation, including their mode of action, location, duration, and direction of action. These include static loads, dynamic loads, impact loads, cyclic loads, and temperature loads, describing the mechanical state of the structure under actual working conditions and forming the basis for structural optimization and analysis. The unfailable region is the area within the design domain that must be protected from damage, failure, or material removal. This region typically includes areas where force transmission paths must remain intact, critical load-bearing components, assembly connections, and safety-critical parts that cannot fail under dynamic loads. Designating this region as unfailable ensures that the structure maintains necessary safety redundancy and functional integrity even in the event of local failures or material reduction due to topology optimization. Topology optimization determines which regions within the design domain retain material and which remove material, ultimately forming the optimal structural form. Dynamic load response refers to the response of a structure to displacement, velocity, acceleration, stress, strain, etc., as a function of time under time-varying loads. Common dynamic loads include periodic loads (such as sinusoidal loads), random vibrations, impact loads, and fatigue loads. Dynamic response analysis is usually more complex than static analysis, requiring consideration of the effects of changes in mass, damping, and stiffness over time on structural behavior. Failure-safe topology optimization is a method that introduces redundant design and local failure considerations into topology optimization, enabling the structure to maintain necessary load-bearing capacity even after damage occurs in certain parts. Its core ideas include simulating local material failure (such as removing parts of the structure) during optimization, ensuring that the structure still meets stiffness, strength, or stability requirements under such failure states, and obtaining an optimal structural layout with "fault tolerance."

[0078] Specifically, step 1 includes:

[0079] Step 1.1: Based on the load-bearing form, main deformation mode and typical stress characteristics of the fiber-reinforced composite material structure to be optimized in practical applications, the three-point bending beam is selected as the structure to be optimized.

[0080] This step analyzes the target structure, which primarily bears load through bending deformation in practical use and is prone to typical composite material failure modes such as bending fatigue and interlaminar shear failure under dynamic load conditions. Therefore, a structure capable of characterizing bending-dominated forces needs to be selected as the optimization target. The three-point bending beam is a standard structural type in materials mechanics and structural mechanics used to characterize the bending resistance and force transmission path characteristics of materials, generating a clear bending moment distribution, shear force distribution, and local stress concentration. The three-point bending beam structure can significantly amplify the influence of fiber orientation differences on structural performance, making it suitable for verifying the synergistic effect of fiber orientation optimization and material distribution optimization. Therefore, the three-point bending beam is selected as the basic structure for the fiber-reinforced composite material structure to be optimized in this application.

[0081] Step 1.2: Based on the geometric and force symmetry of the three-point bending beam, the right half of the three-point bending beam is selected as the design domain to construct a semi-three-point bending beam model.

[0082] The loading point of the three-point bending beam is located at the center of the structure. The load-bearing paths, deflection distribution, and stress distribution on both sides are centrally symmetrical about the central axis, thus simplifying the model to half its original size without altering the structural stress mechanism. Since topology optimization requires repeated solving of the finite element dynamic response equations, resulting in significant computational overhead, the structural symmetry allows the computational domain to be reduced to half its original size, effectively reducing the computational burden and improving optimization efficiency. This application extracts the right half of the three-point bending beam, retaining the typical loading path, boundary constraints, and structural deformation characteristics of this half, and sets necessary symmetry constraints based on the symmetry relationship. Therefore, the right half of the three-point bending beam is used as the design domain of this application, i.e., the semi-three-point bending beam structural domain.

[0083] The design domain area of ​​this application for:

[0084] (1)

[0085] In the formula: L is the length of the design domain; H is the height of the design domain; x and y are spatial coordinate variables in the structural coordinate system where the design domain is located, used to describe the position of any point in the design domain, where x is the horizontal coordinate (along the beam length direction), x=0 corresponds to the left end face, x=L corresponds to the right end face, and y is the vertical coordinate (beam height direction), describing the position of each point in the beam thickness or beam height direction.

[0086] Step 1.3: Determine the design parameters of the design domain based on its geometric range, dynamic load application mode, boundary constraint arrangement, and failure safety design requirements.

[0087] Step 1.3 involves determining the design parameters of the design domain, which serve as the input conditions for topology optimization and fiber orientation optimization. The design parameters of the design domain include geometric dimensions, dynamic load parameters, boundary constraint parameters, and parameters of non-failure regions.

[0088] Engineering practice and structural mechanics research generally agree that when the span-to-depth ratio L / H ≥ 3, the main deformation mode of the structure is bending deformation, with a relatively small shear deformation component; when the span-to-depth ratio L / H < 3, the proportion of shear deformation increases, no longer conforming to the typical stress mode of a bending beam. Therefore, choosing a length-to-width ratio of approximately 3:1 ensures that the structure mainly exhibits a bending-dominated response, clear and stable bending moment and stress distribution patterns, and displays typical bending dynamic characteristics under dynamic loads. Therefore, the geometric dimensions of the design domain in this application are determined to be a rectangular area with a length-to-width ratio of 3:1, for example, a rectangular design domain with a length of 180 and a width of 60, based on the span-to-depth ratio requirement, stress distribution range, and engineering scale of a semi-three-point bending beam, ensuring that the design domain can cover the complete stress path and the free optimization region. (See attached specification) Figure 3 As shown, Figure 3 The embodiment of this application shows that the fiber-reinforced composite material structure adopts a semi-three-point bending beam structure as the optimization analysis object.

[0089] Dynamic load parameters include dynamic load type (sinusoidal periodic load), dynamic load amplitude, dynamic load frequency and period, dynamic load direction, and dynamic load location (top left corner or node) to induce typical bending response and local stress fluctuations. Load types for bending structures include, but are not limited to, periodic dynamic loads (such as periodic vibration of rotating components), impact loads (such as external impacts, start-stop impacts), random loads (such as wind loads, environmental excitations), cyclic fatigue loads, static dead loads, or static operating loads. In this embodiment, to simulate the typical behavior of a composite beam in a vibration environment, the load type is determined to be a time-varying periodic sinusoidal dynamic load:

[0090] (2)

[0091] In the formula: The magnitude of the dynamic load applied to the structure to be optimized; t is a sine function used to describe the periodic load variation characteristics; t is a time variable used to express the load variation over time. The angular frequency coefficient of the sine function; Let be the period parameter of the load, indicating the period of action of the dynamic load. Second; The phase angle of the dynamic load reflects the periodic characteristics of the applied load changing over time.

[0092] The direction of load application is determined based on mechanical properties: Combining the bending characteristics and stress distribution law of composite beam structures, the dynamic load is set to act vertically downward to induce significant bending deformation and principal stress changes, thereby more effectively driving fiber orientation optimization.

[0093] Determining the load application location based on structural function: By analyzing the force transmission path, key load-bearing components, and stress concentration areas of the structure, it is determined that the dynamic load should be applied at the location that can produce a typical bending response. To simulate the typical working state of a semi-three-point bending beam under dynamic load, this embodiment sets the following stress conditions: 1. Apply a fixed constraint at the lower right corner of the design domain; 2. Apply a time-varying dynamic load at the upper left corner of the design domain.

[0094] Specifically, the lower right corner is one of the beam's support points. According to the definition of a semi-three-point bending beam, its load-bearing method is that one end is fixed while the other end is subjected to a dynamic load. In this application, the basis for choosing to apply a fixed constraint at the lower right corner of the design domain includes: 1. The fixed constraint restricts both rotation and displacement at this location, forming a significant bending stress concentration area and a typical bending failure mode. 2. The fixed support can serve as a defined boundary for topology optimization, preventing accidental deletion by the optimization algorithm and satisfying the requirement for clear boundaries in topology optimization. 3. Composite beams in engineering often use end-fixed supports, consistent with the actual installation method of composite beams in engineering.

[0095] The boundary conditions for the lower right corner of the design domain, which is fixed, are as follows:

[0096] (3)

[0097] In the formula: The displacement vector of the fixed support at the lower right corner; The velocity vector is the fixed support at the lower right corner.

[0098] The upper left corner of the design domain conforms to the characteristic of "loading point biased to one end" in a semi-three-point bending beam, and it adopts a sinusoidal time-varying load. In this application, the rationale for applying a time-varying dynamic load at the upper left corner of the design domain includes: 1. Typical vibration sources include rotating machinery vibration, periodic impact, and wind-induced vibration; sinusoidal dynamic loads can reflect common vibration patterns in engineering structures. 2. The upper left corner is a position where the span deviates from the support point, allowing for significant bending and local stress peaks; loading at the upper left corner can create a unilateral bending effect. 3. Composite material orientation design is highly sensitive to dynamic loads; this loading method can fully demonstrate the optimization advantages and thoroughly test the performance of composite material orientation optimization under dynamic conditions.

[0099] Based on the unidirectional bending characteristics of the three-point bending beam structure, lateral movement is restricted. Therefore, a horizontal degree-of-freedom constraint is applied to the left end face of the design domain to prevent lateral drift of the three-point bending beam under dynamic load. This ensures that the main response caused by the dynamic load is concentrated in the vertical bending direction, which helps control the problem dimension and stabilizes the optimization process. Therefore, the boundary constraint parameter of the design domain is a horizontal degree-of-freedom constraint applied to the left end face.

[0100] The horizontal constraint on the left end face of the design domain is:

[0101] (4)

[0102] In the formula: Horizontal movement is prohibited on the left end face, but vertical movement is permitted.

[0103] Because a dynamic load F is applied to the upper left corner of the left boundary of the design domain, and a single-point fixation is applied to the lower right corner of the right boundary, the left and right ends are the constraint and load application locations, belonging to the critical load-bearing areas. Failure in these areas will lead to overall structural instability. Therefore, areas of a preset unit width (e.g., 30 units wide) on both sides of the design domain are designated as non-failure-prone areas, as shown in the appendix to the specification. Figure 3 As shown in the blue area, the left and right boundary regions will not fail in the failure-safe topology optimization of fiber-reinforced composite materials. This ensures that key functional areas will not be removed by the material during the topology optimization process, and avoids dangerous designs caused by the optimization algorithm in pursuit of lightweighting.

[0104] The non-failable region (e.g., unit width b=30) is:

[0105] (5)

[0106] The material density of the non-failure region is fixed and it does not participate in topology optimization.

[0107] (6)

[0108] In the formula: The matrix material density remains constant in the non-failure region or the initial design domain.

[0109] The material property field of the design domain is:

[0110] (7)

[0111] In the formula: The fiber direction is the initial angle. The directional stiffness matrix of the composite material.

[0112] The original topological model of the design domain is represented by finite element discretization:

[0113] (8)

[0114] In the formula: The design domain is the overall finite element discrete region, which contains finite element solution elements such as nodes, elements, material properties, boundaries, loads, and degrees of freedom. The region of the e-th finite element element; To design domain The total number of finite element elements obtained after mesh generation, for example: design domain 180×60, element size 1×1, then ; For all units When pieced together, they form the overall design domain. .

[0115] It should be noted that formula (1) defines a continuous geometric region (geometric domain). The design domain is a continuous rectangular region in a two-dimensional coordinate system, which is the initial geometric space range used for topology optimization, material optimization, load application, etc. It has no mesh and no elements, and is a region in a mathematical sense. Formula (8) defines the set of finite element elements (discrete domain) after discretizing the continuous geometric region. Divided into several finite element elements, each element is All units The set is recombined to obtain the discrete domain in the finite element sense, thus forming the original topological model for dynamic analysis and topology optimization.

[0116] Step 2: Based on the design domain constructed in Step 1, perform dynamic mechanical analysis on the fiber-reinforced composite structure and obtain the equivalent static load for topology optimization.

[0117] The purpose of step 2 is to transform the time-varying dynamic load into an equivalent static load that can be used for static topology optimization by solving the transient dynamic response of the structure under dynamic loads. This allows the subsequent optimization process to achieve topology optimization design of the structure in a static form while preserving the influence of dynamic characteristics.

[0118] Step 2 includes:

[0119] Step 2.1: Based on the geometric dimensions, boundary constraint arrangement, material properties and mesh generation of the design domain determined in Step 1, the design domain is discretized by finite element method, and the mass matrix M, damping matrix C and stiffness matrix K of the structure are generated to construct the transient dynamic control equations.

[0120] The transient dynamic governing equations for fiber-reinforced composite structures are:

[0121] (9)

[0122] In the formula: The structural mass matrix; Here is the structural damping matrix; Here is the stiffness matrix of the structure; It is the displacement vector; It is the velocity vector; It is the acceleration vector; Let t be the dynamic load on the structure at time t.

[0123] Structural mass matrix M, structural damping matrix Stiffness matrix of the structure All of these are related to design variables and change with the design variables. The design domain is a 3:1 rectangular area, such as 180×60. The size determines the number of meshes, and the number of meshes determines the DOF (degrees of freedom). The dimensions (matrix order) of M, C, and K are directly determined by the number of DOFs. The element length and width affect the element stiffness matrix k. e Mass matrix m e Therefore, the values ​​of the overall K and M matrices after assembly will be affected. Since the boundary constraints in step 1 define the lower right corner fixed support (vertical / horizontal / rotational full constraints) and the left end face horizontal degree of freedom constraint, these constraints determine that the constrained degrees of freedom of the stiffness matrix K of the structure must be eliminated or boundary conditions must be applied, affecting the dynamic modal characteristics and transient response. The material properties of fiber-reinforced composite materials: material density ρ, elastic modulus E, shear modulus G, and layup direction, these affect the element local stiffness matrix k(θ), element mass matrix m, and overall damping matrix C (if proportional damping is used, C=αM+βK), so the material parameters determine the values ​​of M, C, and K. In addition, the design domain dynamic load in step 1 is applied to the "upper left corner node", the direction is vertically downward, and the load type is a periodic sinusoidal dynamic load that varies with time, all of which will determine the vector dimension and degree of freedom of action of F(t).

[0124] For fiber-reinforced composite structures discretized using the finite element method, the overall mass matrix M can be expressed as the mass matrix m of each element. e The set of elements, whose expression is:

[0125] (10)

[0126] In the formula: The number of grid cells in the design domain in step 1 is determined, for example, by the design domain size of 180×60 and the cell size; The unit mass matrix depends on the material density. , unit length Unit area And beam types.

[0127] In formula (10), a two-dimensional beam element (Euler-Bernoulli) is used, and the element mass matrix is... This can be expressed as:

[0128] (11)

[0129] In the formula: The density of the material is determined by the material properties of the composite material in step 1; The cross-sectional area or equivalent area of ​​the unit is determined by the unit size, thickness, and thickness of the laminated structure. The element length is determined by the design domain size and mesh generation. The contribution of inertia related to the node itself to the element mass; This refers to the coupling inertial effect of the internal mass of the unit on the two nodes.

[0130] For example: The mass distribution of the beam element is continuous along its length, therefore each node bears approximately [a certain percentage] of the mass. The total mass. The mass at a certain point inside the beam will affect the inertial coupling between the left and right nodes, thus causing... .

[0131] The overall stiffness matrix K is derived from the stiffness matrices of each element. Obtained by assembling unit by unit:

[0132]

[0133] In the formula: The number of grid cells in the design domain in step 1 is determined, for example, by the design domain size of 180×60 and the cell size; The element stiffness matrix is ​​formed by the material stiffness and geometry, and is a constituent element of the overall K.

[0134] In formula (12), the element stiffness matrix The directional stiffness of the corresponding fiber-reinforced composite material is expressed as:

[0135]

[0136] In the formula: is the spatial region of the e-th finite element element, used for the integration required to calculate the element stiffness; The strain-displacement matrix of the element is obtained by differentiating the element shape function N with respect to spatial coordinates. It is used to convert nodal displacements into strains within the element. It is the transpose of the strain-displacement matrix B of the element, used to transfer the relationship between material stiffness and strain back to the displacement space; For fiber-reinforced composites at fiber orientation angle The directional stiffness matrix is ​​derived from the material constitutive matrix. and coordinate transformation matrix We obtain the mechanical properties used to describe the anisotropy of composite materials; This represents the strain energy density generated by the element under a unit displacement. In the cell region The energy term above Integrating, we obtain the stiffness matrix of the corresponding element. .

[0137] Formula (13) yields the element stiffness by combining the strain matrix, the material directional stiffness matrix, and the derivative of the shape function, and by integrating the strain energy density over the element region to obtain the overall resistance to deformation of the element. This stiffness matrix also reflects the anisotropic properties of the composite material as the fiber direction changes.

[0138] In formula (13), the directional stiffness matrix of the composite material is... for:

[0139]

[0140] In the formula: Transformation matrix The transpose of the matrix; The constitutive matrix of the material is determined by the material parameters: elastic modulus in the fiber direction, elastic modulus in the matrix direction, shear modulus, and Poisson's ratio. This is the coordinate transformation matrix, which rotates the material coordinate system. To the structural coordinate system; This represents the rotation angle between the material coordinate system and the structural coordinate system, such as ±45°, 0°, 90°, etc.

[0141] The damping matrix C adopts the proportional damping (Rayleigh damping) form:

[0142] (15)

[0143] In the formula: α and β are damping proportionality coefficients, respectively; M is the structural mass matrix; K is the structural stiffness matrix.

[0144] In summary, combining steps 1 and 2.1, the matrix generation process is as follows: Mesh the matrix according to the design domain size (e.g., 180×60) to obtain... For each element, generate an element mass matrix based on the material parameters (E, G, ρ, θ). Element stiffness matrix Element damping matrix Apply the boundary constraints defined in step 1 to K, M, and C: delete the rows and columns corresponding to the fixed nodes and restrict the horizontal degrees of freedom on the left.

[0145] Step 2.2: Based on the transient dynamic control equations, solve for the displacement vector of the structure at time m under dynamic load. .

[0146] The Newmark-β method, typical of structural dynamics, is adopted, which discretizes the time domain [0,T] into N (N=T / T is the load period in step 1, and the time variable in step 1 is t = m × ) time steps:

[0147] (16)

[0148] In the formula: The time spent walking is long; This represents the m-th time point; m is the time step index, m∈0,1,2,3,...,N. For example, T=1 second. =0.001 seconds, then N=1000 time steps, each corresponding to a different value of m.

[0149] Assuming the structure is initially in static equilibrium, with initial displacement, initial velocity, and initial acceleration all taken as zero, i.e. , , During the time integration process, the predicted displacement and predicted velocity at step m are first calculated, and then discretely solved using the Newmark-β time integration method. The displacement discrete update formula is as follows:

[0150]

[0151] (17)

[0152] The predicted displacement is: (18)

[0153] Furthermore, the displacement correction formula is obtained as follows: (19)

[0154] The velocity discrete update formula is:

[0155]

[0156] (20)

[0157] The prediction speed is: (twenty one)

[0158] Furthermore, the velocity correction formula is obtained as follows: (twenty two)

[0159] Substitute equations (18) and (20) discretized into the dynamic equation. +C +K = ,get:

[0160] K

[0161] C

[0162] Adding the two items above together, we get:

[0163] (twenty three)

[0164] All containing Move the terms to the left, and keep the remaining known terms on the right, resulting in:

[0165] (twenty four)

[0166] The equivalent stiffness matrix is: (25)

[0167] The acceleration at time m is calculated using formulas (24) and (25). :

[0168] (26)

[0169] Furthermore, using the correction formulas (19) and (22), the true displacement at time m is obtained. With speed The displacement vector obtained from the structural dynamics solution. The degree of freedom ordering, node coordinate system, and numbering method are consistent with the design domain finite element model established in step 1. Therefore, the actual structural displacement vector at time m is... This is the nodal displacement vector used for subsequent equivalent static load calculations. Therefore .

[0170] Step 2.3: Based on the displacement at time m Calculate the equivalent static load at that moment.

[0171] Equivalent static load of dynamic load at time m for:

[0172] (27)

[0173] In the formula: Let be the equivalent static load vector at time m, which represents the dynamic response of the structure at time step m. Under the action of dynamic displacement, the static load form is equivalent to the dynamic displacement. It can convert the displacement caused by dynamic load into a static load form that can be directly used in topology optimization, so that the topology optimization process can reflect dynamic effects in a static framework. The superscript m is the mth time step or the mth moment, and the subscript eq is the abbreviation of equivalent, which means equivalent and is used to refer to the static load that is equivalent to the dynamic displacement response. This is the stiffness matrix of the structure, i.e., the overall stiffness matrix of the structure under finite element discretization; Let be the nodal displacement vector of the structure at time m, obtained by solving Newmark-β in step 2.2. .

[0174] Step 2.4: Equivalent static load at the m-th time step obtained in Step 2.3 The equivalent static load data for n time steps are obtained, and the n equivalent static load vectors are weighted, summed, or averaged to obtain the final equivalent static load for topology optimization. .

[0175] Final equivalent static load for:

[0176] (28)

[0177] In the formula: This is the final equivalent static load; Let be the equivalent static load vector at time step m; K is the overall structural stiffness matrix. Let be the nodal displacement vector of the structure at time step m; n is the number of time steps used to construct the final equivalent static load.

[0178] Based on steps 1 and 2.1, construct the original topology model of the design domain:

[0179] (29)

[0180] In the formula: This is the initial design domain; For displacement boundaries (fixed supports and horizontal constraints); The load boundary (dynamic load at the top left corner). These are the mass, damping, and stiffness matrices, respectively. This is the initial material stiffness matrix (initial fiber direction).

[0181] Based on formula (29), the obtained equivalent static load It is applied to the original topology model of the design domain for subsequent failure-safe topology optimization design.

[0182] Step 3: Replace the dynamic load of the original structure with the equivalent static load. Take the model with the equivalent static load as the analysis object. Use the material density variable, angle sub-interval selection variable and fiber angle variable in the above model as design variables. Construct multiple local failure conditions and, based on the equivalent directional stiffness matrix that is continuously related to the fiber angle and the normalized KS objective function, combine the moving asymptote method to perform failure-safe topology optimization iteration on the design variables until the inner loop converges.

[0183] The purpose of step 3 is to transform the time history response of the structure under dynamic loads into an equivalent static load, thereby converting the complex working condition that originally required dynamic nonlinear solution into a static optimization problem. This allows for the failure-safe topology optimization design of fiber-reinforced composite structures under local failure scenarios while ensuring dynamic consistency and equivalence. Through topology optimization driven by equivalent static loads, the structure can maintain the necessary load-bearing capacity and safety redundancy after local damage occurs, achieving an optimal fiber-reinforced composite structure layout with fault tolerance, resistance to local failure, and minimum flexibility.

[0184] Step 3 includes:

[0185] Step 3.1: Establish a failure safety topology optimization model for fiber-reinforced composite materials with material density variable, angle sub-interval selection variable and fiber angle variable as design variables, and use the maximum value of the compliance under each damage condition as the optimization objective.

[0186] This optimization model is used to simultaneously solve for material density variables, angle sub-interval selection variables, and fiber layup angle variables. The failure-safe topology optimization model for fiber-reinforced composite materials is constructed under the following constraints:

[0187] (30)

[0188] In the formula: The set of variables that the solver needs to solve; The density of the material; Choose variables for the angle sub-intervals, n sub-intervals; For fiber angle variables; To minimize the objective function f in all m damage scenarios and improve the stiffness under the weakest condition; For the first The flexibility of a damage scenario is expressed as follows: the greater the flexibility, the worse the stiffness; m represents the number of damage scenarios. To find the worst softness in all damage scenarios; In the first Damage scenario, displacement Mechanical equilibrium must be satisfied; This is the equivalent static load; Here is the stiffness matrix for the l-th damage scenario; The displacement vector is obtained under the l-th damage scenario; Material usage (volume fraction) constraints are used to limit the total amount of material used across the entire design domain; The sum of all cell densities. The larger the value, the more material is used; N is the total number of elements (e.g., the number of meshes); f is the maximum allowed material volume fraction. The variable chosen for the angle interval is also in [0,1]. To limit the fiber laying angle Scope; The fiber orientation angle determines the directional stiffness matrix D of the composite material. ); To limit the range of angle variation, the larger n is, the smaller the angle range, and the more precise the control of fiber direction.

[0189] Ensure that the variable being solved is the optimal value. Ensure that the structure retains sufficient stiffness even under the most severe damage conditions, thus achieving fault-tolerant design. Ensuring that every damage scenario satisfies physical and mechanical equilibrium during the optimization process is the core constraint of topology optimization. Achieve lightweight design and avoid full-material structure. Ensure that all design variables fall within a legal and continuously optimizable range. Ensure the fiber orientation conforms to engineering manufacturing constraints and avoid excessive angular jumps. By solving the above optimization model, the material distribution, fiber angle range, and fiber laying direction that meet the preset failure safety requirements are obtained.

[0190] Step 3.2: Based on the finite element mesh, the candidate local damage blocks are laid out and translated according to the preset mesh size to generate multiple local failure regions. The finite element elements inside and outside each region are assigned the minimum weakening coefficient and 1 respectively to form the corresponding weakening coefficient field, so as to construct the local failure scenario under multiple damage conditions required in Step 3.1.

[0191] In step 3.1, for each of the multiple damage scenarios that need to be considered, a corresponding local failure region is automatically generated. With weakening coefficient field This allows the material loss of fiber-reinforced composite structures under different local damage conditions to be introduced into the failure-safe topology optimization model in step 3.1 in the form of a unified "local failure scenario".

[0192] Specifically, to avoid confusion with time step m in step 2, this step uses subscripts. Indicates the first Each local failure scenario. The weakening coefficient acting on the e-th finite element is defined as:

[0193] (31)

[0194] In the formula: The design domain is given in step 1; For the first There are local failure regions; e is the e-th finite element element within the design domain; For this unit in a partial failure scenario The weakening coefficient below; As the preset minimum weakening coefficient, this embodiment takes... .

[0195] As can be seen from formula (31), when element e is located outside the local failure region in the design domain, This indicates that the material at that location is intact and undamaged; when element e is located in a local failure region... At that time, This value is a local minimum, which in finite element analysis is equivalent to the material stiffness being approximately zero at that point, meaning the material is completely weakened. This corresponds to local failure in that region, leading to material loss. The physical meaning of this local failure is illustrated as follows: Figure 4 As shown, Figure 4 This diagram illustrates the distribution of a single local failure region and weakening coefficient in a fiber-reinforced composite structure. Within the design domain, material loss occurs only in the local area within the red box, while the material remains intact in the remaining areas. A set of weakening coefficients can be used with formula (31). To characterize the first Under what specific damage conditions will the failure occur at a particular location and to what extent?

[0196] Furthermore, to systematically and uniformly consider potential local damage locations within the design domain, the damage group method is used to classify the local failure scenarios of the semi-three-point bending beam: First, the preset size of the local damage is determined, and the preset side length of the local failure region is set as L, that is, each local failure region... The design domain is a square region of size L×L, where the value of L is selected based on the design domain size and the finite element mesh size. This ensures that a localized damage block contains several finite element elements, reflecting the localized damage without being overly fragmented. Then, within the design domain... Within the area (excluding the non-failable region defined in step 1), square damage blocks are arranged in a tiling manner: the square damage blocks do not overlap, and are laid row by row and column by column in the horizontal and vertical directions with grid alignment, until the entire potential failure region is covered. Each square damage block corresponds to a candidate local failure region. The overall schematic diagram of the locally damaged area obtained by the above-mentioned flat tiling is shown in the attached instruction manual. Figure 5 As shown in the left figure, Figure 5 A schematic diagram of the division and translation of local failure regions in the design domain of a three-point bending beam (schematic diagram of multi-damage conditions). Figure 5 In the left image, the square damage subdomain fills the entire design domain, while the long rectangular area remains the entire design domain. The design domain is composed of a checkerboard of small squares of different colors. Each small square represents a candidate local damage area with a size of L×L. These small squares do not overlap with each other and fill the entire design domain in the longitudinal and transverse directions. During the initial division, each colored square can be considered to correspond to the position of a "damage group". In subsequent optimization, one or more of them can be selected as the actual local failure conditions.

[0197] To consider more potential local damage locations without increasing the number of meshes, the damaged blocks are translated to generate additional candidate failure regions based on the initial tiling: First, the square damage blocks are translated diagonally, with a translation step of L / 2 in both the horizontal and vertical directions, meaning the center point of the damage block moves by L / 2 in the x and y directions respectively. Then, the translated regions are screened. For each translation position, a new set of candidate local failure regions is formed. Local damage blocks that exceed the design domain boundary or overlap with non-failable regions are deleted, and all damage blocks located within the failable regions are retained as new candidate failure regions. After the above translation and filtering processes, a series of denser local damage candidate regions covering the design domain are obtained. The division and distribution of these regions are illustrated in the appendix to the specification. Figure 5 As shown in the right figure.

[0198] Step 3.3: Divide the original design interval of the fiber direction into several sub-intervals according to parameter n, calculate the directional elastic matrix corresponding to the central angle of each sub-interval, and introduce discrete-continuous weight variables and small deflection angles within the sub-intervals. Weight and rotate the directional elastic matrix to obtain an equivalent directional stiffness matrix that is continuously related to the fiber angle design variable. Substitute the equivalent directional stiffness matrix into the overall stiffness matrix expression to form the stiffness matrix used in Step 3.1.

[0199] Step 3.3 introduces the fiber angle design variable from Step 3.1 into the directional stiffness matrix of the finite element, allowing the fiber orientation to vary continuously within a given angle range while satisfying the manufacturing constraint that the angle can only take values ​​in a finite number of sub-intervals. This enables the fiber orientation of each element to be updated in a differentiable manner during topology optimization iterations. This step provides the fiber angle variable and the directional stiffness matrix D(…). The specific mapping relationship between the element stiffness matrix and the element stiffness matrix in step 2 is determined. Connected.

[0200] Specifically, firstly, the original design range of the fiber direction... Divide the data into n sub-intervals according to the parameter n in step 3.1, with each sub-interval having a width of n. The central angle of the i-th subinterval is denoted as . Its expression is:

[0201] (32)

[0202] In formula (32), each For a typical layup angle (such as -45°, -30°, ..., 45°), the fiber direction of any unit in the structure is restricted to a small offset range near these typical angles, thereby ensuring that the fiber layup angle meets the process discretization requirements.

[0203] In the principal direction coordinate system of the material, the initial elastic matrix of the composite laminate is denoted as: :

[0204] (33)

[0205] Corresponding central angle stress-strain transformation matrix for:

[0206] (34)

[0207] Then, in the structural coordinate system, the directional elasticity matrix corresponding to the i-th central angle is... for:

[0208] (35)

[0209] Each Characterizing the strict orientation of fibers The stiffness characteristics of composite materials at that time.

[0210] Furthermore, to select a subinterval in a differentiable manner during optimization, n discrete-continuous weight variables are introduced. And its corresponding weight function. Where, sub-weights... Defined as:

[0211] (36)

[0212] In the formula: The design variable associated with the i-th sub-interval can take continuous values ​​in [0,1]; q is a penalty parameter that slowly increases from 1 to 3 every few iterations during the optimization process in order to gradually approach the 0 / 1 solution.

[0213] Normalized weights for:

[0214] (37)

[0215] Normalized weights satisfy:

[0216]

[0217] After the penalty parameter q is increased sufficiently, the optimization result will cause a certain ,the remaining This achieves the discrete effect of "selecting 1 from n sub-intervals", but it is still continuously differentiable for optimization algorithms.

[0218] Furthermore, by weighting and summing the elasticity matrices of each central angle, we obtain the formula for the equivalent directional stiffness matrix after discrete-continuous parameterization:

[0219] (38)

[0220] at this time, This is equivalent to selecting a principal direction stiffness matrix from n candidate central angles.

[0221] Based on the selected master sub-interval, to further improve the precision of fiber orientation optimization, continuous offset angles within the sub-interval are introduced. Its value range is consistent with the angle constraint in step 3.1:

[0222]

[0223] This reflects a small deflection of the fiber direction relative to a certain central angle, representing continuous fine-tuning within a sub-interval. The corresponding rotation matrix is ​​denoted as... Form and They are the same, just from different angles.

[0224] Finally, the directional elasticity matrix at the location of the element. Expressed as:

[0225] (39)

[0226] At this point, the actual angle of the fiber can be understood as a certain central angle. Add offset angle ,and It is the equivalent stiffness matrix corresponding to that actual angle.

[0227] In step 2, the element stiffness matrix of formula (13) has been written as follows: ,in The fiber orientation angle is The directional stiffness matrix at time. Formula (39) obtained in this step. Substitute the values ​​to obtain the element stiffness matrix, which depends on the fiber angle design variables. Then Substituting into formula (12), the overall stiffness matrix K used in step 3.1 is finally obtained for failure-safe topology optimization. This realizes the mapping from "design variables to element stiffness matrices", providing updated stiffness information for each topology optimization iteration.

[0228] Step 3.4: Based on the failure-safe topology optimization model established in Step 3.1 with the maximum flexibility of each local failure condition as the objective, and according to the structural flexibility under each local failure scenario obtained in Step 3.2, the maximum flexibility objective is transformed into a continuously differentiable objective function using the normalized KS aggregation function.

[0229] This step, based on the failure-safe topology optimization model established in step 3.1, which aims to maximize the compliance of each local failure condition, will further refine the target value. Transforming the non-smooth discrete form into a differentiable continuous function form facilitates the use of gradient-based optimization algorithms to solve subsequent topology optimization problems.

[0230] As described in step 3.2, the structure may experience various local failure scenarios, and corresponding compliance values ​​are obtained under m local failure conditions. In this embodiment, to improve the load-bearing capacity of the structure under the most unfavorable local failure conditions, the optimization objective is to minimize the maximum flexibility of the structure under all local failure conditions.

[0231] Under the l-th local failure condition, the equivalent static load obtained in step 2 is applied. The structural displacement vector is It satisfies the equilibrium equation:

[0232] (40)

[0233] In the formula: In the case of partial failure Below, based on the material density variable, attenuation coefficient field, and fiber angular directional elastic matrix. The overall stiffness matrix obtained by joint assembly.

[0234] Corresponding compliance is defined as:

[0235] (41)

[0236] Since the objective function f=max( (This belongs to a piecewise linear discrete function, in multiple...) Since the function is not differentiable when approximating the maximum value, it is difficult to solve directly through optimization. Therefore, it needs to be transformed into a continuously differentiable smooth function. To this end, the Kreisselmeier–Steinhauser (KS) aggregation function is introduced to approximate the maximum value. Using the normalized KS function, the objective function is written as:

[0237] (42)

[0238] In the formula: This is the objective function value after aggregation by the KS function. The smaller the value, the smaller the maximum compliance under each local failure scenario. is the natural logarithm function; m is the total number of local failure scenarios described in step 3.2; This represents the sequence number of the partial failure scenario. =1-m; For the first The flexibility of a structure under equivalent static load in a local failure scenario is defined as follows: a higher flexibility indicates a more flexible structure and poorer stiffness under that failure scenario; e is the base of the natural logarithm, representing the exponential function. This is different from the unit number symbol e in this application; This is the regularization parameter, a positive number, used to control how closely the KS function approximates the maximum compliance; The normalization coefficient is used to scale the compliance index term for each local failure scenario to improve numerical stability.

[0239] To ensure numerical stability while also considering the accuracy of maximum value approximation, this embodiment selects γ=30 / based on experience. ,in To match various flexibility A reference compliance value of the same order of magnitude could be used, for example, the compliance value under the first type of local failure scenario. Or the average value of compliance under several operating conditions.

[0240] Through the above processing, while maintaining the focus on the most unfavorable local failure condition, the originally difficult-to-differentiate maximum compliance objective can be transformed into a continuously differentiable KS objective function, providing a foundation for subsequent sensitivity analysis and optimization iteration.

[0241] Step 3.5: Based on the KS objective function constructed in Step 3.4 and the explicit derivative relationship of compliance with design variables, calculate the sensitivity of the objective function to material density variables, angle sub-interval selection variables, and fiber angle variables, and use the moving asymptote method to iteratively update the three types of design variables based on the sensitivity until the inner loop converges.

[0242] Based on the KS objective function obtained in step 3.4, the sensitivity expressions of the objective function to each design variable are derived, and the moving asymptote method (MMA) is used to evaluate the material density variable. Angle sub-interval selection variables and fiber angle variables The process is iteratively updated until the inner loop converges, thereby obtaining the optimal topology that meets the failure safety requirements.

[0243] Specifically, since step 3.4 defines the KS objective function Each of them It is the first The compliance of a partial failure condition. Each compliance condition... Another design variable: material density variable Angle sub-interval selection variables and fiber angle variables The function, because And stiffness matrix Due to material density variables Angle sub-interval selection variables and fiber angle variables Decide.

[0244] Furthermore, the sensitivity of the objective function f to three types of design variables, including material density variable, is given. The sensitivity is:

[0245] (43)

[0246] In the formula: Design variables for the material density of the e-th element for the objective function f. Sensitivity; The value of the objective function after aggregation by the KS function is smaller, indicating a smaller maximum compliance under each local failure scenario; m is the total number of local failure conditions. This is the sequence number of the partial failure condition; For the first The structural flexibility under equivalent static load in a partial failure condition; For softness For material density variables The partial derivatives; Design variables (relative density) for the volume distribution of the e-th unit. Right now The value is between 0 and 1, and is used to indicate whether the unit is equipped with fiber reinforcement material; This refers to the regularization parameter in the KS function; For the first Weighting factors corresponding to each local failure condition Softness The larger, the corresponding The larger the value, the greater its weight in the weighted sum; e is the base of the natural logarithm, representing the exponential function. ; To address all partial failure conditions ( =1,……,m) is the sign for summation.

[0247] Angle sub-interval selection variable The sensitivity is:

[0248] (44)

[0249] In the formula: Choose variables for the objective function f for the i-th angle sub-interval in the e-th unit. Sensitivity; The design variable associated with the i-th angle sub-interval of the e-th unit can take continuous values ​​in [0,1] and is used to characterize whether the unit selects the angle sub-interval. The value of the objective function after aggregation by the KS function is smaller, indicating a smaller maximum compliance under each local failure scenario; m is the total number of local failure conditions. This is the sequence number of the partial failure condition; For the first The structural flexibility under equivalent static load in a partial failure condition; For softness Choose variables for angle sub-intervals The partial derivatives; This refers to the regularization parameter in the KS function; For the first Weighting factors corresponding to each local failure condition Softness The larger, the corresponding The larger the value, the greater its weight in the weighted sum; e is the base of the natural logarithm, representing the exponential function. ; To address all partial failure conditions ( =1,……,m) is the sign for summation.

[0250] Unit fiber angle variable The sensitivity is:

[0251] (45)

[0252] In the formula: The objective function f is the fiber angle variable in the e-th element. Sensitivity; The fiber angle design variable (fiber orientation angle within the element) for the e-th finite element is varied through the directional stiffness matrix. Affects the stiffness of the unit and the whole; The value of the objective function after aggregation by the KS function is smaller, indicating a smaller maximum compliance under each local failure scenario; m is the total number of local failure conditions. This is the sequence number of the partial failure condition; For the first The structural flexibility under equivalent static load in a partial failure condition; In the first Compliance under partial failure conditions Design variables for unit fiber angle The partial derivative represents the effect of a small change in the fiber angle of the unit on the flexibility under this condition; This refers to the regularization parameter in the KS function; For the first Weighting factors corresponding to each local failure condition Softness The larger, the corresponding The larger the value, the greater its weight in the weighted sum; e is the base of the natural logarithm, representing the exponential function. ; To address all partial failure conditions ( =1,……,m) is the sign for summation.

[0253] In formulas (43)-(45), , For the first Weighting factors for local failure conditions The regularization parameters are given in step 3.4. For the first The compliance of a partial failure condition; the greater the compliance, the corresponding... The larger the value, the greater its proportion in the weighted sum, thus ensuring that the objective function focuses more on the maximum compliance condition.

[0254] In formula (43), the first Flexibility under partial failure conditions For unit volume distribution variables derivative for:

[0255] (46)

[0256] In the formula: For the first The structural flexibility under partial failure conditions; The material density design variable for the e-th element, with a value between 0 and 1, is used to indicate whether fiber-reinforced material is arranged in this element; In the first The column vector of nodal displacements related to element e under a local failure condition is the global displacement vector. Extraction of unit degrees of freedom; For the first The equivalent stiffness matrix of element e under local failure conditions; The element stiffness matrix, ; For the first The weakening coefficient of element e under a certain local failure scenario, with a value range of 0.001 or 1, is used to characterize whether the element fails in this failure scenario; The penalty exponent for volume distribution interpolation is used to enhance design variables. The 0-1 characteristic.

[0257] In formula (44), the flexibility Choose variables for angle sub-intervals derivative for:

[0258] (47)

[0259] In the formula: For the first Compliance under partial failure conditions Choose a variable for the i-th angle sub-interval in the e-th unit. The derivative; The design variable associated with the i-th angle sub-interval of the e-th unit can take continuous values ​​in [0,1] and is used to characterize whether the unit selects the angle sub-interval. In the first The column vector of nodal displacements under local failure conditions is the overall displacement vector. Extraction of the unit's degrees of freedom; For the first The equivalent stiffness matrix of element e under local failure conditions; Equivalent stiffness matrix Choose variables for angle sub-intervals The partial derivatives are used to characterize the design variables. The effect of variations on element stiffness and flexibility; subscript The subscript 'e' represents the local failure condition number; the subscript 'e' represents the element number; and the subscript 'i' represents the angle sub-interval number. For the first The weakening coefficient of element e under a certain local failure scenario; The equivalent material density of element e is determined by the material density design variable. After processing, the value range is [0,1], which is used to describe the proportion of material actually involved in the load-bearing of the unit; The penalty exponent for volume distribution interpolation; Let e ​​be the stiffness matrix of element e. Choose variables for angle sub-intervals The partial derivatives are used to characterize the effect of changes in the design variable on the element stiffness.

[0260] In formula (45), the flexibility For angle variables derivative for:

[0261] (48)

[0262] In the formula: For the first Compliance under partial failure conditions For the element fiber angle variable in the e-th element The derivative; Design variables for the fiber angle of the e-th finite element; In the first The column vector of nodal displacements under local failure conditions is the overall displacement vector. Extraction of the unit's degrees of freedom; For the first The equivalent stiffness matrix of element e under local failure conditions; Equivalent stiffness matrix For unit fiber angle variables The partial derivatives are used to characterize the design variables. The effect of variations on element stiffness and flexibility; subscript The subscript 'e' represents the local failure condition number; the subscript 'e' represents the element number; and the subscript 'i' represents the angle sub-interval number. For the first The weakening coefficient of element e under a certain local failure scenario; The equivalent material density of element e is determined by the material density design variable. After processing, the value range is [0,1], which is used to describe the proportion of material actually involved in the load-bearing of the unit; The penalty exponent for volume distribution interpolation; Let e ​​be the stiffness matrix of element e. For unit fiber angle variables The partial derivatives are used to characterize the effect of changes in the design variable on the element stiffness.

[0263] Furthermore, in formulas (44) and (45), and Further derivation:

[0264] (49)

[0265] In the formula: Let e ​​be the stiffness matrix of element e. Choose variables for angle sub-intervals The partial derivatives of, where Let be the stiffness matrix of the e-th finite element. Selecting a variable for the angle sub-interval corresponding to the i-th angle sub-interval on unit e is a design variable that controls which sub-interval the fiber orientation of the unit falls into. In the cell region The area element on the surface is used to represent the integration over that region; The element strain-displacement matrix is ​​obtained by differentiating the element shape function with respect to spatial coordinates, and is used to convert nodal displacements into strains within the element. for The transpose of the matrix; The equivalent directional stiffness matrix of element e Choose a variable for the i-th angle sub-interval The partial derivative matrix reflects the rate of change of material stiffness when the selected variable for this subinterval is changed, where... The equivalent directional stiffness matrix of element e is determined by considering the choice of variables in the angular sub-interval. It was obtained later.

[0266] (50)

[0267] In the formula: Let e ​​be the stiffness matrix of element e. For unit fiber angle variables The partial derivatives of, where Let be the stiffness matrix of the e-th finite element. Design variables for the fiber angle of unit e to describe the actual fiber laying direction in the unit; In the cell region The area element on the surface is used to represent the integration over that region; The element strain-displacement matrix is ​​obtained by differentiating the element shape function with respect to spatial coordinates, and is used to convert nodal displacements into strains within the element. for The transpose of the matrix; The equivalent directional stiffness matrix of element e For the unit fiber angle variable The partial derivative matrix reflects the effect of minute changes in fiber angle on the stiffness of the element material, where... The equivalent directional stiffness matrix of element e is determined by considering fiber angle variables. It was obtained later.

[0268] In formulas (49) and (50), and The following conclusions can be drawn:

[0269] (51)

[0270] In the formula: The equivalent directional stiffness matrix of element e Choose a variable for the i-th angle sub-interval The partial derivative matrix; Design variables for fiber angles of unit e. A defined stress-strain coordinate transformation matrix is ​​used to rotate the material coordinate system to the structural coordinate system; for The transpose of the matrix; Equivalent stiffness matrix in principal direction Angle sub-interval selection variable The partial derivative reflects the rate of change of the principal direction stiffness when the weight of the subinterval is changed, where The principal direction equivalent stiffness matrix is ​​obtained after discrete-continuous parameterization, and is determined by the central angle of each subinterval. The corresponding directional elasticity matrix By normalized weights The weighted average is obtained.

[0271] (52)

[0272] In the formula: The equivalent directional stiffness matrix of element e For the unit fiber angle variable The partial derivative matrix; Transformation matrix Design variables for fiber angle The derivative of is used to characterize the rate of change of the strain coordinate transformation relationship when the angle changes slightly, where Design variables for fiber angles of unit e. The determined stress-strain coordinate transformation matrix; Transpose transformation matrix Design variables for fiber angle The derivative of is used to characterize the rate of change of the stress coordinate transformation relationship when the angle changes slightly; is the principal direction equivalent stiffness matrix after discrete-continuous parameterization; for The transpose of the matrix; Because the stress coordinate transformation matrix varies Changes in stiffness caused by variations; Because the strain coordinate transformation matrix follows Changes in stiffness caused by variations.

[0273] In formula (51):

[0274] (53)

[0275] In the formula: Equivalent stiffness matrix in principal direction Angle sub-interval selection variable The partial derivatives; For normalized weights For design variables The partial derivatives describe the choice of variables when the angle subinterval is changed in cell e. At that time, the i-th weight The rate of change of, where The normalized weight coefficient for the i-th angle sub-interval (Formula 37); For the i-th central angle The corresponding directional elasticity matrix (Formula 35); n is the total number of angle sub-intervals obtained by equally dividing the original fiber angle interval [−π / 2,π / 2].

[0276] In formula (52):

[0277] (54)

[0278] In the formula: Transformation matrix Design variables for fiber angle The derivative; , Angle 2 respectively The sine and cosine functions are the transformation matrix of the original transformation matrix. (See Formula 34) , , The bi-angle trigonometric function terms that appear after differentiating the trigonometric terms.

[0279] The optimization iteration and convergence determination process of the Moving Asymptote Method (MMA) includes: Based on the KS objective function from step 3.4 and the sensitivity information given in this step, constructing local approximate subproblems of the MMA, incorporating volume constraints, etc., into the subproblems. Solving the subproblems yields updated design variables, and based on this, reassembling the stiffness matrix under each local failure condition, solving for the displacement vector, and then calculating the compliance and objective function values. Comparing the change in structural compliance between two adjacent iterations: if the relative change in structural compliance after this iteration is less than a preset threshold (e.g., 0.1%), applying the updated design variables to the original structure yields the optimized new structure, and recalculating the structural compliance. The structural compliance is... ,in For displacement, If the structural stiffness matrix is ​​related to the design variables, then the inner loop is considered to have converged, and the failure-safe topology optimization process for fiber-reinforced composite materials is terminated. If this condition is not met, the MMA iteration is continued with the new design variables as initial values ​​until the convergence criterion is met.

[0280] By using step 3.5 above, combined with the KS objective function processing in step 3.4, we can obtain the fiber-reinforced composite material topology with the minimum flexibility under the most unfavorable local failure condition, given a local failure scenario and equivalent static load conditions, thus realizing the failure safety design of the structure under local failure conditions.

[0281] Step 4: Determine whether the external loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load and the structural flexibility when constructing the equivalent static load is less than 0.1%. If it has converged, output the failure safety topology optimization result of the fiber-reinforced composite material at this time. Otherwise, perform dynamic analysis again on the updated structure and construct a new equivalent static load before returning to step 3 to continue optimization.

[0282] The purpose of step 4 is that after completing the inner loop topology optimization for the given equivalent static load in step 3, the dynamic characteristics of the structure change due to changes in the structural topology and stiffness distribution, rendering the original equivalent static load inaccurate. Step 4 aims to determine, based on the updated design variables after optimization, whether a new dynamic analysis and construction of a new equivalent static load are necessary. Through outer loop iteration, the error of "equivalent static load replacing dynamic load" is controlled within a preset range, resulting in a convergent failure-safe topology optimization result under the combined action of dynamic and static loads.

[0283] Specifically, the design variables (material density variables) obtained after the inner loop converges in step 3 are... Angle sub-interval selection variables and fiber angle variables Mapping this onto the original structural model, updating the element material stiffness, and obtaining the optimized new structure. This involves constructing the current equivalent static load. The corresponding compliance has been calculated on the "reference structure" used at that time. For the optimized new structure, under the same equivalent static load... Perform static analysis again under the action, solve for the displacement vector, and calculate the new compliance. :

[0284] (55)

[0285] In the formula: To optimize the flexibility of the structure under a given load; This is the column vector of nodal displacements of the optimized structure under a given equivalent load. displacement vector transpose; This is the overall stiffness matrix of the optimized structure.

[0286] Furthermore, the relative change in flexibility is calculated and compared with a preset convergence threshold (0.1% in this embodiment). The relative change in flexibility is:

[0287] (56)

[0288] In the formula: This represents the relative change in flexibility, expressed as a percentage. To optimize the flexibility of the initial structure; To optimize the flexibility of the structure.

[0289] External circulation convergence condition: If the change in flexibility satisfies This indicates that the reference structure used to construct the equivalent static load has very little difference in stiffness from the current optimized structure, and the equivalent static load is sufficiently accurate in replacing the actual dynamic load, thus confirming the convergence of the external circulation. If the flexibility change satisfies... If the structural topology changes have a significant impact on the dynamic characteristics, the original equivalent static load can no longer accurately represent the actual dynamic effect, the external circulation has not converged, and the equivalent static load needs to be updated and the iteration continued.

[0290] The dynamic-static coupling iterative process before convergence: when Then, using the new structure obtained from step 3 as the current structure, a new dynamic analysis is performed under actual dynamic loads (or a given excitation) to obtain updated response results; following the method described in step 2, a new equivalent static load is reconstructed based on the new dynamic response. (new) will (new) As a new external load, return to step 3 and perform failure-safe topology optimization on the fiber-reinforced composite structure again to complete a new round of internal loop iteration.

[0291] External loop termination condition and result output: When the compliance change is first satisfied Once the dynamic-static load coupling external loop converges, the iteration ends. The structural topology, local failure layout, and fiber angle distribution obtained at this time are the optimal fiber-reinforced composite material failure safety topology optimization results under the combined action of a given dynamic load and equivalent static load.

[0292] Example:

[0293] This embodiment uses the aforementioned failure safety topology optimization method for fiber-reinforced composite materials under dynamic load to optimize the design of a semi-three-point bending beam structure, thereby illustrating the effectiveness of the method of the present invention under local failure conditions.

[0294] First, the structural model and boundary conditions are set: the design domain is the right half of the three-point bending beam, a rectangular region with a length of 180 and a height of 60. The mesh is divided into 180×60 planar elements using 1×1 element sizes. Regions with a width of 30 at each end are designated as non-failure zones, with a material density of ρ0. Unidirectional fiber-reinforced composite material is used, with the initial fiber direction set to 0°. Its anisotropic elastic matrix is ​​given by the aforementioned formula. The loads and constraints are set as follows: a fixed constraint is applied at the lower right corner of the design domain, a horizontal displacement constraint is applied at the left end face, and a time-varying sinusoidal dynamic load F(t) = sin(2πt / 4×10⁻³) is applied at the upper left node, with the load direction vertically downwards.

[0295] Following step 2, solve for the transient response of the structure over one load cycle. Use the Newmark-β integral to obtain the displacement vector at each time step, and construct the equivalent static load vector F accordingly. eq The equivalent static load F eq The original topology model of the design domain is applied as an external load for subsequent static topology optimization. The material volume fraction constraint is set to 0.4, meaning that the material usage of the optimized structure does not exceed 40% of the design domain volume.

[0296] Furthermore, a fiber orientation design variable is set: To illustrate the influence of the number of discrete sub-intervals n in the fiber orientation, this embodiment compares n=1, 2, and 3 respectively. The original design interval for the fiber orientation is [-π / 2, π / 2], which is divided into n sub-intervals according to equation (32), and each sub-interval corresponds to a central angle. Using discrete-continuous weight variables And the deflection angle θ, construct the equivalent directional stiffness matrix The overall stiffness matrix K is obtained by assembling the matrix, and then the topology optimization iteration is completed.

[0297] The topology optimization results and analysis are as shown in the appendix to the instruction manual. Figure 6 , Figure 7 As shown. Figure 6 (a1)-(a3) show the topology optimization results under dynamic load when the number of subintervals n is 1, 2 and 3, respectively, without pre-set damage patches. Figure 6 Figures (b1)-(b3) show the results of failure-safe topology optimization for the same working condition after a damage patch of length L=30 is pre-placed on the upper edge of the beam. The direction of the short lines in the figures indicates the optimized local fiber orientation. To facilitate identification of fiber orientation in different sub-intervals, Figure 6 Different colors are used for marking: when n=1, there is only one sub-interval, and all fiber directions are represented by white; when n=2, white and yellow correspond to the direction intervals [-π / 2,0] and [0,π / 2] respectively; when n=3, three colors are used to represent the direction intervals [-π / 2,-π / 6], [-π / 6,π / 6] and [π / 6,π / 2] respectively.

[0298] Depend on Figure 6 As can be seen from (a1)-(a3), under the condition of a target volume fraction of 40% without considering pre-existing damage, the optimization results provide a clear main load-bearing path. (Comparison) Figure 6 (a1) (a2) Figure 6 In (a3), when n=1, the final fiber direction is approximately perpendicular to the direction of the web member; when n=2 or 3, the fiber direction tends to be arranged along the member direction, making the material's principal direction more consistent with the force transmission path of the bending problem. Correspondingly, the compliance value decreases from approximately 66.71 to approximately 62.37 and 62, indicating that increasing the number of subintervals can achieve a more reasonable fiber orientation and lower compliance.

[0299] With pre-applied damage patches Figure 6 As can be observed in (b1)-(b3), multiple redundant load paths are formed within the structure obtained by the failure-safe topology optimization. When local failure occurs, the load can bypass the damaged area through these alternative paths, thereby improving the overall resistance to failure. From the final compliance index, the compliance of schemes (b) and (c) is approximately 88.98 and 84.55 / 84.15, respectively, which is significantly better than the corresponding scheme (a). It can be determined that the structure obtained when n=1 is a locally optimal solution.

[0300] To further illustrate the effectiveness of the method of the present invention under localized sudden damage, in Figure 6 Based on the structures (a2) and (b2) in the diagram, a post-damage patch is applied at the upper edge of the mid-span to obtain... Figure 7 The results are shown. From Figure 7As can be seen, even under such severe damage conditions, the failure-safe structure obtained by this method still retains multiple load transfer paths, and the overall structure does not experience catastrophic failure. Using the flexibility increment as the evaluation index, when... Figure 6 The area applied to the structure of (a2) is approximately 30. 2 When applying a posterior damage patch, the flexibility value increases to approximately 2.15 × 10⁻⁶. 9 ; and Figure 6 When the same damage patch is applied to the structure in (b2), the compliance only increases to about 287.21, a difference of several orders of magnitude. This indicates that the failure-safe topology optimization method proposed in this invention can improve the safety margin and robustness of the structure under dynamic loads and local failure conditions without increasing the amount of material used.

[0301] In summary, this embodiment demonstrates that by considering fiber orientation optimization and the arrangement of damage patches, the failure safety topology optimization method for fiber-reinforced composite materials under dynamic load of the present invention can effectively suppress the sudden drop in load-bearing capacity caused by local damage, reduce the risk of the structure falling into local optima, and thus obtain a composite material topology with high safety and reliability.

[0302] This invention also provides a failure-safe topology optimization design system for fiber composites under dynamic loads, the system comprising:

[0303] The design parameter definition module is used to define the design parameters of the design domain of fiber-reinforced composite material structures based on the stress conditions and geometric characteristics of the structure to be optimized.

[0304] The Dynamic Analysis and Equivalent Static Load Generation Module is used to perform dynamic mechanical analysis on fiber-reinforced composite structures based on the design domain constructed by the Design Parameter Definition Module, and obtain the equivalent static load for topology optimization.

[0305] The failure-safe topology optimization iteration module is used to replace the dynamic load of the original structure with the equivalent static load obtained by the dynamic analysis and equivalent static load generation module. The model with the applied equivalent static load is used as the analysis object. The material density variable, angle sub-interval selection variable and fiber angle variable in the above model are used as design variables to construct multiple local failure conditions. Based on the equivalent directional stiffness matrix and normalized KS objective function that are continuously related to the fiber angle, the failure-safe topology optimization iteration is performed on the design variables in combination with the moving asymptote method until the inner loop converges.

[0306] The outer loop convergence judgment and result output module is used to determine whether the outer loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load obtained by the failure-safe topology optimization iteration module is less than 0.1%. When converged, the failure-safe topology optimization result of the fiber-reinforced composite material is output. Otherwise, the dynamic analysis is re-performed on the updated structure and a new equivalent static load is constructed before returning to the failure-safe topology optimization iteration module to continue optimization.

[0307] The above descriptions are merely embodiments of this application, and common knowledge regarding specific structures and characteristics in the solutions is not described in detail here. It will be apparent to those skilled in the art that this application is not limited to the details of the above exemplary embodiments, and that this application can be implemented in other specific forms without departing from the spirit or essential characteristics of this application. Therefore, the embodiments should be considered exemplary and non-limiting in all respects, and the scope of this application is defined by the appended claims rather than the foregoing description. Therefore, it is intended that all variations falling within the meaning and scope of equivalents of the claims be included within this application. No reference numerals in the claims should be construed as limiting the scope of the claims.

Claims

1. A method for topology optimization design to ensure failure safety of fiber composites under dynamic loads, characterized in that, The method includes: Step 1: Based on the stress conditions and geometric characteristics of the structure to be optimized, define the design parameters of the design domain for the fiber-reinforced composite structure; Step 2: Based on the design domain constructed in Step 1, perform dynamic mechanical analysis on the fiber-reinforced composite structure and obtain the equivalent static load for topology optimization; Step 3: Replace the dynamic load of the original structure with the equivalent static load. Take the model with the equivalent static load as the analysis object. Use the material density variable, angle sub-interval selection variable and fiber angle variable in the above model as design variables. Construct multiple local failure conditions and, based on the equivalent directional stiffness matrix and normalized KS objective function that are continuously related to the fiber angle, perform failure-safe topology optimization iteration on the design variables using the moving asymptote method until the inner loop converges. Step 4: Determine whether the external loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load and the structural flexibility when constructing the equivalent static load is less than 0.1%. If it has converged, output the failure safety topology optimization result of the fiber-reinforced composite material at this time. Otherwise, perform dynamic analysis again on the updated structure and construct a new equivalent static load before returning to step 3 to continue optimization.

2. The method for topology optimization design for failure safety of fiber composites under dynamic load as described in claim 1, characterized in that, Step 1 includes: Step 1.1: Based on the load-bearing form, main deformation mode and typical stress characteristics of the fiber-reinforced composite structure to be optimized in practical applications, a three-point bending beam is selected as the structure to be optimized. Step 1.2: Based on the geometric and force symmetry of the three-point bending beam, the right half of the three-point bending beam is selected as the design domain to construct a semi-three-point bending beam model; Step 1.3 determines the design parameters of the design domain to establish the input conditions for topology optimization and fiber orientation optimization.

3. The method for topology optimization design for failure safety of fiber composites under dynamic load as described in claim 1, characterized in that, Step 2 includes: Step 2.1: Based on the geometric dimensions, boundary constraint arrangement, material properties and mesh generation of the design domain determined in Step 1, the design domain is discretized by finite element method, and the mass matrix M, damping matrix C and stiffness matrix K of the structure are generated to construct the transient dynamic control equations. The transient dynamic governing equations for fiber-reinforced composite structures are: ； In the formula: The structural mass matrix; Here is the structural damping matrix; Here is the stiffness matrix of the structure; It is the displacement vector; It is the velocity vector; It is the acceleration vector; The dynamic load on the structure at time t; Step 2.2: Based on the transient dynamic control equations, solve for the displacement vector of the structure at time m under dynamic load. ; Step 2.3: Based on the displacement at time m Calculate the equivalent static load at that moment; Step 2.4: Equivalent static load at the m-th time step obtained in Step 2.3 The equivalent static load data for n time steps are obtained, and the n equivalent static load vectors are weighted, summed, or averaged to obtain the final equivalent static load for topology optimization. ; Final equivalent static load for: ； In the formula: This is the final equivalent static load; Let be the equivalent static load vector at time step m; K is the overall structural stiffness matrix. Let be the nodal displacement vector of the structure at time step m; n is the number of time steps used to construct the final equivalent static load.

4. The method for topology optimization design for failure safety of fiber composites under dynamic load as described in claim 1, characterized in that, Step 3 includes: Step 3.1: Establish a failure safety topology optimization model for fiber-reinforced composite materials with material density variable, angle sub-interval selection variable and fiber angle variable as design variables, and use the maximum value of the compliance under each damage condition as the optimization objective; Step 3.2: Based on the finite element mesh, the candidate local damage blocks are laid out and translated according to the preset mesh size to generate multiple local failure regions. The finite element elements inside and outside each region are assigned the minimum weakening coefficient and 1 respectively to form the corresponding weakening coefficient field, so as to construct the local failure scenario under multiple damage conditions required in Step 3.

1. Step 3.3: Divide the original design interval of the fiber direction into several sub-intervals according to parameter n, calculate the directional elastic matrix corresponding to the central angle of each sub-interval, and introduce discrete-continuous weight variables and small deflection angles within the sub-intervals. Weight and rotate the directional elastic matrix to obtain the equivalent directional stiffness matrix that is continuously related to the fiber angle design variable. Substitute the equivalent directional stiffness matrix into the overall stiffness matrix expression to form the stiffness matrix used in Step 3.

1. Step 3.4: Based on the failure-safe topology optimization model established in Step 3.1 with the maximum flexibility of each local failure condition as the objective, and according to the structural flexibility under each local failure scenario obtained in Step 3.2, the maximum flexibility objective is transformed into a continuously differentiable objective function using the normalized KS aggregation function. Step 3.5: Based on the KS objective function constructed in Step 3.4 and the explicit derivative relationship of compliance with design variables, calculate the sensitivity of the objective function to material density variables, angle sub-interval selection variables, and fiber angle variables, and use the moving asymptote method to iteratively update the three types of design variables based on the sensitivity until the inner loop converges.

5. The method for topology optimization design for failure safety of fiber composites under dynamic load according to claim 4, characterized in that, In step 3.1, the failure-safe topology optimization model for fiber-reinforced composite materials is constructed under the following constraints: ； In the formula: The set of variables that the solver needs to solve; The density of the material; Choose variables for the angle sub-intervals, n sub-intervals; For fiber angle variables; To minimize the objective function f in all m damage scenarios and improve the stiffness under the weakest condition; For the first Flexibility under various damage scenarios; m represents the number of damage scenarios; To find the worst softness in all damage scenarios; In the first Damage scenario, displacement Mechanical equilibrium must be satisfied; This is the equivalent static load; Here is the stiffness matrix for the l-th damage scenario; The displacement vector is obtained under the l-th damage scenario; Material usage constraints; The sum of the densities of all elements; N is the total number of elements; f is the maximum allowable material volume fraction; ; The variable chosen for the angle interval is also in [0,1]. To limit the fiber laying angle Scope; It is the fiber orientation angle; To limit the range of angle change.

6. The method for topology optimization design for failure safety of fiber composites under dynamic load according to claim 4, characterized in that, In step 3.2, the weakening coefficient of the e-th finite element is: ； In the formula: The design domain is given in step 1; For the first There are local failure regions; e is the e-th finite element element within the design domain; For this unit in a partial failure scenario The weakening coefficient below; This is the preset minimum weakening coefficient.

7. The method for topology optimization design of fiber composites under dynamic load failure as described in claim 4, characterized in that, In step 3.4, the maximum compliance objective is transformed into a continuously differentiable objective function using the normalized KS aggregation function: ； In the formula: The objective function value after aggregation by the KS function; It is the natural logarithm function; m represents the total number of local failure scenarios described in step 3.2; This represents the sequence number of the partial failure scenario. =1-m; For the first The structural compliance under equivalent static load in a local failure scenario; e is the base of the natural logarithm, representing the exponential function. ; For regularization parameters; This is the normalization coefficient.

8. The method for topology optimization design for failure safety of fiber composites under dynamic load according to claim 4, characterized in that, In step 3.5, the sensitivity of the objective function to the material density variable is: ； The sensitivity of the angle sub-interval selection variable is: ； The sensitivity of the fiber angle variable is: ； In the formula: Design variables for the material density of the e-th element for the objective function f. Sensitivity; is the objective function value after aggregation by the KS function; m is the total number of local failure conditions; This is the sequence number of the partial failure condition; For the first The structural flexibility under equivalent static load in a partial failure condition; For softness For material density variables The partial derivatives; Design variables for the volume distribution of the e-th unit. Right now ; This refers to the regularization parameter in the KS function; For the first The weighting factor corresponding to each local failure condition; e is the base of the natural logarithm, representing the exponential function. ; The sign for summing over all partial failure conditions; Choose variables for the objective function f for the i-th angle sub-interval in the e-th unit. Sensitivity; Design variables associated with the i-th angular sub-interval of the e-th unit; For softness Choose variables for angle sub-intervals The partial derivatives; The objective function f is the fiber angle variable in the e-th element. Sensitivity; Design variables for the fiber angle of the e-th finite element; In the first Compliance under partial failure conditions Design variables for unit fiber angle The partial derivatives of .

9. The method for topology optimization design of fiber composites under dynamic load failure safety according to claim 2, characterized in that, In step 1.3, the stress condition of the semi-three-point bending beam under dynamic load is set as follows: a fixed constraint is applied to the lower right corner of the design domain, a time-varying dynamic load is applied to the upper left corner of the design domain, and a horizontal degree of freedom constraint is applied to the left end face of the design domain; the areas of a preset unit width on both sides of the design domain are set as non-failure regions. : x and y are spatial coordinate variables in the structural coordinate system where the design domain is located, where x is the coordinate in the horizontal direction along the length of the semi-three-point curved beam, and y is the coordinate in the vertical direction along the height of the semi-three-point curved beam; L is the length of the design domain along the beam length; b is the preset unit width of the non-failure regions set on the left and right sides of the design domain along the beam length.

10. A topology optimization design system for the failure safety of fiber composites under dynamic loads, characterized in that, The system includes: The design parameter definition module is used to define the design parameters of the design domain of the fiber-reinforced composite structure based on the stress conditions and geometric characteristics of the structure to be optimized. The dynamic analysis and equivalent static load generation module is used to perform dynamic mechanical analysis on fiber-reinforced composite structures based on the design domain constructed by the design parameter definition module, and obtain the equivalent static load for topology optimization. The failure-safe topology optimization iteration module is used to replace the dynamic load of the original structure with the equivalent static load obtained by the dynamic analysis and equivalent static load generation module. The model with the applied equivalent static load is used as the analysis object. The material density variable, angle sub-interval selection variable and fiber angle variable in the above model are used as design variables to construct multiple local failure conditions. Based on the equivalent directional stiffness matrix and normalized KS objective function that are continuously related to the fiber angle, the failure-safe topology optimization iteration is performed on the design variables in combination with the moving asymptote method until the inner loop converges. The outer loop convergence judgment and result output module is used to determine whether the outer loop has converged based on whether the relative change in the flexibility of the optimized structure under the current equivalent static load obtained by the failure-safe topology optimization iteration module is less than 0.1%. When converged, the failure-safe topology optimization result of the fiber-reinforced composite material is output. Otherwise, the dynamic analysis is re-performed on the updated structure and a new equivalent static load is constructed before returning to the failure-safe topology optimization iteration module to continue optimization.

Images