A curve fiber wing intelligent optimization algorithm fusing physical information bending agent

Abstract

This invention discloses an intelligent optimization algorithm for curved fiber wing that integrates physical information buckling proxy, relating to the field of aircraft design technology. The method includes: equivalently planarizing the spatially twisted quadrilateral; deriving the Jacobian matrix through a double-triangle affine mapping to construct a variable stiffness total potential energy functional; embedding the mechanical control equations using a physical information neural network (PINN) to quickly solve for the buckling critical load factor and modes, replacing traditional finite element calculations; parameterizing the fiber path and thickness to construct design variables; and employing a multi-stage adaptive differential evolution algorithm for global optimization, integrating process, static / dynamic aeroelasticity, and strength constraints. This invention achieves a closed-loop modeling-evaluation-optimization process, significantly reducing computational costs while maintaining mechanical consistency, solving the convergence problem of high-dimensional variable coupling, and greatly improving the automation and engineering efficiency of curved fiber wing design.

Description

Technical Field

[0001] This invention relates to the field of aircraft design technology, and in particular to a smart optimization algorithm for curved fiber wings that integrates physical information buckling agents. Background Technology

[0002] In the modern aerospace field, wing skin structures are one of the most important load-bearing components of aircraft, and their structural weight, stiffness distribution, and aeroelastic properties directly affect flight performance and structural safety. With the widespread application of advanced composite materials in wing skin structures, by rationally designing the layup angle, thickness distribution, and fiber path parameters, it is possible to significantly reduce the structural weight while meeting engineering requirements such as strength, stiffness, and stability. This has become an important technical approach to improving the overall performance of aircraft.

[0003] In the existing technology, the high computational cost and low efficiency of the optimization design of curved fiber composite wing skin structure are usually improved in two ways: first, by introducing a surrogate model to reduce the use of high-cost numerical simulations; and second, by introducing an optimization algorithm to achieve automatic search of design variables.

[0004] Regarding surrogate models, existing methods mainly employ traditional approximation models such as multinomial response surface models, Kriging models, and radial basis function models, or data-driven neural network models to establish a nonlinear mapping relationship between design variables and structural responses. However, these methods generally suffer from problems such as strong dependence on high-quality sample data, difficulty in ensuring physical consistency, and insufficient generalization ability in high-dimensional complex problems, making it difficult to meet the prediction accuracy requirements of high-dimensional design problems for curved fiber optic wing skins.

[0005] In terms of optimization algorithms, current engineering practices mainly employ sensitivity optimization algorithms for local optimization or swarm intelligence optimization algorithms such as genetic algorithms for global search. The former is typically suitable for low-dimensional continuous variable problems, but it is prone to getting trapped in local optima when facing highly nonlinear and highly coupled design variables such as curved fiber layups. While the latter has global search capabilities, its optimization process relies on a large number of individual evaluations, and each evaluation usually requires calling high-cost finite element and aeroelastic analysis processes, resulting in low overall computational efficiency and making it difficult to support effective search in large-scale, high-dimensional design spaces. Summary of the Invention

[0006] The purpose of this invention is to provide a smart optimization algorithm for curved fiber wing that integrates physical information buckling agent, aiming to solve or improve at least one of the above-mentioned technical problems.

[0007] To achieve the above objectives, the present invention provides the following solution: A smart optimization algorithm for curved fiber optic wings that integrates physical information buckling proxies includes: The quadrilateral region with spatial distortion is equivalently planarized and projected to obtain the two-dimensional physical domain Ω; Define reference domain And for the physical domain Ω and the reference domain Perform a bi-triangular affine mapping to derive the Jacobian matrix; Based on the physical domain Ω and the Jacobian matrix, combined with the curve fiber path parameters, the equivalent stiffness matrix that varies with the space is derived, and the energy functional of the governing equation is mapped from the physical domain to the triangular reference domain to construct the total potential energy functional of the laminate. Based on the total potential energy functional of the laminate, the buckling critical load factor is solved using the physical information neural network PINN to generate the first-order buckling factor and buckling mode. Establish a mapping mechanism from engineering manufacturing constraints to mathematical optimization variables, and construct a global design variable vector; The process can construct constraints, static aeroelastic constraints, dynamic aeroelastic constraints, and strength constraints, and uses an improved adaptive differential evolution algorithm to perform global optimization on the global design variable vector to generate the optimal design variable vector.

[0008] Furthermore, the spatially distorted quadrilateral region is equivalently planarized, including: Calculate the geometric center O of the four vertices of the quadrilateral region Q. Based on the coordinate deviation of each vertex relative to the geometric center O, construct the covariance matrix C, expressed as: In the formula, Let be the global spatial coordinate vector of the i-th vertex of region Q, arranged in counterclockwise order; Let O be the coordinate vector of the geometric center. Eigenvalue decomposition is performed on the covariance matrix C to find the minimum eigenvalue and the corresponding unit eigenvector n. This eigenvalue is then used as the normal to the best-fit plane of region Q to determine the equivalent planar region Q'. The expression is as follows: In the formula, n is the normal to the best-fit plane; v is any eigenvector of the covariance matrix C.

[0009] Furthermore, the projection yields a two-dimensional physical domain Ω, which includes: The direction of the first side of the quadrilateral is selected as the first basis vector. And normalize it, the expression is: In the formula, Let it be the first basis vector of the local coordinate system; Let be the Euclidean length of the first edge; The second basis vector is generated using the cross product operation. The expression is: First basis vectors Second basis vector By combining the normal to the best-fit plane n, a local right-handed orthogonal basis is generated. ; Construct a local rotation matrix based on the local right-handed orthogonal basis, and translate the vertex coordinates in the global coordinate system to the equivalent geometric center O' to obtain the vertex coordinates of the two-dimensional physical domain Ω. The expression is as follows: In the formula, Let be the two-dimensional local coordinates of the i-th vertex in the physical domain Ω; Let Q' be the coordinates of the i-th vertex of the equivalent planar region Q'. The coordinates of the equivalent center point O'; It is a local rotation matrix.

[0010] Furthermore, define the reference domain. And for the physical domain Ω and the reference domain Perform a bi-triangular affine mapping to derive the Jacobian matrix, including: Define the reference field of the standard unit square Along the diagonal, for the physical domain Ω and the reference domain Perform synchronous partitioning to generate two triangular physical domains and two triangular reference domains; In the reference domain Introducing piecewise local parameters A linear transformation is performed on the triangular reference domain to map the irregular triangle to the standard parameter domain. The expression is as follows: , In the formula, The lower triangular reference domain; The upper triangular reference domain; For global coordinate variables of the reference domain; Construct an affine mapping function from the triangular reference domain to the triangular physical domain, and perform linear interpolation using the vertex coordinates of the physical domain. The expression is as follows: In the formula, Let be the vertex coordinate vector of the physical domain Ω; Let be the coordinate vector of any point within the physical domain Ω; Within the triangular reference domain, the derivation of local parameters is performed respectively. The Jacobian matrix to the physical domain Ω is expressed as: In the formula, Let be the Jacobian matrix of the i-th triangular reference field.

[0011] Furthermore, based on the physical domain Ω and the Jacobian matrix, combined with the curved fiber path parameters, the equivalent stiffness matrix that varies with space is derived. The energy functional of the governing equations is then mapped from the physical domain to the triangular reference domain, constructing the total potential energy functional of the laminate, including: Using the local coordinate system of the physical domain Ω as a basis, the translation method is used to assume that the curve fiber follows... Extending in direction, define the fiber angle function for the k-th layer layup; Based on material engineering constants, a single-layer stiffness matrix of the material principal axis is constructed, and the direction cosine and sine of the local fiber angle are calculated according to the fiber angle function. The single-layer stiffness matrix is ​​transformed to the local coordinate system of the physical domain Ω using the tensor transformation rule to obtain the off-axis stiffness matrix. Integrating the off-axis stiffness matrix along the thickness direction of the laminate yields the in-plane stiffness matrix, the coupling stiffness matrix, and the bending stiffness matrix. Based on Kirchhoff's thin plate theory, the mid-surface displacement field is defined, the relationship between the displacement field and the strain vector and curvature vector is established, and the strain energy density and geometric potential energy density per unit area are constructed by combining the in-plane stiffness matrix, the coupling stiffness matrix and the bending stiffness matrix. Based on the Jacobian matrix and the corresponding inverse matrix, the differential operator on the physical domain Ω is transformed into an operator on the reference domain; Based on the operator transformation and the strain energy density and geometric potential energy density per unit area in the physical domain Ω, the total potential energy functional in the physical domain Ω is mapped to the triangular reference domain to obtain the total potential energy functional of the laminate.

[0012] Furthermore, based on the total potential energy functional of the laminate, the physical information neural network PINN is used to solve for the buckling critical load factor, generating the first-order buckling factor and buckling modes, including: Based on the mapping relationship between the physical domain Ω and the triangular reference domain, the integration domain is transformed to the triangular reference domain, and the approximate integral at the sampling points is obtained using the rational calculation method. The expression is as follows: In the formula, N is the total number of sampling points; Let be the local coordinates of the i-th sampling point within the reference domain; The absolute value of the row and column of the Jacobian matrix at the i-th sampling point is used to convert the density values ​​in the reference domain into energy contributions in the physical domain using weighted averages. Let be the strain energy density at the i-th sampling point; Let be the geometric potential energy density at the i-th sampling point; Combining the Rayleigh quotient objective and normalization constraints, the total loss function of the neural network is constructed as follows: In the formula, This is the total loss function of the physical information neural network; This is the set of weights and bias parameters for the neural network. The strain energy density is calculated based on the output displacement field and automatic differential derivative of the neural network. The geometric potential energy density is calculated based on the first derivative of the displacement field output by the neural network. β These are the component weights of the loss function; This is the normalized loss term; The normalized loss term for the buckling modal displacement field is expressed as follows: In the formula, This is the normalized loss term; Let be the lateral deflection value predicted by the neural network at the i-th sampling point; These are the trainable parameters of the neural network; Area weighting factor; By minimizing the total loss function through an optimization algorithm, and after training convergence, the first-order buckling factor and buckling mode are output.

[0013] Furthermore, a mapping mechanism is established from engineering manufacturing constraints to mathematical optimization variables, and a global design variable vector is constructed, including: For the upper and lower skin of the wing structure, design domains are constructed separately. These design domains are then divided into non-overlapping discrete sub-regions along the spanwise and chordwise directions of the wing, expressed as: , In the formula, and These are the design domains for the upper skin and the lower skin, respectively. and These are the r-th sub-regions of the upper and lower skin, respectively; and These represent the number of sub-regions in the upper and lower skin, respectively. For any subregion, the thickness design variable is defined as the total thickness of the subregion. Assuming that the thicknesses of different plies are equal, the expression for the thickness of the k-th ply is: In the formula, The total thickness of the r-th sub-region is a design variable. Let r be the total number of layers in the r-th subdomain; The thickness of a single layer in the k-th layer within the r-th sub-region; The fiber path for each layup is set by The control points are defined, and the expression for the fiber path angle sequence of the k-th layer of the upper and lower skin is as follows: , In the formula, This is the fiber path angle sequence of the k-th layer of the upper skin; The path angle sequence of the k-th layer of the lower skin; The fiber angle at any position y is obtained by linear interpolation of the angles of adjacent control points, and the expression is: In the formula, The number of control points for a single fiber path; Design variables for the fiber angle of the k-th layer at the k-th control point; It is an interpolation function; Based on the thickness variables of all sub-regions and the fiber angle control point variables of each layup, a unified global design variable vector is constructed, expressed as: In the formula, Design a global variable vector; Design a variable vector for thickness; Design a variable vector for fiber angles; Thickness design variable vector and fiber angle design variable vector The expression is: In the formula, For the r-th sub-region, the k-th layer layup, and the fiber angle at the m-th control point; This represents the total number of skinned sub-regions; Number of layers for the skin; This represents the number of skin control points.

[0014] Furthermore, the process can manufacture constraints, static aeroelastic constraints, dynamic aeroelastic constraints, and strength constraints, which respectively include: Manufacturability constraints include: Thickness process constraints, including: For spanwise adjacent sub-regions, the ratio of layer drop thickness to spanwise length must satisfy the spanwise thickness change rate constraint, expressed as: In the formula, Let be the total ply thickness of the i-th and j-th sub-regions, where i is the spanwise partition number and j is the chordwise partition number; The spanwise distance between adjacent spanwise intervals; This represents the maximum allowable rate of change of the spanwise angle in the process. For spanwise adjacent subregions, the ratio of layer drop thickness to spanwise length satisfies the chordal thickness variation rate constraint, expressed as: In the formula, The chordal distance between adjacent chordal intervals; This represents the maximum allowable rate of change of the chordal angle in the process. Angular process constraints, including: For adjacent layers within the same subregion, the fiber angles satisfy the interlayer angle continuity constraint, expressed as: In the formula, For the i-th spanwise region, the j-th chordwise region, the k-th layup, and the fiber angle of the m-th grid column; For adjacent grids within the same subregion, the fiber angle satisfies the spanwise angle change rate constraint, expressed as: In the formula, The distance to the center of the neighboring column grid; This represents the maximum allowable rate of angular change in the process. Static aeroelastic constraints include: The wingtip deflection constraint is expressed as follows: In the formula, Let x be the deflection penalty value for individual x; This is the upper limit of the deflection penalty value; This is the normalized deflection over-limit ratio; The nonlinear exponent is used to penalize deflection. This represents the maximum threshold for wingtip deflection. Let x be the wingtip deflection of individual x; The wingtip twist angle constraint is expressed as follows: In the formula, The torsion angle penalty value for individual x; This represents the maximum penalty value for the torsion angle. The normalized torsion angle over-limit ratio; The nonlinear exponent is penalized for the torsion angle; This represents the maximum threshold value for the wingtip twist angle. Let x be the wingtip twist angle of individual x; The derivative of the control surface and the stiffness ratio constraint are expressed as follows: In the formula, The penalty value is the ratio of rudder effect to spring stiffness. For individual x, the rudder-elasticity ratio is given. For in the interval The rudder-elasticity ratio after linear normalization. The nonlinear exponent of the rudder-efficiency stiffness ratio penalty value; The upper limit of the penalty value for the rudder effect-to-elasticity ratio; This is the maximum penalty threshold; No penalty threshold; Aeroelastic constraints include: The flutter velocity constraint is expressed as follows: In the formula, This is the flutter speed penalty value; Let x be the flutter velocity of individual x; This is the maximum penalty threshold; No penalty threshold; This represents the upper limit of the flutter speed penalty value; The nonlinear exponent for the flutter velocity penalty value; Strength constraints include: The expression for composite strain constraint is: In the formula, Let x be the strain penalty value for individual x; This represents the upper limit of the strain penalty value; and For control parameters; The total number of units that violate the constraints; To assess the risk level; For design domain; It is a non-linear exponent; The element strain exceeds the limit; The number of layers in a unit; Let be the allowable strain value for the j-th strain; For the first The k-th layer of a unit, and the j-th strain component; Local buckling constraint, expressed as: In the formula, Let x be the buckling penalty value for individual x; This represents the maximum value for the buckling penalty. To assess the degree of risk of buckling; and It is a nonlinear parameter; The total number of regions where buckling occurs; Let be the first-order buckling factor of the j-th region; This is the safety threshold for the buckling factor.

[0015] Furthermore, an improved adaptive differential evolution algorithm is used to globally optimize the global design variable vector, generating the optimal design variable vector, including: Minimizing individual fitness as the optimization objective is expressed as: In the formula, For individual fitness; This is the total penalty value; For quality adaptability; The expression for the total penalty value is: In the formula, This is the deflection penalty value; This is the torsion angle penalty value; The penalty value is the ratio of rudder effect to spring stiffness. This is the flutter speed penalty value; This is the strain penalty value; An improved linear population size reduction differential evolution algorithm framework is adopted, which introduces a linear population decay strategy and a three-stage differentiation evolution mechanism. By dynamically adjusting the population size and control parameters, the update strategy is iteratively executed, and the optimal design variable vector is obtained after convergence.

[0016] Furthermore, an improved linear population size reduction differential evolution algorithm framework is adopted, introducing a linear population decay strategy and a three-stage differentiation evolution mechanism. By dynamically adjusting the population size and control parameters, an update strategy is iteratively executed, including: Population linear decay strategies include: The expression for the population in generation g is: In the formula, This is the g-th generation population; This is the i-th individual in the g-th generation population; The total number of individuals in the g-th generation population; This refers to the number of individuals in the first generation population. The number of individuals in the Gth generation population; The maximum number of generations is preset. Using the population mutation strategy under the L-SHADE framework, mutated individual vectors are generated, expressed as follows: In the formula, Let be the mutated individual vector of the i-th individual in the g-th generation population; This is the i-th individual in the g-th generation population; Individuals that rank in the top p% of the g-th generation population; and Two individuals are randomly selected from the population in the g-th generation. This is the difference scaling factor; Vector of mutated individuals With the target vector Perform binomial crossover to generate the test vector, expressed as: In the formula, Let j be the j-th component of the test vector; Let be the crossover probability of the i-th individual; This is a random number corresponding to the weight j of the i-th individual in the g-th generation population; Calculate the test vector fitness value If greater than or equal to the target vector The fitness value, then the trial vector Enter the next generation of the population, and record the differential scaling factor at the same time. and crossover probability Store in a temporary successful set; otherwise, retain the original individual. ; After each generation, the historical mean of the success scaling factor in the historical memory bank is updated using the weighted Lehmer mean and the arithmetic mean. Historical average of the probability of successful crossover .

[0017] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects: This invention discloses an intelligent optimization algorithm for curved fiber optic wings that integrates physical information buckling proxies, with the following beneficial effects: 1. This invention employs a buckling assessment method based on a physical information neural network. The mechanical governing equations and energy functionals of composite laminates are directly introduced into the neural network training process, enabling rapid prediction of the first-order buckling factor. Compared to traditional buckling analysis methods that rely on finite element eigenvalue solutions, this method eliminates the need for establishing a refined finite element model and a global stiffness matrix. While maintaining mechanical consistency, it significantly reduces computational costs and provides a feasible path for the frequent invocation of buckling constraints in high-dimensional optimization processes.

[0018] 2. This invention preserves the explicit mapping relationship between layup angle and thickness on the stiffness characteristics of laminates during the modeling process. Curved fiber paths directly participate in structural performance evaluation through control point parameterization, ensuring that adjustments to design variables during optimization are accurately reflected in the structural stiffness distribution, stability characteristics, and aeroelastic response. Compared to traditional optimization methods that only use equivalent parameters or simplified models, this invention can more realistically characterize the engineering effects of curved fiber designs.

[0019] 3. To address the challenges of high-dimensional design variables, strong variable coupling, and a narrow feasible solution space in curved fiber optic wing skin design, this invention introduces a multi-stage population optimization strategy. This strategy adaptively adjusts the search strategy according to different stages of the optimization process, achieving a smooth transition from global exploration to local fine-grained search. Compared to traditional single genetic algorithms or gradient-based optimization algorithms, this method exhibits better search stability and convergence efficiency in high-dimensional design spaces.

[0020] 4. In the optimization process, this invention comprehensively introduces process constraints, stability constraints, strength constraints, aeroelastic response constraints, and aeroelastic stability constraints. It can directly use commonly used evaluation indicators in engineering design (such as buckling factor, strain limit, rudder effect stiffness ratio, wingtip deflection, flutter velocity, etc.) as the judgment basis, so that the method can be directly connected to the real engineering design process and has good engineering interpretability and applicability.

[0021] 5. This invention integrates the curved fiber parametric modeling module, the PINN mechanical evaluation module, and the multi-stage group optimization module to form a closed-loop process of "design variable update—performance evaluation—constraint determination—optimization iteration". This process can automatically complete large-scale scheme search and performance screening without repeated manual intervention, effectively improving the automation level and design efficiency of curved fiber wing structure optimization design.

[0022] 6. This invention introduces a rapid buckling assessment method based on PINN and uses it as the core performance evaluation module in the optimization process. This significantly reduces the number of calls to traditional finite element buckling solutions and costly simulation analyses during the evaluation of a large number of candidate design schemes. Compared to optimization processes entirely driven by finite elements, this method greatly reduces the overall computational cost while maintaining engineering accuracy, making it more suitable for engineering optimization design scenarios involving complex structures. Attached Figure Description

[0023] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0024] Figure 1 This is a schematic flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the wing optimization process in this embodiment. Detailed Implementation

[0025] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0026] The purpose of this invention is to provide a smart optimization algorithm for curved fiber wing that integrates physical information buckling agent, aiming to solve or improve at least one of the above-mentioned technical problems.

[0027] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0028] like Figure 1 As shown, this invention provides a smart optimization algorithm for curved fiber optic wings that integrates physical information buckling proxies, including: Step 1: The spatially distorted quadrilateral region is equivalently planarized and projected to obtain a two-dimensional physical domain Ω, including: Step 11: Solve for the best-fit plane based on the covariance matrix of the quadrilateral region Q to generate the equivalent planar region Q', including: Calculate the geometric center O of the four vertices of the quadrilateral region Q. Based on the coordinate deviation of each vertex relative to the geometric center O, construct the covariance matrix C, expressed as: In the formula, Let be the global spatial coordinate vector of the i-th vertex of region Q, arranged in counterclockwise order; Let O be the coordinate vector of the geometric center. Eigenvalue decomposition is performed on the covariance matrix C to find the minimum eigenvalue and the corresponding unit eigenvector n. This eigenvalue is then used as the normal to the best-fit plane of region Q to determine the equivalent planar region Q'. The expression is as follows: In the formula, n is the normal to the best-fit plane; v is any eigenvector of the covariance matrix C.

[0029] Step 12, using the equivalent geometric center O' of the equivalent planar region Q' as the origin, construct a local right-handed orthogonal basis, including: The direction of the first side of the quadrilateral is selected as the first basis vector. And normalize it, the expression is: In the formula, Let it be the first basis vector of the local coordinate system; Let be the Euclidean length of the first edge; The second basis vector is generated using the cross product operation. The expression is: First basis vectors Second basis vector By combining the normal to the best-fit plane n, a local right-handed orthogonal basis is generated. .

[0030] Step 13: Construct a local rotation matrix based on the local right-handed orthogonal basis, and translate the vertex coordinates in the global coordinate system to the equivalent geometric center O' to obtain the vertex coordinates of the two-dimensional physical domain Ω. The expression is: In the formula, Let be the two-dimensional local coordinates of the i-th vertex in the physical domain Ω; Let Q' be the coordinates of the i-th vertex of the equivalent planar region Q'. The coordinates of the equivalent center point O'; It is a local rotation matrix.

[0031] Step 2, Define the reference domain And for the physical domain Ω and the reference domain Perform a bi-triangular affine mapping to derive the Jacobian matrix, including: Step 21, Define the reference domain of the standard unit square. The coordinates of the reference domain vertices, in counter-clockwise order, are as follows: , , , ; Along the diagonal, for the physical domain Ω and the reference domain Perform synchronous partitioning to generate two triangular physical domains and two triangular reference domains; Step 22, in the reference domain Introducing piecewise local parameters A linear transformation is performed on the triangular reference domain to map the irregular triangle to the standard parameter domain. The expression is as follows: , In the formula, The lower triangular reference domain; The upper triangular reference domain; For global coordinate variables of the reference domain; Step 23: Construct an affine mapping function from the triangular reference domain to the triangular physical domain, using linear interpolation based on the vertex coordinates of the physical domain. The expression is: In the formula, Let be the vertex coordinate vector of the physical domain Ω; Let be the coordinate vector of any point within the physical domain Ω; Within the triangular reference domain, the derivation of local parameters is performed respectively. The Jacobian matrix to the physical domain Ω is expressed as: In the formula, Let be the Jacobian matrix of the i-th triangular reference field.

[0032] As mentioned above, the elements of the Jacobian matrix are independent of local parameters and depend only on the coordinates of the vertices in the physical domain; therefore, its determinant is a constant. The Jacobian matrix is ​​composed of the coordinate differences of the vertices in the physical domain and represents the scaling and shearing deformation relationship from the reference parameter space to the actual physical space.

[0033] Step 3: Based on the physical domain Ω and the Jacobian matrix, and combined with the curved fiber path parameters, derive the equivalent stiffness matrix that varies with space. Map the energy functional of the governing equations from the physical domain to the triangular reference domain to construct the total potential energy functional of the laminate, including: Step 31: Using the local coordinate system of the physical domain Ω as a basis, the translation method is used to assume that the curve fiber follows... Extending in direction, the fiber angle function of the k-th layer is defined as follows: In the formula, The fiber angle function for the k-th layer layup; These are local coordinate variables; The starting fiber angle for the nth interval; The termination fiber angle for the nth interval; The length of the nth interval; n is the interval index; Let j be the length of the j-th interval; Step 32: Construct the single-layer stiffness matrix of the material principal axis based on material engineering constants, and calculate the direction cosine and sine of the local fiber angles according to the fiber angle function. Use the tensor transformation rule to transform the single-layer stiffness matrix to the local coordinate system of the physical domain Ω to obtain the off-axis stiffness matrix. ,include: The expression for the single-layer stiffness matrix is: In the formula, The stiffness matrix is ​​for a single layer. The tensile modulus in the first direction; The tensile modulus in two directions; Shear modulus in 12 directions; Poisson's ratio in 12 directions; The Poisson's ratio in direction 21; The expression for the off-axis stiffness matrix is: In the formula, Let be the rotational stiffness coefficient of the k-th layer; The direction cosine of the local fiber angle; The sine of the local fiber angle; Let be the stiffness coefficient of the k-th layer in the principal axis coordinate system of the material; These represent the normal stress and shear direction in the corresponding local coordinate system, respectively. In this embodiment, 1 represents... ,2 represents 6 represents ,For example Indicates in When the direction produces unit strain, in Stress is generated in the direction.

[0034] Step 33: Integrate the off-axis stiffness matrix along the thickness direction of the laminate to obtain the in-plane stiffness matrix, coupling stiffness matrix, and bending stiffness matrix, including: Define the total thickness of the laminate and the coordinates of each layer interface as follows: , , , In the formula, This represents the total thickness of the laminate. This represents the total number of layers in the laminate. Let the thickness be the k-th layer; The middle surface of the laminate coordinate; The upper surface of the k-th layer coordinate; The lower surface of the k-th layer coordinate; The mid-surface of the k-th layer coordinate; The expressions for the in-plane stiffness matrix, coupling stiffness matrix, and bending stiffness matrix are as follows: In the formula, Here is the in-plane stiffness matrix; Here is the coupling stiffness matrix; Here is the bending stiffness matrix; Construct the variable stiffness constitutive relation, expressed as: , , , In the formula, The resultant external force matrix; It is the external torque vector; Here is the in-plane stiffness matrix; Here is the coupling stiffness matrix; Here is the out-of-plane stiffness matrix; The mid-surface strain vector; It is a curvature vector; , and They are respectively , And the in-plane and out-of-plane force components in the shear direction; , and External torque , as well as Towards components; , and The middle surface of the laminate is respectively , and shear strain; , and The middle surface of the laminate is respectively , Curvature and twist rate.

[0035] Step 34: Based on Kirchhoff's thin plate theory, define the mid-surface displacement field, establish the relationship between the displacement field and the strain vector and curvature vector, and construct the strain energy density and geometric potential energy density per unit area, including: The expression for the surface displacement field in a laminated plate is: In the formula, The displacement of the mid-surface of the laminate along axis 1 of the material coordinate system; The displacement of the mid-surface of the laminate along the 2-axis of the material coordinate system; This refers to the out-of-plane displacement of the laminate. Mid-surface strain vector and curvature vector The expression is: , Calculate the strain of the laminate The relationship between thickness and other properties is as follows: In the formula, In the thickness direction; The strain energy density and geometric potential energy density per unit area within the physical domain Ω are expressed as follows: In the formula, Strain energy density; The geometric potential energy density; , and The preload is the in-plane force obtained from global aeroelastic analysis.

[0036] Step 35, based on the Jacobian matrix and its corresponding inverse matrix, transform the differential operator in the physical domain Ω into an operator in the reference domain, including: Let the inverse of the Jacobian matrix be Then for any scalar field The first derivative transformation is expressed as: , In the formula, Inverse matrix The elements are constants, representing the linear coefficients of the coordinate transformation; The second derivative transformation is expressed as: Step 36: Based on the operator transformation and the strain energy density and geometric potential energy density per unit area within the physical domain Ω, the total potential energy functional over the physical domain Ω is mapped to the triangular reference domain to obtain the total potential energy functional of the laminate, expressed as: In the formula, The total potential energy functional of the laminate; Let i be the triangular reference domain; The strain energy density function mapped to the reference domain; This is the load proportionality factor; Let be the geometric potential energy density function mapped to the reference domain; For coordinate variables within the reference domain; A constant area scaling factor; Let be the Jacobian matrix from the reference domain to the physical domain for the i-th triangular reference domain.

[0037] Step 4: Based on the total potential energy functional of the laminate, the physical information neural network PINN is used to solve for the buckling critical load factor, generating the first-order buckling factor and buckling modes, including: According to the minimum potential energy primitive, the first variation of the total potential energy functional of the structure in equilibrium is zero; when the structure is in the critical buckling state, its stability is determined by the second variation, that is, the second variation of the total potential energy functional is zero.

[0038] Step 41, the energy variational condition expression under the critical state is: , In the formula, For the second variation of strain energy, denoted as the restoring stiffness of the structure; The second variation of the geometric potential energy represents the softening stiffness caused by the preload; As the loading factor, when Reaching the minimum value At this time, the structure undergoes first-order buckling.

[0039] The above expression shows that the buckling problem is essentially a generalized eigenvalue problem, and its physical solution corresponds to the minimum eigenvalue of the system's total stiffness degradation.

[0040] Step 42, construct the buckling factor objective function based on Rayleigh quotient, including: The expression for the first-order buckling factor is: In the formula, It is a first-order buckling factor; Let be the displacement field function of the mid-surface of the laminated plate, where Lateral deflection plays a dominant role; Strain energy density; The geometric potential energy density; Based on the mapping relationship between the physical domain Ω and the triangular reference domain, the integration domain is transformed to the triangular reference domain, and the approximate integral at the sampling points is obtained using the rational calculation method. The expression is as follows: In the formula, N is the total number of sampling points; Let be the local coordinates of the i-th sampling point within the reference domain; The absolute value of the row and column of the Jacobian matrix at the i-th sampling point is used to convert the density values ​​in the reference domain into energy contributions in the physical domain using weighted averages. Let be the strain energy density at the i-th sampling point; Let be the geometric potential energy density at the i-th sampling point.

[0041] Since the generalized eigenvalue problem has a zero solution (i.e., when the displacement field is all zero, both the numerator and denominator are zero, and the ratio is meaningless), and the eigenvector has arbitrary scaling, directly optimizing the Rayleigh quotient may cause the network to converge to the zero solution or a numerically unstable maximum / minimum.

[0042] Step 43, construct the total loss function containing the Rayleigh quotient objective term and the modal normalization constraint term, including: Combining the Rayleigh quotient objective and normalization constraints, a total loss function for the neural network is constructed to minimize the loss function, thereby simultaneously obtaining the optimal buckling factor and buckling mode. The expression is: In the formula, This is the total loss function of the physical information neural network; This is the set of weights and bias parameters for the neural network. The strain energy density is calculated based on the output displacement field and automatic differential derivative of the neural network. The geometric potential energy density is calculated based on the first derivative of the displacement field output by the neural network. β These are the component weights of the loss function; This is the normalized loss term; Introducing a normalized loss term for the buckling modal displacement field, the root mean square of the forced displacement field is on the order of unit, expressed as: In the formula, For the normalized loss term, when the weighted mean square value of the displacement modes is close to 1, the loss term is close to 0; Let be the lateral deflection value predicted by the neural network at the i-th sampling point; These are the trainable parameters of the neural network; This is an area weighting factor to ensure that normalization is performed in the sense of the physical domain area, rather than the reference domain.

[0043] Step 44: Minimize the total loss function through the optimization algorithm. After training converges, output the first-order buckling factor and buckling mode.

[0044] Step 5: Establish a mapping mechanism from engineering manufacturing constraints to mathematical optimization variables, and construct a global design variable vector, including: Step 51: For the upper and lower skins of the wing structure, design domains are constructed separately. To accommodate the local characteristics of variable stiffness design, the design domains are divided into non-overlapping discrete sub-regions along the spanwise and chordwise directions of the wing, expressed as: , In the formula, and These are the design domains for the upper skin and the lower skin, respectively. and These are the r-th sub-regions of the upper and lower skin, respectively; and These represent the number of sub-regions for the upper and lower skin, respectively.

[0045] Step 52: For any sub-region, define the thickness design variable as the total thickness of the sub-region. Assuming that the thicknesses of different plies are equal, the expression for the thickness of the k-th ply is: In the formula, The total thickness of the r-th sub-region is a design variable. Let r be the total number of layers in the r-th subdomain; Let be the thickness of the k-th layer within the r-th sub-region.

[0046] To ensure the manufacturability of curved fiber paths, the Shift Method assumption from Automated Fiber Placement (AFP) technology is adopted: that is, within the same sub-region, the fiber paths of all layups have similar shapes, only the angle amplitudes differ, and the angles along the direction perpendicular to the reference path (such as span y) remain constant or change according to a specific pattern.

[0047] Step 53, set the fiber path for each layup layer. The control points are defined, and the expression for the fiber path angle sequence of the k-th layer of the upper and lower skin is as follows: , In the formula, This is the fiber path angle sequence of the k-th layer of the upper skin; The path angle sequence of the k-th layer of the lower skin; The fiber angle at any position y is obtained by linear interpolation of the angles of adjacent control points, and the expression is: In the formula, The number of control points for a single fiber path determines the complexity and degrees of freedom of the path; Design variables for the fiber angle of the k-th layer at the k-th control point; This is an interpolation function used to generate a continuous and smooth fiber angle field.

[0048] Step 54: Based on the thickness variables of all sub-regions and the fiber angle control point variables of each layup, construct a unified global design variable vector, expressed as: In the formula, Design a global variable vector; Design a variable vector for thickness; Design a variable vector for fiber angles; Thickness design variable vector and fiber angle design variable vector The expression is: In the formula, For the r-th sub-region, the k-th layer layup, and the fiber angle at the m-th control point; This represents the total number of skinned sub-regions; Number of layers for the skin; This represents the number of skin control points.

[0049] Step 6: Construct process manufacturability constraints, static aeroelastic constraints, dynamic aeroelastic constraints, and strength constraints, and use an improved adaptive differential evolutionary algorithm to globally optimize the global design variable vector, generating the optimal design variable vector, including: Manufacturability constraints include: Thickness process constraints, including: For spanwise adjacent sub-regions, the ratio of layer drop thickness to spanwise length must satisfy the spanwise thickness change rate constraint, expressed as: In the formula, Let be the total ply thickness of the i-th and j-th sub-regions, where i is the spanwise partition number and j is the chordwise partition number; The spanwise distance between adjacent spanwise intervals; This represents the maximum allowable rate of change in spanwise angle during the process.

[0050] For spanwise adjacent subregions, the ratio of layer drop thickness to spanwise length satisfies the chordal thickness variation rate constraint, expressed as: In the formula, The chordal distance between adjacent chordal intervals; This represents the maximum rate of change of the chordal angle allowed by the process.

[0051] Angular process constraints, including: For adjacent layers within the same subregion, the fiber angles satisfy the interlayer angle continuity constraint, expressed as: In the formula, For the i-th spanwise region, the j-th chordwise region, the k-th layup, and the fiber angle of the m-th grid column; For adjacent grids within the same subregion, the fiber angle satisfies the spanwise angle change rate constraint, expressed as: In the formula, The distance to the center of the neighboring column grid; This represents the maximum allowable rate of angular change in the process.

[0052] Static aeroelastic constraints include: The wingtip deflection constraint is expressed as follows: In the formula, Let x be the deflection penalty value for individual x; This is the upper limit of the deflection penalty value; This is the normalized deflection over-limit ratio; The nonlinear exponent is used to penalize deflection. This represents the maximum threshold for wingtip deflection. Let x be the wingtip deflection of individual x; The wingtip twist angle constraint is expressed as follows: In the formula, The torsion angle penalty value for individual x; This represents the maximum penalty value for the torsion angle. The normalized torsion angle over-limit ratio; The nonlinear exponent is penalized for the torsion angle; This represents the maximum threshold value for the wingtip twist angle. Let x be the wingtip twist angle of individual x; The derivative of the control surface and the stiffness ratio constraint are expressed as follows: In the formula, The penalty value is the ratio of rudder effect to spring stiffness. For individual x, the rudder-elasticity ratio is given. For in the interval The rudder-elasticity ratio after linear normalization. The nonlinear exponent of the rudder-efficiency stiffness ratio penalty value; The upper limit of the penalty value for the rudder effect-to-elasticity ratio; This is the maximum penalty threshold; This is the threshold with no penalty.

[0053] Aeroelastic constraints include: The flutter velocity constraint is expressed as follows: In the formula, This is the flutter speed penalty value; Let x be the flutter velocity of individual x; This is the maximum penalty threshold; No penalty threshold; This represents the upper limit of the flutter speed penalty value; The nonlinear exponent is the flutter velocity penalty value.

[0054] Strength constraints include: The expression for composite strain constraint is: In the formula, Let x be the strain penalty value for individual x; This represents the upper limit of the strain penalty value; and For control parameters; The total number of units that violate the constraints; To assess the risk level; For design domain; It is a non-linear exponent; The element strain exceeds the limit; The number of layers in a unit; Let be the allowable strain value for the j-th strain; For the first The k-th layer of a unit, and the j-th strain component.

[0055] Local buckling constraint, expressed as: In the formula, Let x be the buckling penalty value for individual x; This represents the maximum value for the buckling penalty. To assess the degree of risk of buckling; and It is a nonlinear parameter; The total number of regions where buckling occurs; Let be the first-order buckling factor of the j-th region; The buckling factor safety threshold; Minimizing individual fitness as the optimization objective is expressed as: In the formula, For individual fitness; This is the total penalty value; For quality adaptability; The expression for the total penalty value is: In the formula, This is the deflection penalty value; This is the torsion angle penalty value; The penalty value is the ratio of rudder effect to spring stiffness. This is the flutter speed penalty value; This is the strain penalty value.

[0056] An improved linear population size reduction differential evolution algorithm (L-SHADE) framework is adopted, introducing a linear population decay strategy and a three-stage differentiation evolution mechanism. By dynamically adjusting the population size and control parameters, the update strategy is iteratively executed, and the optimal design variable vector is obtained after convergence, including: Population linear decay strategies include: The expression for the population in generation g is: In the formula, This is the g-th generation population; This is the i-th individual in the g-th generation population; The total number of individuals in the g-th generation population; This refers to the number of individuals in the first generation population. The number of individuals in the Gth generation population; This is the preset maximum number of generations.

[0057] Using the population mutation strategy under the L-SHADE framework, mutated individual vectors are generated, expressed as follows: In the formula, Let be the mutated individual vector of the i-th individual in the g-th generation population; This is the i-th individual in the g-th generation population; Individuals that rank in the top p% of the g-th generation population; and Two individuals are randomly selected from the population in the g-th generation. This is the differential scaling factor, used to control the amplitude of the differential perturbation; its value directly determines the search step size. When the size is large, the algorithm tends to explore on a larger scale, which is beneficial for escaping local optima; when When the step size is smaller, the search step size decreases, which is more conducive to performing fine local searches near the current good solution.

[0058] Vector of mutated individuals With the target vector Perform binomial crossover to generate the test vector, expressed as: In the formula, Let j be the j-th component of the test vector; Let be the crossover probability of the i-th individual; is a random number corresponding to the weight j of the i-th individual in the g-th generation population.

[0059] The above steps differ from the traditional differential evolution algorithm's practice of manually fixing F and CR. This paper employs an adaptive parameter adjustment mechanism based on successful history. The F and CR of each individual are not fixed constants but are randomly sampled from a historical memory distribution, which is continuously updated based on the parameters corresponding to "successfully generating candidate individuals superior to their parent's solutions." This mechanism enables the algorithm to automatically favor larger perturbation parameters in the early stages of the search to enhance exploratory behavior, and gradually converge to smaller parameters in the later stages to improve local search accuracy, thus achieving an adaptive balance between exploration and development.

[0060] The three-stage population renewal strategy includes: The current search phase is defined using a periodic phase division rule, expressed as: In the formula, These correspond to the exploration phase, the transition phase, and the development phase, respectively. The above-mentioned stage discrimination method is independent of the total number of algebras and is determined only by the current algebra, so that the algorithm presents a stable "exploration-convergence-refinement" cycle rhythm throughout the optimization process.

[0061] The core idea of ​​the three-stage regulation is reflected in the following: The phased adjustment of F and CR.

[0062] Specifically, during the exploration phase, to enhance global search capabilities, the algorithm employs a larger... The parameters F and CR make the differential direction sources more dispersed and the step size larger, which helps to cover a wider design space and avoid premature convergence; during the development phase, the parameter settings tend to be conservative. The perturbation range of F and CR is significantly reduced, allowing the search to mainly revolve around the current optimal region, thereby improving the accuracy of local optimization. The parameter settings in the transition phase are between the two, playing a role in smoothly transitioning from global exploration to local development.

[0063] Calculate the test vector fitness value If greater than or equal to the target vector The fitness value, then the trial vector Enter the next generation of the population, and record the differential scaling factor at the same time. and crossover probability Store in a temporary successful set; otherwise, retain the original individual. ; After each generation, the historical mean of the success scaling factor in the historical memory bank is updated using the weighted Lehmer mean and the arithmetic mean. Historical average of the probability of successful crossover .

[0064] The aforementioned update mechanism enables the algorithm to "learn" the most effective combination of parameters for the current evolutionary stage (exploration, transition, or development), thereby achieving true adaptive three-stage evolution.

[0065] The beneficial effects of this invention include: 1. It can achieve efficient evaluation of buckling stability indices while ensuring physical consistency. This invention employs a buckling assessment method based on a physical information neural network. The mechanical governing equations and energy functionals of composite laminates are directly incorporated into the neural network training process, enabling rapid prediction of the first-order buckling factor. Compared to traditional buckling analysis methods that rely on finite element eigenvalue solutions, this method eliminates the need for a refined finite element model and a global stiffness matrix. While maintaining mechanical consistency, it significantly reduces computational costs and provides a feasible path for the frequent invocation of buckling constraints in high-dimensional optimization processes.

[0066] 2. It can accurately reflect the influence of curved fiber layup parameters on structural performance. This invention preserves the explicit mapping relationship between layup angle and thickness on the stiffness characteristics of laminates during the modeling process. Curved fiber paths directly participate in structural performance evaluation through control point parameterization, ensuring that adjustments to design variables during optimization are accurately reflected in the structural stiffness distribution, stability characteristics, and aeroelastic response. Compared to traditional optimization methods that only use equivalent parameters or simplified models, this invention can more realistically characterize the engineering effects of curved fiber designs.

[0067] 3. It can effectively adapt to high-dimensional design variable spaces, improving the feasibility of complex optimization problems. To address the challenges of high-dimensional design variables, strong variable coupling, and a narrow feasible solution space in curved fiber optic wing skin design, this invention introduces a multi-stage population optimization strategy. This strategy adaptively adjusts the search strategy according to different stages of the optimization process, achieving a smooth transition from global exploration to local fine-grained search. Compared to traditional single genetic algorithms or gradient-based optimization algorithms, this method exhibits better search stability and convergence efficiency in high-dimensional design spaces.

[0068] 4. Capable of handling multiple types of engineering constraints simultaneously, exhibiting good engineering applicability. This invention comprehensively incorporates process constraints, stability constraints, strength constraints, aeroelastic response constraints, and aeroelastic stability constraints during the optimization process. It can directly use commonly used evaluation indicators in engineering design (such as buckling factor, strain limit, rudder effect stiffness ratio, wingtip deflection, flutter velocity, etc.) as the judgment criteria, enabling the method to directly connect with the real engineering design process and possessing good engineering interpretability and applicability.

[0069] 5. It can form a complete closed-loop process from evaluation to optimization, improving the degree of design automation. This invention integrates a curved fiber parametric modeling module, a PINN mechanical evaluation module, and a multi-stage group optimization module to form a closed-loop process of "design variable update—performance evaluation—constraint determination—optimization iteration." This process can automatically complete large-scale scheme search and performance screening without repeated manual intervention, effectively improving the automation level and design efficiency of curved fiber wing structure optimization design.

[0070] 6. It can significantly reduce reliance on traditional, high-cost finite element simulation. Because this invention introduces a rapid buckling assessment method based on PINN and uses it as the core performance evaluation module in the optimization process, it can significantly reduce the number of calls to traditional finite element buckling solutions and high-cost simulation analyses during the evaluation of a large number of candidate design schemes. Compared with optimization processes driven entirely by finite elements, this method significantly reduces the overall computational cost while maintaining engineering accuracy, making it more suitable for engineering optimization design scenarios involving complex structures.

[0071] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0072] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A smart optimization algorithm for curved fiber optic wings that integrates physical information buckling proxy, characterized in that, include: The quadrilateral region with spatial distortion is equivalently planarized and projected to obtain the two-dimensional physical domain Ω; Define reference domain And for the physical domain Ω and the reference domain Perform a bi-triangular affine mapping to derive the Jacobian matrix; Based on the physical domain Ω and the Jacobian matrix, combined with the curve fiber path parameters, the equivalent stiffness matrix that varies with the space is derived, and the energy functional of the governing equation is mapped from the physical domain to the triangular reference domain to construct the total potential energy functional of the laminate. Based on the total potential energy functional of the laminate, the buckling critical load factor is solved using the physical information neural network PINN to generate the first-order buckling factor and buckling mode. Establish a mapping mechanism from engineering manufacturing constraints to mathematical optimization variables, and construct a global design variable vector; The process can construct constraints, static aeroelastic constraints, dynamic aeroelastic constraints, and strength constraints, and uses an improved adaptive differential evolution algorithm to perform global optimization on the global design variable vector to generate the optimal design variable vector.

2. The intelligent optimization algorithm for curved fiber optic wings that integrates physical information buckling proxy as described in claim 1, characterized in that, The equivalent planarization of the spatially distorted quadrilateral region includes: Calculate the geometric center O of the four vertices of the quadrilateral region Q. Based on the coordinate deviation of each vertex relative to the geometric center O, construct the covariance matrix C, expressed as: In the formula, Let be the global spatial coordinate vector of the i-th vertex of region Q, arranged in counterclockwise order; Let O be the coordinate vector of the geometric center. Eigenvalue decomposition is performed on the covariance matrix C to find the minimum eigenvalue and the corresponding unit eigenvector n. This eigenvalue is then used as the normal to the best-fit plane of region Q to determine the equivalent planar region Q'. The expression is as follows: In the formula, n is the normal to the best-fit plane; v is any eigenvector of the covariance matrix C.

3. The intelligent optimization algorithm for curved fiber optic wings based on physical information buckling proxy as described in claim 1, characterized in that, The projection yields a two-dimensional physical domain Ω, including: The direction of the first side of the quadrilateral is selected as the first basis vector. And normalize it, the expression is: In the formula, Let it be the first basis vector of the local coordinate system; Let be the Euclidean length of the first edge; The second basis vector is generated using the cross product operation. The expression is: First basis vectors Second basis vector By combining the normal to the best-fit plane n, a local right-handed orthogonal basis is generated. ; Construct a local rotation matrix based on the local right-handed orthogonal basis, and translate the vertex coordinates in the global coordinate system to the equivalent geometric center O' to obtain the vertex coordinates of the two-dimensional physical domain Ω. The expression is as follows: In the formula, Let be the two-dimensional local coordinates of the i-th vertex in the physical domain Ω; Let Q' be the coordinates of the i-th vertex of the equivalent planar region Q'. The coordinates of the equivalent center point O'; It is a local rotation matrix.

4. The intelligent optimization algorithm for curved fiber optic wings that integrates physical information buckling proxy as described in claim 1, characterized in that, The defined reference domain And for the physical domain Ω and the reference domain Perform a bi-triangular affine mapping to derive the Jacobian matrix, including: Define the reference field of the standard unit square Along the diagonal, for the physical domain Ω and the reference domain Perform synchronous partitioning to generate two triangular physical domains and two triangular reference domains; In the reference domain Introducing piecewise local parameters A linear transformation is performed on the triangular reference domain to map the irregular triangle to the standard parameter domain. The expression is as follows: ， In the formula, The lower triangular reference domain; The upper triangular reference domain; For global coordinate variables of the reference domain; Construct an affine mapping function from the triangular reference domain to the triangular physical domain, and perform linear interpolation using the vertex coordinates of the physical domain. The expression is as follows: In the formula, Let be the vertex coordinate vector of the physical domain Ω; Let be the coordinate vector of any point within the physical domain Ω; Within the triangular reference domain, the derivation of local parameters is performed respectively. The Jacobian matrix to the physical domain Ω is expressed as: In the formula, Let be the Jacobian matrix of the i-th triangular reference field.

5. The intelligent optimization algorithm for curved fiber optic wings that integrates physical information buckling proxy according to claim 1, characterized in that, The process involves deriving the equivalent stiffness matrix that varies spatially based on the physical domain Ω and the Jacobian matrix, combined with the curved fiber path parameters. This maps the energy functional of the governing equations from the physical domain to the triangular reference domain, constructing the total potential energy functional of the laminate, including: Using the local coordinate system of the physical domain Ω as a basis, the translation method is used to assume that the curve fiber follows... Extending in direction, define the fiber angle function for the k-th layer layup; Based on material engineering constants, a single-layer stiffness matrix of the material principal axis is constructed, and the direction cosine and sine of the local fiber angle are calculated according to the fiber angle function. The single-layer stiffness matrix is ​​transformed to the local coordinate system of the physical domain Ω using the tensor transformation rule to obtain the off-axis stiffness matrix. Integrating the off-axis stiffness matrix along the thickness direction of the laminate yields the in-plane stiffness matrix, the coupling stiffness matrix, and the bending stiffness matrix. Based on Kirchhoff's thin plate theory, the mid-surface displacement field is defined, the relationship between the displacement field and the strain vector and curvature vector is established, and the strain energy density and geometric potential energy density per unit area are constructed by combining the in-plane stiffness matrix, the coupled stiffness matrix and the bending stiffness matrix. Based on the Jacobian matrix and the corresponding inverse matrix, the differential operator on the physical domain Ω is transformed into an operator on the reference domain; Based on the operator transformation and the strain energy density and geometric potential energy density per unit area in the physical domain Ω, the total potential energy functional in the physical domain Ω is mapped to the triangular reference domain to obtain the total potential energy functional of the laminate.

6. The intelligent optimization algorithm for curved fiber optic wings based on physical information buckling proxy as described in claim 1, characterized in that, The method involves using a Physical Information Neural Network (PINN) to solve for the buckling critical load factor based on the total potential energy functional of the laminate, generating a first-order buckling factor and buckling modes, including: Based on the mapping relationship between the physical domain Ω and the triangular reference domain, the integration domain is transformed to the triangular reference domain, and the approximate integral at the sampling points is obtained using the rational calculation method. The expression is as follows: In the formula, N is the total number of sampling points; Let be the local coordinates of the i-th sampling point within the reference domain; The absolute value of the row and column of the Jacobian matrix at the i-th sampling point is used to convert the density values ​​in the reference domain into energy contributions in the physical domain using weighted averages. Let be the strain energy density at the i-th sampling point; Let be the geometric potential energy density at the i-th sampling point; Combining the Rayleigh quotient objective and normalization constraints, the total loss function of the neural network is constructed as follows: In the formula, This is the total loss function of the physical information neural network; This is the set of weights and bias parameters for the neural network. The strain energy density is calculated based on the output displacement field and automatic differential derivative of the neural network. The geometric potential energy density is calculated based on the first derivative of the displacement field output by the neural network. β These are the component weights of the loss function; This is the normalized loss term; The normalized loss term for the buckling modal displacement field is expressed as follows: In the formula, This is the normalized loss term; Let be the lateral deflection value predicted by the neural network at the i-th sampling point; These are the trainable parameters of the neural network; Area weighting factor; By minimizing the total loss function through an optimization algorithm, and after training convergence, the first-order buckling factor and buckling mode are output.

7. The intelligent optimization algorithm for curved fiber optic wings that integrates physical information buckling proxy according to claim 1, characterized in that, The establishment of a mapping mechanism from engineering manufacturing constraints to mathematical optimization variables, and the construction of a global design variable vector, includes: For the upper and lower skin of the wing structure, design domains are constructed separately. These design domains are then divided into non-overlapping discrete sub-regions along the spanwise and chordwise directions of the wing, expressed as: ， In the formula, and These are the design domains for the upper skin and the lower skin, respectively. and These are the r-th sub-regions of the upper and lower skin, respectively; and These represent the number of sub-regions in the upper and lower skin, respectively. For any subregion, the thickness design variable is defined as the total thickness of the subregion. Assuming that the thicknesses of different plies are equal, the expression for the thickness of the k-th ply is: In the formula, The total thickness of the r-th sub-region is a design variable; Let r be the total number of layers in the r-th subdomain; The thickness of a single layer in the k-th layer within the r-th sub-region; The fiber path for each layup is set by The control points are defined, and the expression for the fiber path angle sequence of the k-th layer of the upper and lower skin is as follows: ， In the formula, This is the fiber path angle sequence of the k-th layer of the upper skin; The path angle sequence of the k-th layer of the lower skin; The fiber angle at any position y is obtained by linear interpolation of the angles of adjacent control points, and the expression is: In the formula, The number of control points for a single fiber path; Design variables for the fiber angle of the k-th layer at the k-th control point; It is an interpolation function; Based on the thickness variables of all sub-regions and the fiber angle control point variables of each layup, a unified global design variable vector is constructed, expressed as: In the formula, Design a global variable vector; Design a variable vector for thickness; Design a variable vector for fiber angles; Thickness design variable vector and fiber angle design variable vector The expression is: In the formula, For the r-th sub-region, the k-th layer layup, and the fiber angle at the m-th control point; This represents the total number of skinned sub-regions; Number of layers for the skin; This represents the number of skin control points.

8. The intelligent optimization algorithm for curved fiber optic wings based on physical information buckling proxy as described in claim 1, characterized in that, The process can manufacture constraints, static aeroelastic constraints, dynamic aeroelastic constraints, and strength constraints, respectively including: Manufacturability constraints include: Thickness process constraints, including: For spanwise adjacent sub-regions, the ratio of layer drop thickness to spanwise length must satisfy the spanwise thickness change rate constraint, expressed as: In the formula, Let be the total ply thickness of the i-th and j-th sub-regions, where i is the spanwise partition number and j is the chordwise partition number; The spanwise distance between adjacent spanwise intervals; This represents the maximum allowable rate of change of the spanwise angle in the process. For spanwise adjacent subregions, the ratio of layer drop thickness to spanwise length satisfies the chordal thickness variation rate constraint, expressed as: In the formula, The chordal distance between adjacent chordal intervals; This represents the maximum allowable rate of change of the chordal angle in the process. Angular process constraints, including: For adjacent layers within the same subregion, the fiber angles satisfy the interlayer angle continuity constraint, expressed as: In the formula, For the i-th spanwise region, the j-th tangential region, the k-th layup, and the fiber angle of the m-th grid column; For adjacent grids within the same subregion, the fiber angle satisfies the spanwise angle change rate constraint, expressed as: In the formula, The distance to the center of the neighboring column grid; This represents the maximum allowable rate of angular change in the process. Static aeroelastic constraints include: The wingtip deflection constraint is expressed as follows: In the formula, Let x be the deflection penalty value for individual x; This is the upper limit of the deflection penalty value; This is the normalized deflection over-limit ratio; The nonlinear exponent is used to penalize deflection. This represents the maximum threshold for wingtip deflection. Let x be the wingtip deflection of individual x; The wingtip twist angle constraint is expressed as follows: In the formula, The torsion angle penalty value for individual x; This represents the maximum penalty value for the torsion angle. The normalized torsion angle over-limit ratio; The nonlinear exponent is penalized for the torsion angle; This represents the maximum threshold value for the wingtip twist angle. Let x be the wingtip twist angle of individual x; The derivative of the control surface and the stiffness ratio constraint are expressed as follows: In the formula, The penalty value is the ratio of rudder effect to spring stiffness. For individual x, the rudder-elasticity ratio is given. In the interval The rudder-elasticity ratio after linear normalization. The nonlinear exponent of the rudder-efficiency stiffness ratio penalty value; The upper limit of the penalty value for the rudder effect-to-elasticity ratio; This is the maximum penalty threshold; No penalty threshold; Aeroelastic constraints include: The flutter velocity constraint is expressed as follows: In the formula, This is the flutter speed penalty value; Let x be the flutter velocity of individual x; This is the maximum penalty threshold; No penalty threshold; This represents the upper limit of the flutter speed penalty value; The nonlinear exponent for the flutter velocity penalty value; Strength constraints include: The expression for composite strain constraint is: In the formula, Let x be the strain penalty value for individual x; This represents the upper limit of the strain penalty value; and For control parameters; The total number of units that violate the constraints; To assess the risk level; For design domain; It is a non-linear exponent; The element strain exceeds the limit; The number of layers in a unit; Let be the allowable strain value for the j-th strain; For the first The k-th layer of a unit, and the j-th strain component; Local buckling constraint, expressed as: In the formula, Let x be the buckling penalty value for individual x; This represents the maximum value for the buckling penalty. To assess the degree of risk of buckling; and It is a nonlinear parameter; The total number of regions where buckling occurs; Let be the first-order buckling factor of the j-th region; This is the safety threshold for the buckling factor.

9. The intelligent optimization algorithm for curved fiber optic wings based on physical information buckling proxy as described in claim 1, characterized in that, The step of using an improved adaptive differential evolution algorithm to globally optimize the global design variable vector and generate the optimal design variable vector includes: Minimizing individual fitness as the optimization objective is expressed as: In the formula, For individual fitness; This is the total penalty value; For quality adaptability; The expression for the total penalty value is: In the formula, This is the deflection penalty value; This is the torsion angle penalty value; The penalty value is the ratio of rudder effect to spring stiffness. This is the flutter speed penalty value; This is the strain penalty value; An improved linear population size reduction differential evolution algorithm framework is adopted, which introduces a linear population decay strategy and a three-stage differentiation evolution mechanism. By dynamically adjusting the population size and control parameters, the update strategy is iteratively executed, and the optimal design variable vector is obtained after convergence.

10. The intelligent optimization algorithm for curved fiber optic wings based on physical information buckling proxy as described in claim 9, characterized in that, The improved linear population size reduction differential evolution algorithm framework introduces a linear population decay strategy and a three-stage differentiation evolution mechanism. It iteratively executes an update strategy by dynamically adjusting the population size and control parameters, including: Population linear decay strategies include: The expression for the population in generation g is: In the formula, This is the g-th generation population; This is the i-th individual in the g-th generation population; The total number of individuals in the g-th generation population; This refers to the number of individuals in the first generation population. The number of individuals in the Gth generation population; The maximum number of generations is preset. Using the population mutation strategy under the L-SHADE framework, mutated individual vectors are generated, expressed as follows: In the formula, Let be the mutated individual vector of the i-th individual in the g-th generation population; This is the i-th individual in the g-th generation population; Individuals that rank in the top p% of the g-th generation population; and Two individuals are randomly selected from the population in the g-th generation. This is the difference scaling factor; Vector of mutated individuals With the target vector Perform binomial crossover to generate the test vector, expressed as: In the formula, Let j be the j-th component of the test vector; Let be the crossover probability of the i-th individual; This is a random number corresponding to the weight j of the i-th individual in the g-th generation population; Calculate the test vector fitness value If greater than or equal to the target vector The fitness value, then the trial vector Enter the next generation of the population, and record the differential scaling factor at the same time. and crossover probability Store in a temporary successful set; otherwise, retain the original individual. ; After each generation, the historical mean of the success scaling factor in the historical memory bank is updated using the weighted Lehmer mean and the arithmetic mean. Historical average of the probability of successful crossover .

Images