Gdq adaptive grid point generation method for composite laminate buckling problem

Abstract

The application discloses a GDQ adaptive grid point generation method for a composite laminate buckling problem, and is based on a classical laminate theory to obtain a buckling differential control equation of a symmetrical laminate under complex load, a new discrete point combination form is constructed by introducing a combination parameter alpha, a GDQ adaptive grid node distribution is obtained to solve the differential control equation, and stable and efficient solving of the GDQ method for the composite laminate buckling problem is realized. Through generation of the GDQ adaptive grid point distribution form, the convergence of the critical buckling load calculation result of the composite laminate under the boundary condition containing a free angle point is greatly improved, fast and efficient solving of the GDQ method is realized, the convergence of the calculation result under the boundary condition containing the free angle point is greatly improved, and the research possibility of the GDQ method for solving other types of composite laminate shells in the buckling problem field under complex load and containing a free edge is further expanded.

Description

Technical Field

[0001] This invention relates to buckling treatment technology for composite laminate shells, and particularly to a GDQ adaptive mesh point generation method for buckling problems of composite laminates. Background Technology

[0002] Composite materials possess advantages such as high specific strength, high specific modulus, and designability, and have been widely used in aerospace, rail transportation, wind turbine blades, construction, and bridges. However, under complex service environments, composite laminates are prone to buckling, leading to structural instability. Therefore, the stability of laminates has become a focus of engineering attention. Currently, researchers have conducted detailed studies on the buckling problem of composite laminates under in-plane linearly varying loads using the Rayleigh-Ritz method, the finite difference method, and the generalized differential quadrature (GDQ) method.

[0003] The Generalized Differential Quadrature (GDQ) method is a more stable, efficient, and accurate numerical calculation method (see "Proceedings of the Third International Conference on Advances in Numerical Methods: Engineering: Theory and Applications"—Swansea, UK, 1990: 978-985 and "Solving Two-Dimensional Incompressible Navier-Stokes Equations Using the Generalized Differential Quadrature Method"—International Journal of Numerical Methods for Fluids, 1992, 15(15): 791-798). Its advantages include the ease of calculating weight coefficients without any restrictions, ease of implementation, and its fast convergence and high accuracy, which have led to its widespread use in the vibration or buckling analysis and optimization of beams, rectangular plates, circular plates, etc. The calculation results are in very good agreement with analytical solutions or other numerical methods. The CBCGE method is a new GDQ grid point distribution method (see International Journal of Solids and Structures, 1997, 34(7): 837-846). This method directly couples the boundary conditions with the governing equations and is applicable to simply supported, fixed, and free boundary conditions. It also provides a simple analysis of the convergence of the calculation results when the grid points are increased. However, this paper only conducts vibration analysis on composite laminates and does not analyze the convergence of solving the critical buckling load of laminates under complex load conditions. In 2016, a load perturbation strategy was proposed to address the calculation oscillation and non-convergence phenomena of the GDQ method in solving the buckling problem of composite laminates under in-plane linearly varying load conditions (see “Vibration and Buckling Optimization of Composite Laminates under In-plane Linearly Variation Loads [J]” -- Sun Shiping, Zeng Qinglong, Hu Zheng. Journal of Composite Materials, 2016, 33(12): 2860-2868.). However, this perturbation strategy is not very effective in solving buckling problems with complex loads such as free edges and shear, and needs further improvement.

[0004] The GDQ and CBCGE methods are characterized by fast convergence, high accuracy, and high efficiency in solving vibration problems of composite laminates. Furthermore, the GDQ and CBCGE methods do not suffer from poor convergence of calculation results when solving vibration problems. However, when solving buckling problems, especially buckling problems with complex in-plane loads and free corners, the convergence is not very good and needs further improvement. Summary of the Invention

[0005] The purpose of this invention is to improve the computational oscillation non-convergence problem existing in the GDQ method and CBCGE method when solving buckling problems under complex load conditions and free boundary conditions. Based on the discrete point perturbation strategy, a new discrete point combination form is proposed, which introduces combination parameters. α This generates an adaptive mesh node distribution, avoiding symmetrical nodes that could lead to very close-to-zero numbers on the diagonal of the load matrix in subsequent calculations. This enables fast and efficient solutions for GDQ, further expanding the possibilities of using the GDQ method to solve buckling problems of other types of composite laminated plates and shells under complex loads and with free edges.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: The GDQ adaptive mesh point generation method for buckling problems of composite laminates comprises the following steps: Step 1: Based on classical laminate theory, the buckling differential governing equations of symmetrical laminates under complex loads are obtained; S01: Buckling differential governing equation of laminated plate: (1) In the formula: w(x,y,t) is the displacement function of the mid-surface. t For time, N x and N y The laminates are respectively in x and y The in-plane unit load borne in the direction, N xy For laminates in xy The unit shear load borne in the plane. For stiffness matrix elements, Its expression is: (2) In the formula: Z k For the first k The distance from the lower surface of the ply to the mid-plane of symmetry. For the first k Conversion stiffness factor for ply; S02: The in-plane linearly varying load borne by the composite laminate can be expressed as: (3) In the formula: N0 is the load amplitude, and the parameter β Describe different types of in-plane loads, such as β= At time 0, the load is uniformly distributed. β= When 1 is a linear load, β= At time 2, the load is a bending load, and the load is a combined load. N x :N y :N xy =1:0:1( β =0) indicates that the laminate is subjected to... x The combined action of a uniformly distributed load and a shear load; by x Taking the edge = 0 as an example, the boundary control equations are as follows: Fixed boundary equations: (4) Simply supported boundary equations: (5) Free boundary equations: (6) For free corner points on two adjacent free edges, additional constraint equations are required: (7) Step 2: Solving the differential equation using the GDQ method: S01: In the GDQ method, the function nth derivative It can be approximated by weighting the function values ​​at N discrete points within the design domain: (8) Although: N x is the number of discrete points in the domain. k ( k =1,..., N ) is the first k discrete points, These are the nth-order weighted coefficients; S02: The GDQ method uses Lagrange polynomials. To define the weighted system The formula for its calculation is: (9) In the formula: Let be the first derivative of the Lagrange polynomial. N The number of discrete points; S03: For a bivariate function g(x,y), its derivative can be expressed by a linear weighted combination as follows: (10) In the formula: N, M Within the domain x , y The number of discrete points in the direction, ( x i , y j ) ( i =1,..., N , j =1,..., M () represents the coordinates of a discrete point. s , q They are respectively x , y Order of partial differential in direction They are respectively x , y Direction s , q Order weighting coefficients; S04: Let and take the displacement function Normalizing the control equation (1) gives: (11) S05: Let the aspect ratio of the laminate be λ=ɑ / b. Substituting equation (10) into equation (11) yields the expression for the discrete equation system: (12) S06: Substituting equation (10) into the boundary constraint equations (4) to (7) and normalizing them, we get: Fixed boundary: (13) Simply supported boundary: (14) Free boundary: (15) Free corner constraint: (16) S07: Record in the plane of the laminate N × M The displacement vector of each discrete point is ,in W b Indicates the boundary and its adjacent (4) N +4 M The known displacement vectors of -16) discrete points W d Indicates the interior ( N -4)×(M -4) unknown displacement vectors of discrete points; S08: Boundary conditions are applied by replacing the corresponding equations in the equation set (12) with boundary equations. When a free corner point appears at an adjacent position of multiple free edges, the equation corresponding to the discrete point where the corner point is located is replaced with the constraint equation (16), resulting in a vector form of the equation set and the corresponding characteristic equation: (17) (18) In the formula: ζ The buckling load factor is... When there is no free edge under load A bb 4 N +4( M -4) Order square matrix A dd and B dd for( N -4)×( M -4) Order square matrix A db and B db for( N × M -4 N -4 M +16)×(4 N +4 M -16) order matrix; S09: Solve the characteristic equation (18) to obtain the minimum eigenvalue and its corresponding critical buckling load. N cr ,Will N cr Normalization was performed to obtain the dimensionless critical buckling load factor. k * : (19) In the formula: E2 is the tensile and compressive elastic modulus of the ply material in the 2-direction, and v is Poisson's ratio; Step 3: Convergence Problem of GDQ Method and its Improvements: S01: GDQ uses the Gauss-Lobatto-Chebyshev expression shown in equation (20) to obtain the discrete point distribution in the two-dimensional design domain: (20) S02: Furthermore, the discrete point position perturbation strategy is used to improve the singularity of the characteristic equation matrix to improve the computational accuracy, as shown in Equation (21): (twenty one) In the formula: δ is the disturbance parameter; S03: The discrete point position perturbation strategy results in non-convergent oscillations for the calculation of boundary combination problems with free corner points, and is sensitive to the discrete point distribution; therefore, a new discrete point distribution combination form is constructed as shown in equation (22): (twenty two) In the formula: α When there are free corner points, the parameters are combined. α =0.5, otherwise take 1; α The smaller the value, the closer the new discrete points are to the boundary and the denser they are. α The closer the value is to 1, the more similar the new discrete point distribution will be to the Gauss-Lobatto-Chebyshev discrete points. Choosing a suitable... α This allows for adjusting the distance between boundary discrete points and their adjacent discrete points, reducing computational instability caused by excessive deformation differences.

[0007] This invention significantly improves the convergence of calculation results under boundary conditions containing free corners by generating an adaptive mesh point distribution, achieving fast and efficient solutions for GDQ. It also greatly improves the convergence of calculation results under boundary conditions containing free corners and further expands the research possibilities of the GDQ method in solving buckling problems of other types of composite laminated plates and shells under complex loads and with free edges. Attached Figure Description

[0008] Figure 1a This is a schematic diagram of a composite laminate plate with symmetrical composite load in an embodiment of the present invention; Figure 1b This is a schematic diagram of the symmetrical layup structure of the composite material laminate in an embodiment of the present invention; Figure 2 This is the in-plane load of the edge x=0 in this embodiment of the invention. N x Schematic diagram; Figure 3 This is a distribution diagram of discrete points within the two-dimensional design domain of the GDQ method in this embodiment of the invention; Figure 4 This is an embodiment of the present invention. N A schematic diagram of the discrete point distribution obtained by the three methods when taking 15; Figure 5 This is an embodiment of the present invention. α Discrete points for different values X 2 takes values ​​as N The trend of increase; Figure 6Under the CSCS boundary conditions in the embodiments of the present invention N Convergence curves of critical buckling load for laminated plates with different values; Figure 7 Under the CSCF boundary conditions in the embodiments of the present invention N Convergence curves of critical buckling load for laminated plates with different values; Figure 8 Under the CFSF boundary conditions in the embodiments of the present invention N Convergence curves of critical buckling load for laminated plates with different values; Figure 9 Under the CSFF boundary conditions in the embodiments of the present invention N Convergence curves of critical buckling load for laminated plates with different values; Figure 10 Under the CSFF boundary conditions in the embodiments of the present invention N The relative error diagram between the critical buckling load of the laminated plate obtained by the three methods and the FEM calculation results for different values; Figure 11 Under the CFFF boundary conditions in the embodiments of the present invention N Convergence curves of critical buckling load for laminated plates with different values; Figure 12 Under the CFFF boundary conditions in the embodiments of the present invention N Error diagram of the critical buckling load of laminated plates calculated by the three methods and the FEM calculation results when different values ​​are taken. Detailed Implementation

[0009] The present invention will be further described below with reference to the accompanying drawings and embodiments. See also Figures 1a to 12 A GDQ adaptive mesh point generation method for buckling problems of composite laminates, the steps of which are as follows: Step 1: Based on classical laminate theory, the buckling differential governing equations of symmetrical laminates under complex loads are obtained.

[0010] S01: Composite laminates subjected to complex loads with symmetrical properties (e.g., Figure 1a and Figure 1b (As shown) Length a ,Width b ,thick h θ k For the first k Fiber angle of laminated layers. Based on classical laminate theory, the buckling differential governing equation of laminate is: (1) In the formula, w(x,y,t) is the displacement function of the mid-surface. t For time, N x and N y The laminates are respectively inx and y The in-plane unit load borne in the direction, N xy For laminates in xy The unit shear load borne in the plane. The elements of the stiffness matrix are expressed as follows: (2) In the formula Z k For the first k The distance from the lower surface of the ply to the mid-plane of symmetry. For the first k Conversion stiffness factor for ply layup.

[0011] S02: The in-plane linearly varying load borne by the composite laminate can be expressed as: (3) In the formula, N0 is the load amplitude, and the parameter β Describe different types of in-plane loads, such as: β= At time 0, the load is uniformly distributed. β= When 1 is a linear load, β= At time 2, it is a bending load; combined load. N x :N y :N xy =1:0:1( β =0) indicates that the laminate is subjected to... x Combined action of uniformly distributed load and shear load (e.g.) Figure 2 (As shown).

[0012] S03: with x Taking the edge = 0 as an example, the boundary control equations are as follows: Fixed boundary equations: (4) Simply supported boundary equations: (5) Free boundary equations: (6) For free corner points on two adjacent free edges, additional constraint equations are required: (7) Step 2: Solve the differential equation using the GDQ method.

[0013] S01: In the GDQ method, the function nth derivative Within the design domain NThe function values ​​at discrete points are weighted to approximate the expression: (8) In the formula N x is the number of discrete points in the domain. k ( k =1,..., N ) is the first k discrete points, for n The order weighting coefficient.

[0014] S02: The GDQ method uses Lagrange polynomials. Let's define the weighted system, and its calculation formula is: (9) In the formula: Let be the first derivative of the Lagrange polynomial. N The number of discrete points.

[0015] S03: For a bivariate function g(x,y), its derivative can be expressed by a linear weighted combination as follows: (10) In the formula N, M Within the domain x , y The number of discrete points in the direction, ( x i , y j ) ( i =1,..., N , j =1,..., M () represents the coordinates of a discrete point. s , q They are respectively x , y Order of partial differential in direction They are respectively x , y Direction s , q The order weighting coefficient.

[0016] S04: Let and take the displacement function Normalizing the control equation (1) gives: (11) S05: Let the aspect ratio of the laminate be λ=ɑ / b. Substituting equation (10) into equation (11) yields the expression for the discrete equation system: (12) S06: Substituting equation (10) into the boundary constraint equations (4) to (7) and normalizing them, we get: Fixed boundary: (13) Simply supported boundary: (14) Free boundary: (15) Free corner constraint: (16)

[0017] S07: Record in the plane of the laminate N × M The displacement vector of each discrete point is W b Indicates the boundary and its adjacent (4) N +4 M The known displacement vectors of -16) discrete points, W d Indicates the interior ( N -4)×( M -4) unknown displacement vectors of discrete points. For the case of a free edge under load, to ensure that the boundary deformation is compatible with the load, the displacements of discrete points adjacent to the free edge are included as unknowns in W. d When free corners exist, the displacement of discrete points at the free corner positions is included as an unknown quantity in W. d .

[0018] S08: Boundary conditions are applied by replacing the corresponding equations in the equation set (12) with boundary equations. When a free corner point appears at an adjacent position of multiple free edges, the equation corresponding to the discrete point where the corner point is located is replaced with the constraint equation (16), resulting in a vector form of the equation set and the corresponding characteristic equation: (17) (18) In the formula ζ The buckling load factor is... When there is no free edge under load A bb 4 N +4( M -4) Order square matrix A dd and B dd for( N -4)×( M -4) Order square matrix A db and Bdb for( N × M -4 N -4 M +16)×(4 N +4 M -16) order matrix; S09: Solve the characteristic equation (18) to obtain the minimum eigenvalue and its corresponding critical buckling load. N cr ,Will N cr Normalization is performed to obtain the dimensionless critical buckling load factor k. * : (19) In the formula, , E 2 represents the tensile and compressive modulus of elasticity of the layup material in direction 2. v It is Poisson's ratio.

[0019] Step 3: Convergence Problem of GDQ Method and Its Improvements S01: The GDQ method uses the Gauss-Lobatto-Chebyshev expression shown in equation (20) to obtain the discrete point distribution in the two-dimensional design domain. The discrete points are denser closer to the boundary and sparser farther away from the boundary, which can obtain more accurate calculation results with fewer grid nodes. However, the GDQ method has poor calculation accuracy and non-convergence of oscillations when solving the buckling problem of plates under complex loads and the vibration problem of plates with free edges.

[0020] (20) S02: Furthermore, the discrete point position perturbation strategy is used to improve the singularity of the characteristic equation matrix to improve the computational accuracy, as shown in Equation (21): (twenty one) In the formula: δ=10 -4 These are the disturbance parameters.

[0021] S03: The discrete point perturbation strategy can stably and quickly obtain the convergent solution of the boundary combination problem without free corners, but the calculation results for the boundary combination problem with free corners oscillate and do not converge, and are sensitive to the discrete point distribution. The reason is that the two constraint equations applied to the free corners have excessive calculation errors. To this end, a new discrete point distribution combination scheme as shown in Equation (22) is constructed to generate an adaptive grid node distribution (e.g., Figure 3 As shown in the figure, this method reduces the computational error caused by the multi-constraint equations of discrete points by reducing the distance between the boundary discrete points and their adjacent discrete points.

[0022] (twenty two) In the formula, α These are the combination coefficients. When free corner points exist... α= 0.5, otherwise take 1. N =15、 α Discrete points plotted with =0.5 X i Distribution (e.g.) Figure 4 Show), Figure 4 "×" represents Discrete points are indicated by "+" for Gauss-Lobatto-Chebyshev discrete points and "o" for new adaptive discrete points. Adjustment... α and N Obtain discrete points X The distribution of 2 (e.g.) Figure 5 (shown), with α As the number of discrete points decreases, the new adaptive discrete point distribution becomes closer to the boundary and denser. By selecting appropriate... α Adjust the distances between the boundary discrete points and their adjacent discrete points to avoid computational instability.

[0023] Example: The following is a specific example of the present invention: With length a Taking a composite laminate with a diameter of 500mm and an aspect ratio of λ=1 as an example, its layup sequence is [30 / -45 / 60 / -75]. s The layup material is T300 / 5208 epoxy resin-based carbon fiber composite material, and the layup thickness is... t 0=0.125mm, the properties of the layup material are shown in Table 1. E 1. E 2. E 3 and G 12 , G 23 , G 13 These are the elastic modulus and shear modulus of the composite material, respectively. v 12 , v 23 , v 13 Poisson's ratio, ρ Density. The load-bearing mode of the laminate is... N x :N y :N xy =1:1:1( β =2), where 1 indicates that the load on the laminate is a unit load (-1). NThe buckling convergence of composite laminates under complex loads and different boundary conditions was calculated using three methods: GDQ, CBCGE, and IGDQ. CBCGE uses the following formula to obtain the distribution of discrete points: (twenty three) In the formula, .

[0024] The relative error is calculated using equation (24): (twenty four) ; Step 1: Consider the composite load borne by the symmetrical laminate [30 / -45 / 60 / -75]s: Nx : Ny : Nxy =1:1:1( β =2), analyzed the perturbation quantity under boundary conditions without free edges. δ =10 -4 The error between the solutions from the three methods and the FEM solution.

[0025] S01: Proceed N x :N y :N xy =1:1:1( β =2) Calculation and solution of buckling performance of composite laminates under CSCS boundary conditions under load, yielding three methods. N Convergence of calculations for the buckling performance of laminates with different values ​​(e.g.) Figure 6 As shown in the figure, the error between the result and the finite element calculation result is calculated.

[0026] ; According to Table 2 and Figure 6 It can be concluded that for the convergent solution of the boundary combination problem without free edges under combined load, as the number of grid points increases, the error between the solutions of GDQ and IGDQ and the solution of finite element method gradually decreases, showing better convergence performance, while the solution of CBCGE oscillates and does not converge.

[0027] Step 2: Consider the composite load borne by the symmetrical laminate [30 / -45 / 60 / -75]s: Nx : Ny : Nxy =1:1:1( β =2), analyzed the perturbation quantity under the boundary condition containing free edges but no free corners. δ =10 -4 The error between the solutions from the three methods and the FEM solution.

[0028] S01: Proceed N x :N y :N xy =1:1:1( β =2) Calculation and solution of buckling performance of composite laminates under CSCF boundary conditions with load, obtaining three methods in... N Convergence of calculations for the buckling performance of laminates with different values ​​(e.g.) Figure 7 As shown in the figure, the error between the result and the finite element calculation result is calculated.

[0029] ; According to Table 3 Figure 7 It is found that for the convergent solution of the boundary combination problem with one free edge and no free corner points under the action of composite load, as the number of grid points increases, the error between the solution of GDQ and IGDQ and the solution of finite element method gradually decreases, showing better convergence performance, while the solution of CBCGE oscillates and does not converge.

[0030] S02: Proceed N x :N y :N xy =1:1:1( β =2) Calculation and solution of buckling performance of composite laminates under CFSF boundary conditions under load, obtaining three methods in... N Convergence of calculations for the buckling performance of laminates with different values ​​(e.g.) Figure 8 As shown in the figure, the error between the result and the finite element calculation result is calculated.

[0031] ; According to Table 4 Figure 8 It is found that for the convergent solution of the boundary combination problem with two free edges but no free corners under combined load, as the number of grid points increases, the error between the solutions of GDQ and IGDQ and the solution of finite element method gradually decreases, showing better convergence performance, while the solution of CBCGE oscillates and does not converge.

[0032] Step 3: Consider the composite load borne by the symmetrical laminate [30 / -45 / 60 / -75]s: N x : N y : N xy =1:1:1( β=2), and analyzed the perturbation quantity under the boundary conditions containing free edges and free corners. δ =10 -4 The error between the solutions from the three methods and the FEM solution.

[0033] S01: Proceed N x :N y :N xy =1:1:1( β =2) Calculation and solution of buckling performance of composite laminates under CSFF boundary conditions under load, obtaining three methods in... N Convergence of calculations for the buckling performance of laminates with different values ​​(e.g.) Figure 9 As shown), and the error between the result and the finite element calculation result is calculated (e.g. Figure 10 (As shown).

[0034] ; According to Table 5 and Figure 9 It can be derived that the convergent solution for the boundary combination problem with one free corner point under combined loads is obtained. However, as the number of grid points increases, the convergence of the critical buckling load solution of the GDQ becomes poor. According to... Figure 10 It can be concluded that the critical buckling load solution obtained by the IGDQ method has higher accuracy and better convergence than the solutions obtained by GDQ and CBCGE.

[0035] S02: Proceed N x :N y :N xy =1:1:1( β =2) Calculation and solution of buckling performance of composite laminates under CFFF boundary conditions with load, obtaining three methods in... N Convergence of calculations for the buckling performance of laminates with different values ​​(e.g.) Figure 11 As shown in the figure, the error between the result and the finite element calculation result is calculated.

[0036] ; According to Table 6 and Figure 11 It can be derived that for the convergent solution of the boundary combination problem with two free corners under combined loads, the convergence of the critical buckling load solutions of GDQ and CBCGE is poor as the number of grid points increases. Figure 12 It is evident that the critical buckling load solution obtained by the IGDQ method has higher accuracy and better convergence compared to the solutions obtained by GDQ and CBCGE.

Claims

1. A GDQ adaptive mesh point generation method for buckling problems of composite laminates, characterized in that, The steps are as follows: Step 1: Based on classical laminate theory, the buckling differential governing equations of symmetrical laminates under complex loads are obtained; S01: Buckling differential governing equation of laminated plate: （1） In the formula: w(x,y,t) is the displacement function of the mid-surface. t For time, N x and N y These are the in-plane unit loads, N, borne by the laminate in the x and y directions, respectively. xy The unit shear load borne by the laminate in the xy plane The elements of the stiffness matrix are expressed as follows: （2） In the formula: z k Let be the distance from the lower surface of the k-th ply to the mid-plane of symmetry. The conversion stiffness factor for the k-th ply; S02: The in-plane linearly varying load borne by the composite laminate can be expressed as: （3） In the formula: N0 is the load amplitude, and the parameter β describes different types of in-plane loads. β When the load is 0, it is a uniformly distributed load. β When =1, it is a linear load. β When the value is 2, it is a bending load, a combined load. N x :N y :N xy =1:0:1, β =0 indicates that the laminate is subjected to a combination of uniformly distributed load and shear load in the x-direction. When x=0, the boundary control equations are as follows: Fixed boundary equations: （4） Simply supported boundary equations: （5） Free boundary equations: （6） For free corner points on two adjacent free edges, additional constraint equations are required: （7） Step 2: Solving the differential equation using the GDQ method: S01: In the GDQ method, the function nth derivative It can be approximated by weighting the function values ​​at N discrete points within the design domain: （8） In the formula: N is the number of discrete points in the domain, x k For the first k There are discrete points, k=1,...,N. These are the nth-order weighted coefficients; S02: The GDQ method uses Lagrange polynomials. To define the weighting coefficients The formula for its calculation is: （9） In the formula: is the first derivative of the Lagrange polynomial, and N is the number of discrete points; S03: For a bivariate function g(x,y), its derivative can be expressed by a linear weighted combination as follows: （10） In the formula: N , M These represent the number of discrete points in the x and y directions within the domain, respectively. xi , yj () represents the coordinates of a discrete point. i =1,..., N , j =1,..., M,s , q Let be the partial differential orders in the x and y directions, respectively. respectively in the x and y directions s , q Order weighting coefficients; S04: Let and take the displacement function Normalizing the control equation (1) gives: （11） S05: Let the aspect ratio of the laminate be λ=ɑ / b. Substituting equation (10) into equation (11) yields the expression for the discrete equation system: （12） S06: Substituting equation (10) into the boundary constraint equations (4) to (7) and normalizing them, we get: Fixed boundary: （13） Simply supported boundary: （14） Free boundary: （15） Free corner constraint: （16） S07: Let the displacement vector of N×M discrete points in the plane of the laminate be... W b Indicates the boundary and its adjacent (4) N +4 M The known displacement vectors of -16) discrete points, W d Indicates the interior ( N -4)×( M -4) unknown displacement vectors of discrete points; S08: Boundary conditions are applied by replacing the corresponding equations in the equation set (12) with boundary equations. When a free corner point appears at an adjacent position of multiple free edges, the equation corresponding to the discrete point where the corner point is located is replaced with the constraint equation (16), resulting in a vector form of the equation set and the corresponding characteristic equation: （17） （18） In the formula: ζ is the buckling load coefficient. When there is no free edge under load A bb 4 N +4( M -4) Order square matrix A dd and B dd for( N -4)×( M -4) Order square matrix A db and B db for( N × M -4 N -4 M +16)×(4 N +4 M -16) order matrix; S09: Solve the characteristic equation (18) to obtain the minimum eigenvalue and its corresponding critical buckling load. N cr ,Will N cr Normalization is performed to obtain the dimensionless critical buckling load factor k. * : （19） In the formula: E2 is the tensile and compressive elastic modulus of the ply material in the 2-direction, and v is Poisson's ratio; Step 3: Convergence Problem of GDQ Method and its Improvements: S01: GDQ uses the Gauss-Lobatto-Chebyshev expression shown in equation (20) to obtain the discrete point distribution in the two-dimensional design domain: （20） S02: Furthermore, the discrete point position perturbation strategy is used to improve the singularity of the characteristic equation matrix to improve the computational accuracy, as shown in Equation (21): （21） In the formula: ( i =1,2,…, N ), j =1,2,…, M ), δ These are the disturbance parameters; S03: The discrete point position perturbation strategy results in non-convergent oscillations for the calculation of boundary combination problems with free corner points, and is sensitive to the discrete point distribution; therefore, a new discrete point distribution combination form is constructed as shown in equation (22): （22） In the formula: α When there are free corner points, the parameters are combined. α =0.5, otherwise take 1; α The smaller the value, the closer the new discrete points are to the boundary and the denser they are. α The closer the new discrete point distribution is to 1, the more it will converge to the Gauss-Lobatto-Chebyshev discrete point distribution; choosing an appropriate... α This allows for adjusting the distance between boundary discrete points and their adjacent discrete points, reducing computational instability caused by excessive deformation differences.

Images