A method for optimizing parameters of a geopolymer composite foundation based on finite element analysis

Abstract

The present application belongs to the field of geotechnical engineering numerical analysis and ground treatment, and relates to a geopolymer composite foundation parameter optimization method based on finite element analysis. It comprises the following steps: step 1, initial mesh partitioning is performed on the calculation domain of the geopolymer mixing pile-soil composite foundation and recorded as a unified reference mesh, and the response field is mapped to the unified reference mesh to form a unified dimension response snapshot vector; step 2, singular value decomposition is performed on the snapshot matrix composed of all unified dimension response snapshot vectors to extract an intrinsic orthogonal matrix, and an order-reduced fast evaluation model is formed by interpolation mapping and the intrinsic orthogonal matrix; step 3, the order-reduced fast evaluation model is used to sample and evaluate the value range of the parameterized design variables, and the optimal design parameter interval of the key design parameters is output. The Borgonovo and Shapley joint sorting simultaneously captures the fair distribution of the distribution shape influence and the interaction effect, and provides a parameter optimization scheme for the geopolymer composite foundation which takes into account the accuracy and robustness.

Description

Technical Field

[0001] This invention belongs to the field of numerical analysis and foundation treatment technology in geotechnical engineering, specifically relating to a method for optimizing parameters of geopolymer composite foundations based on finite element analysis. Background Technology

[0002] Geopolymer mixing pile composite foundations are a novel reinforcement method in the field of soft soil foundation treatment. They involve mixing geopolymer cementitious materials with in-situ soil to form reinforced piles, which, together with the soil between the piles, bear the load from the superstructure. Compared with traditional cement mixing piles, geopolymer mixing piles have advantages such as lower carbon emissions and faster early strength development, and have received increasing attention in engineering practice in recent years. However, geopolymer composite foundations involve numerous design parameters, including pile length, pile diameter, pile spacing, geopolymer material stiffness, and soil moisture content. These parameters have interdependent effects on bearing capacity and settlement, making it difficult to obtain a globally optimal design scheme by relying solely on empirical formulas or single-factor analysis.

[0003] Numerical simulation is an important tool for studying the mechanical behavior of composite foundations. Currently, a common approach is to use the Cauchy continuum-based finite element method to model and analyze pile-soil systems. However, the Cauchy continuum model has inherent limitations in describing the shear band at the pile-soil interface: the shear band width decreases with mesh refinement and tends towards zero, making the calculation results highly sensitive to mesh size and difficult to obtain objective and reliable numerical solutions. Furthermore, most existing finite element analyses use fixed meshes, which cannot adaptively adjust for stress and deformation concentration areas during loading. This results in either high computational costs due to global mesh refinement or insufficient local accuracy affecting the reliability of the results.

[0004] In parameter sensitivity analysis, commonly used methods include the Sobol index based on variance decomposition and the Morris screening method. The Sobol index relies on the variance decomposition assumption and can only reflect the parameter's contribution to the response variance, failing to capture the parameter's influence on the overall shape of the response probability distribution. The Morris screening method is a qualitative screening tool and cannot provide a quantitative sensitivity ranking. Furthermore, neither of these methods considers the fair allocation of interaction effects between parameters, making it difficult to accurately identify the truly critical design parameters when multiple parameters have strong coupling effects. In addition, the long single-run solution time of a full-order finite element model makes the computational cost prohibitively high in sensitivity analysis and optimization searches that require evaluating a large number of parameter combinations. Summary of the Invention

[0005] The main objective of this invention is to provide a parameter optimization method for geopolymer composite foundations based on finite element analysis. By using Borgonovo and Shapley joint sorting, the method simultaneously captures the fair allocation of the influence of distribution shape and interaction effects, providing a parameter optimization scheme that balances accuracy and robustness for geopolymer composite foundation design.

[0006] To address the aforementioned technical problems, this invention provides a method for optimizing geopolymer composite foundation parameters based on finite element analysis, comprising the following steps: Step 1: Initially mesh the computational domain of the geopolymer mixing pile-soil composite foundation and record it as a unified reference mesh. Set translational degrees of freedom and micro-polar rotational degrees of freedom on the mesh nodes and adopt the Cosserat micro-polar constitutive relation within the element. Set parameterized design variables and generate several sets of parameter combinations. Perform incremental load loading on each set of parameter combinations. After each load increment step converges, perform node redistribution operation based on the rotational gradient index and perform element subdivision operation based on the plastic strain index. Extract the pile top bearing capacity, total settlement of the composite foundation and differential settlement corresponding to each set of parameter combinations, and map the response field to the unified reference mesh to form a unified dimension response snapshot vector. Step 2: Perform singular value decomposition on the snapshot matrix composed of all uniform dimension response snapshot vectors to extract the intrinsic orthogonal basis matrix. Project the uniform dimension response snapshot vector of each parameter combination onto the intrinsic orthogonal basis matrix to obtain the reduced order coefficient vector. Construct an interpolation mapping from parameterized design variables to the reduced order coefficient vector. The reduced order fast evaluation model is composed of the interpolation mapping and the intrinsic orthogonal basis matrix. Step 3: Use the reduced-order rapid evaluation model to sample and evaluate the value range of the parameterized design variables, calculate the Borgonovo sensitivity index and Shapley interaction assignment value of each parameterized design variable, multiply them to obtain the comprehensive sensitivity ranking value, select key design parameters, traverse and evaluate within the value range of the key design parameters and screen feasible solutions that meet the preset bearing capacity constraints and preset settlement constraints, and output the optimal design parameter range of the key design parameters.

[0007] Furthermore, in step 1, the initial mesh is divided into an initial tetrahedral mesh, dividing the geopolymer mixing pile body, the surrounding soil, and the bearing layer at the pile tip into three material subdomains; each mesh node is set with three translational degrees of freedom and three independent micropolar rotational degrees of freedom; the parametric design variables are pile length, pile diameter, pile spacing, geopolymer elastic modulus, and soil moisture content, totaling five; several sets of parameter combinations are generated by Latin hypercube sampling within the value range of the five parametric design variables; each load increment step is solved until convergence through Newton-Raphson iteration.

[0008] Furthermore, in step 1, the Cosserat micropole constitutive relation establishes a mapping between the couple stress tensor and the micropole rotational gradient tensor by introducing an intrinsic characteristic length parameter, so that the numerical solution of the shear zone at the pile-soil interface depends on the intrinsic characteristic length parameter rather than the mesh size.

[0009] Further, in step 1, the node redistribution operation includes: traversing each cell in the current mesh, extracting the micro-polar rotation gradient tensor at each Gaussian integration point of the cell, and taking the maximum value of the Frobenius norm as the rotation gradient index of the cell; calculating the arithmetic mean of the rotation gradient indices of all cells in the global domain; marking cells whose rotation gradient indices exceed the threshold obtained by the product of the arithmetic mean and the preset r-type amplification factor as r-type candidate cells; for each node of each r-type candidate cell, translating the node along the direction of the fastest decrease in the rotation gradient index value, with the translation step size being the product of the shortest cell side length connected to the node and the preset translation ratio factor; after each node translation is completed, checking the Jacobian determinant value of all cells associated with the node, and if the Jacobian determinant value of any cell is less than zero, then the node is reverted to its position before translation.

[0010] Furthermore, in step 1, the element subdivision operation includes: traversing each element in the current mesh, extracting the equivalent plastic strain at each Gaussian integration point of the element and taking the maximum value as the element's plastic strain index; calculating the arithmetic mean of the plastic strain indices of all elements in the global domain; marking elements whose plastic strain indices exceed the threshold obtained by the product of the arithmetic mean and the preset h-type magnification factor as h-type candidate elements; for each h-type candidate element, performing one mid-surface bisection along the principal direction of the equivalent plastic strain gradient within the h-type candidate element to generate two sub-elements; for each suspended node generated by the subdivision, constraining the translational degrees of freedom of the suspended node to the linear interpolation of the translational degrees of freedom of two adjacent non-suspended nodes, and simultaneously constraining the micro-polar rotational degrees of freedom of the suspended node to the linear interpolation of the micro-polar rotational degrees of freedom of two adjacent non-suspended nodes.

[0011] Furthermore, in step 1, the unified reference grid is the node coordinates and topological relationships of the initial grid; the pile top bearing capacity is obtained by summing the vertical reactions of all boundary nodes on the pile top surface; the total settlement of the composite foundation is the vertical displacement of the center node on the pile top surface; the differential settlement is the maximum value of the vertical displacement difference between the center node on the pile top surface and the preset observation node on the soil surface between the piles; the response field is mapped to the unified reference grid by mapping the global node displacement field and micro-polar rotation field obtained by combining each set of parameters on the adaptive grid to the nodes of the unified reference grid through inverse distance interpolation.

[0012] Furthermore, in step 2, the singular value decomposition obtains the left singular vector sequence and the corresponding singular value sequence. Starting from the first order, the squares of the singular values ​​are accumulated step by step until the proportion of the accumulated result to the total sum of the squares of all singular values ​​reaches the preset energy retention threshold for the first time. The corresponding first few left singular vectors are then grouped into an eigenorthogonal basis matrix.

[0013] Furthermore, in step 2, the interpolation mapping is a radial basis function interpolation mapping, constructed with the parameterized design variable values ​​corresponding to all parameter combinations as input nodes and the corresponding reduced-order coefficient vectors as output nodes. For a new combination of parameterized design variables, the predicted reduced-order coefficient vector is obtained through radial basis function interpolation mapping, and then the predicted reduced-order coefficient vector is linearly combined and restored to the global response field on the unified reference grid through the intrinsic orthogonal basis matrix. The predicted value is obtained in the same way as the pile top bearing capacity, total settlement of composite foundation and differential settlement in step 1.

[0014] Furthermore, in step 3, the sampling evaluation involves uniformly sampling within the value range of the five parametric design variables at a preset sampling density, and obtaining the predicted values ​​of pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each sampling point using a reduced-order rapid evaluation model. The calculation of the Borgonovo sensitivity index includes: for the i-th parametric design variable, dividing the value range of the i-th parametric design variable into several intervals according to a preset number of segments; within each interval, collecting the predicted values ​​of the mechanical response corresponding to all sampling points falling within the interval. The conditional probability density curve of the mechanical response within the interval is fitted using the Gaussian kernel density estimation method. Simultaneously, the unconditional probability density curve is fitted to the predicted mechanical response values ​​corresponding to all sampling points using the same method. Within each interval, the smaller of the conditional and unconditional probability density curve values ​​is taken point-by-point along the range of the mechanical response, and the integral is performed. The local offset of the interval is obtained by subtracting 1 from the integral result. The arithmetic mean of the local offsets of all intervals is used to obtain the Borgonovo sensitivity index of the i-th parameterized design variable.

[0015] Furthermore, in step 3, the calculation of the Shapley interaction allocation value includes: treating all parameterized design variables as participants in a cooperative game; for the i-th parameterized design variable, enumerating all parameter subsets that do not contain the i-th parameterized design variable; for each parameter subset, fixing the parameterized design variables within the parameter subset to their respective discrete sampled values ​​in a preset sample set, and iterating through all discrete sampled values ​​of the parameterized design variables outside the parameter subset in the preset sample set, calculating the variance of the predicted mechanical response value as the characteristic function value of the parameter subset; after adding the i-th parameterized design variable to the parameter subset, calculating the characteristic function value of the expanded subset in the same way, and using the difference between the characteristic function value of the expanded subset and the characteristic function value of the parameter subset as the marginal contribution under the parameter subset; determining the corresponding Shapley combination weight according to the size of the parameter subset, and then applying the formula to all parameter subsets. The sum of the marginal contributions after weighted Shapley combination yields the Shapley interaction assignment value of the i-th parameterized design variable. The key design parameters are the top three parameterized design variables selected after arranging all parameterized design variables in descending order of comprehensive sensitivity. The optimal design parameter interval is determined by performing an equally spaced grid traversal within the value range of the three key design parameters, evaluating the pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each grid traversal point using a reduced-order rapid evaluation model, and selecting all feasible traversal points where the pile top bearing capacity is not lower than the preset lower limit and the total settlement and differential settlement of the composite foundation do not exceed the preset upper limit. The closed interval formed by the minimum and maximum values ​​of each key design parameter corresponding to all feasible traversal points is taken as the optimal design parameter interval for the key design parameters.

[0016] The geopolymer composite foundation parameter optimization method based on finite element analysis of the present invention has the following beneficial effects: By introducing the Cosserat micropole constitutive relation into the finite element model and setting independent micropole rotational degrees of freedom, the numerical solution of the shear band at the pile-soil interface depends on the intrinsic characteristic length of the material rather than the mesh size. This fundamentally eliminates the ill-conditioned convergence problem in the traditional Cauchy continuum model where the shear band width tends to zero with mesh refinement, significantly improving the objectivity and reliability of the calculation of the mechanical response at the pile-soil interface.

[0017] In a typical embodiment of the present invention, when the grid size is changed from the initial value... Halved to In this invention, the relative change in the width of the numerical solution of the shear band does not exceed 5%, while the traditional Cauchy continuum model shows a shear band width change of approximately 50% under the same mesh refinement, fundamentally improving the mesh dependency problem. By performing node redistribution operations based on the micropolar rotation gradient index and directional element subdivision operations based on the plastic strain index in each load increment step, adaptive tracking of the computational mesh to high gradient regions is achieved. This effectively improves local computational accuracy without significantly increasing the total number of global degrees of freedom, avoiding the waste of computational resources caused by global refinement in traditional fixed mesh methods.

[0018] By extracting the main modes of the response field based on intrinsic orthogonal decomposition and constructing an interpolation mapping from parameterized design variables to reduced-order coefficients, the computation time for a single parameter evaluation is shortened from the long time required for full-order finite element solutions to the millisecond level. In a typical embodiment, a single solution of a full-order Cosserat micro-finite element method takes approximately 2 to 4 hours, while the reduced-order fast evaluation model of this invention takes approximately 5 to 15 milliseconds, achieving a speedup of [missing information]. This scale allows for a dense evaluation of 3.2 million sampling points to be completed within hours.

[0019] By multiplying the Borgonovo sensitivity index with the Shapley interaction assignment value to construct a comprehensive sensitivity ranking value, this invention captures the influence of parameters on the shape of the response probability distribution and the fair allocation of interaction effects between parameters. This overcomes the shortcomings of traditional sensitivity methods based on variance decomposition, which cannot reflect changes in the overall shape of the distribution and ignore the fairness of parameter interaction contributions. In a typical embodiment, the key parameters identified by the comprehensive ranking of this invention differ from those ranked using the Sobol index alone in the 3rd to 5th positions. Verification examples show that the ranking results of this invention are more consistent with the engineering experience judgment of design sensitivity.

[0020] This invention can more accurately identify key design parameters and output the optimal design parameter range that satisfies bearing capacity constraints and settlement constraints, providing engineering designers with a basis for parameter selection that balances safety and robustness. Attached Figure Description

[0021] Figure 1 This is a schematic diagram of the initial grid division of the geopolymer mixing pile-soil composite foundation provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the r-type node redistribution operation provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the subdivision operation of the h-type orientation unit provided in an embodiment of the present invention; Figure 4 The diagram illustrates the principle of Borgonovo sensitivity analysis provided in this embodiment of the invention. Detailed Implementation

[0022] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0023] A method for optimizing parameters of geopolymer composite foundations based on finite element analysis includes the following steps: Step 1: Initially mesh the computational domain of the geopolymer mixing pile-soil composite foundation and record it as a unified reference mesh. Set translational degrees of freedom and micro-polar rotational degrees of freedom on the mesh nodes and adopt the Cosserat micro-polar constitutive relation within the element. Set parameterized design variables and generate several sets of parameter combinations. Perform incremental load loading on each set of parameter combinations. After each load increment step converges, perform node redistribution operation based on the rotational gradient index and perform element subdivision operation based on the plastic strain index. Extract the pile top bearing capacity, total settlement of the composite foundation and differential settlement corresponding to each set of parameter combinations, and map the response field to the unified reference mesh to form a unified dimension response snapshot vector. Step 2: Perform singular value decomposition on the snapshot matrix composed of all uniform dimension response snapshot vectors to extract the intrinsic orthogonal basis matrix. Project the uniform dimension response snapshot vector of each parameter combination onto the intrinsic orthogonal basis matrix to obtain the reduced order coefficient vector. Construct an interpolation mapping from parameterized design variables to the reduced order coefficient vector. The reduced order fast evaluation model is composed of the interpolation mapping and the intrinsic orthogonal basis matrix. Step 3: Use the reduced-order rapid evaluation model to sample and evaluate the value range of the parameterized design variables, calculate the Borgonovo sensitivity index and Shapley interaction assignment value of each parameterized design variable, multiply them to obtain the comprehensive sensitivity ranking value, select key design parameters, traverse and evaluate within the value range of the key design parameters and screen feasible solutions that meet the preset bearing capacity constraints and preset settlement constraints, and output the optimal design parameter range of the key design parameters.

[0024] In one specific implementation, the computational domain of the geopolymer mixing pile-soil composite foundation consists of three material subdomains: the geopolymer mixing pile body, the surrounding soil, and the pile tip bearing layer. Taking a soft soil foundation reinforcement project as an example, the planar extent of the computational domain extends outward from the pile center to three times the distance between adjacent piles, and the vertical extent extends downward from the ground surface to a depth of at least five times the pile diameter below the pile tip. The computational domain is initially meshed using tetrahedral elements. The element edge length in the pile body region is one-eighth to one-sixth of the pile diameter. In the surrounding soil region, the element edge length within one pile diameter distance from the pile wall is one-quarter of the pile diameter. The element edge length in the far-field soil region gradually increases to one to two times the pile diameter. The element edge length in the pile tip bearing layer region is consistent with that in the pile body region. After the initial meshing is completed, the coordinates of all nodes and the element topology connections are stored as a unified reference mesh. Subsequent snapshot vectors required for reduced-order modeling are all based on this unified reference mesh.

[0025] refer to Figure 1 , Figure 1 The initial tetrahedral meshing results of the computational domain and the partitioning methods of the three types of material subdomains are shown. Figure 1 The central gray-filled area represents the geopolymer mixing pile, the areas on either side of the pile represent the surrounding soil, and the bottom of the computational domain represents the pile tip bearing layer. The mesh density in the pile area is higher than that in the far-field area of ​​the surrounding soil. The element side length in the pile area is one-eighth to one-sixth of the pile diameter, while the element side length in the surrounding soil within a distance of one pile diameter is one-quarter of the pile diameter. The element side length in the far-field soil area gradually increases to one to two times the pile diameter. This differentiated mesh density setting is based on the following considerations: the area near the pile-soil interface is where shear deformation and stress concentration are most significant, requiring a higher mesh resolution to accurately capture the mechanical response gradient changes at the interface; while the stress and displacement distribution in the far-field soil is relatively uniform, and sufficient computational accuracy can be obtained without an excessively dense mesh. After the initial mesh generation is completed, the coordinates of all nodes and the element topology connections are recorded as a unified reference mesh. A unified reference mesh serves as the response field mapping benchmark for all parameter combinations in subsequent steps. This ensures that different parameter combinations, after undergoing varying degrees of adaptive mesh adjustments, can still generate response snapshot vectors with completely consistent dimensions, thus satisfying the requirement of eigenorthogonal decomposition for equal-dimensional column vectors in the snapshot matrix. Each mesh node is assigned 3 translational degrees of freedom and 3 independent micropole rotational degrees of freedom, for a total of 6 degrees of freedom. The 3 translational degrees of freedom correspond to displacement components in the 3 coordinate directions, and the 3 micropole rotational degrees of freedom describe the independent rotational motion of the material element at the node. Within each element, a Cosserat micropole constitutive relation is adopted, using the intrinsic characteristic length parameter... Establish the couple stress tensor With micropole rotational gradient tensor Mapping between them. Intrinsic feature length parameter. The representative dimensions characterizing the microstructure of the material are between 5 mm and 20 mm for geopolymer-mixed pile materials and between 1 mm and 10 mm for the surrounding soil. By introducing the intrinsic characteristic length parameter, the numerical solution width of the shear band at the pile-soil interface depends on the intrinsic characteristic length parameter rather than the mesh size, thus avoiding the ill-conditioned convergence problem in traditional continuum models where the shear band width tends to zero with mesh refinement. Figure 1 The initial mesh in the model forms the starting state for subsequent incremental load loading and adaptive rh-type mesh repartitioning operations.

[0026] In an optional implementation, the initial mesh can also be divided using hexahedral elements or a hexahedral-tetrahedral hybrid mesh, wherein the pile body uses a structured hexahedral mesh, the transition zone of the soil around the pile uses pyramidal elements to connect the hexahedral to tetrahedral elements, and the far-field soil uses a tetrahedral mesh.

[0027] Couple stress tensor With micropole rotational gradient tensor The mapping relationship between them is as follows: ;in For Cosserat coupled shear modulus, Cosserat flexural modulus, The components of the micropole rotation gradient tensor are defined as the micropole rotation vector. Spatial coordinates The partial derivatives of . As a coefficient appearing in the couple stress constitutive model, its physical meaning lies in: when the intrinsic characteristic length parameter As the torque increases, the resistance of the couple stress to the micropolar rotation gradient strengthens, thus increasing the shear band width; conversely, when the torque decreases... As the coefficient of friction approaches zero, the contribution of the couple stress disappears, and the Cosserat micropolar constitutive model degenerates into the traditional Cauchy constitutive model. It is this mechanism that makes the numerical solution of the shear zone at the pile-soil interface depend on the intrinsic characteristic length parameter rather than the mesh size, fundamentally solving the mesh dependency problem.

[0028] refer to Figure 2 , Figure 2 The two blue dotted lines mark the location of the pile-soil interface, with the surrounding soil and the pile itself on the left and right sides, respectively. The gray dashed lines and gray dots represent the original mesh and node positions before translation, while the black solid lines and black dots represent the updated mesh and node positions after translation. Red arrows indicate the translation direction and distance for each node. Figure 2 As can be seen, nodes closer to the pile-soil interface have all translated along the direction pointing towards the interface, resulting in a higher mesh density near the interface compared to the original mesh, while the mesh density in areas farther from the interface decreases accordingly. The driving mechanism for this node migration is the micropolar rotation gradient index. After each load increment step converges, the micropolar rotation gradient tensor at each Gaussian integration point of each element is extracted. And calculate the Frobenius norm. The maximum value of the Frobenius norm at each Gaussian integral point within the element is taken as the element's rotational gradient index. The rotational gradient index reflects the intensity of the micro-polar rotational field changes within the element; a larger value indicates that the element is more likely to be located in the core region of the pile-soil interface shear zone. The arithmetic mean of the rotational gradient indices of all elements in the entire domain is calculated. Elements whose rotational gradient indices exceed the threshold obtained by the product of the arithmetic mean and a preset r-type amplification factor are marked as r-type candidate elements. For each node of each r-type candidate element, the centroid coordinates and corresponding rotational gradient index values ​​of all adjacent elements to which the node belongs are used as discrete samples. A negative gradient estimate is performed on the rotational gradient index value. The node is translated along the negative gradient direction, i.e., the direction in which the rotational gradient index value decreases the fastest. The translation step size is the product of the shortest element side length connected to the node and a preset translation ratio factor. For example... Figure 2 As shown by the red arrow, nodes migrate towards the edge of high-gradient concentration regions, causing the originally uniformly distributed mesh nodes to re-aggregate in the shear band and its transition region. After each node translation, the Jacobian determinant of all elements associated with that node is checked. If the Jacobian determinant of any element is less than zero, the node is rolled back to its pre-translation position. A Jacobian determinant value less than zero means that the node translation caused a geometric flip of the element, and the mesh topology degenerates into an invalid state. The rollback mechanism ensures that the mesh remains valid. The advantage of r-type node redistribution is that it does not increase the total number of nodes or degrees of freedom; it only improves the mesh resolution in key regions by redistributing the spatial positions of existing nodes, thus without increasing the computational scale.

[0029] In an optional embodiment, the mapping relationship between couple stress and micropole rotation gradient can also be in an anisotropic form, that is, different intrinsic characteristic length parameter values ​​are assigned to different coordinate directions to adapt to the material anisotropy caused by the mixing process in the geopolymer mixing pile.

[0030] Cauchy stress tensor The constitutive relation between the strain metric and the strain measure can be written as: ;in Lamé's constant, The volumetric strain is the sum of the three normal strain components. For Kronecker symbol when The value is 1 if it is true and 0 otherwise. For classical shear modulus, For Cosserat coupled shear modulus, The components of the Cosserat strain tensor are defined as the difference between the displacement gradient and the Levi-Civita permutation of the micropole rotation vector, i.e. ,in For displacement components, Let the third-order Levi-Civita permutation tensor be used. At this time, the Cosserat strain tensor degenerates into the classical symmetric strain tensor, and the above constitutive relation degenerates into an isotropic linear elastic constitutive relation.

[0031] Five parametric design variables were set: pile length, pile diameter, pile spacing, geopolymer elastic modulus, and soil moisture content. In one specific implementation scenario, the pile length ranged from 6 to 15 meters, the pile diameter from 0.5 to 1.2 meters, the pile spacing from 1.5 to 3.5 meters, the geopolymer elastic modulus from 200 MPa to 800 MPa, and the soil moisture content from 25% to 55%. The ranges of these five parametric design variables were determined based on site conditions and design specifications. Latin hypercube sampling was used to generate several parameter combinations within these ranges. Latin hypercube sampling divides the range of each parametric design variable into equal intervals equal to the total number of samples, randomly selects one value from each interval, and then randomly arranges and combines the values ​​of each parametric design variable. This sampling method achieves uniform coverage of the parameter space with fewer samples compared to pure random sampling. The total number of samples is usually 10 to 30 times the number of parameterized design variables. In this embodiment, 80 sets of parameter combinations are used.

[0032] In alternative implementations, the parameterized design variables can be expanded to include more parameters, such as the mass ratio of alkali activator to powder and fly ash content in the geopolymer mix design. The sampling method can also use quasi-random, low-difference sequences such as Sobol or Hammersley sequences instead of Latin hypercube sampling.

[0033] For each parameter combination, a new geometric model is generated based on the corresponding pile length, pile diameter, and pile spacing. The pile material stiffness parameters are updated according to the corresponding geopolymer elastic modulus, and the effective stress parameters and permeability coefficient of the soil surrounding the pile are updated using the soil-water characteristic curve based on the corresponding soil moisture content. The regenerated geometric model is then meshed and incrementally loaded. The load is applied through the pile top surface as a uniformly distributed vertical pressure, increasing incrementally from zero to 1.5 to 2 times the design load. The load increment for each increment step is 5% to 10% of the total load.

[0034] Within each load increment step, the nonlinear equilibrium equations are solved using Newton-Raphson iteration. Since the total stiffness matrix of the Cosserat micropole finite element method contains coupled terms of displacement and micropole rotational degrees of freedom, each iteration requires simultaneous updates to both the displacement and micropole rotational fields. The convergence criterion is typically set at either the ratio of the residual force norm to the external load norm not exceeding 0.001, or the ratio of the displacement increment norm to the cumulative displacement norm not exceeding 0.001; convergence is determined by satisfying either of these conditions.

[0035] After the Newton-Raphson iteration converges for each load increment step, an adaptive Rh-type mesh repartition operation is performed on the current mesh, followed by node redistribution and element subdivision operations.

[0036] The process of performing node redistribution is as follows: Traverse each cell in the current mesh, and at each Gaussian integration point in each cell, calculate the micropolar rotation gradient tensor based on the numerical values ​​of the micropolar rotation degrees of freedom. Micropolar rotation gradient tensor It is a 3x3 second-order tensor containing 9 components, each representing the spatial rate of change of a certain component of the micropole rotation vector along a certain coordinate direction. Taking the Frobenius norm of the micropole rotation gradient tensor, which is the square root of the sum of the squares of the 9 components: ;in Let Frobenius norm be the infinitesimal rotational gradient tensor. For the infinitesimal rotational gradient tensor at the th... Line number The components of the column. The maximum value of the Frobenius norm is taken as the rotational gradient index for each element across all Gaussian integral points. The physical meaning of the rotational gradient index lies in reflecting the drastic spatial variation of the micropolar rotational field within the element. A larger rotational gradient index indicates that the element is more likely to be located in the core region of the pile-soil interface shear zone, requiring node migration to improve the mesh's resolution in this region.

[0037] After each node translation is completed, the Jacobian determinant value of all elements associated with that node is immediately checked. The Jacobian determinant value reflects the local volume ratio of the element's mapping from the reference configuration to the physical configuration. A positive Jacobian determinant value indicates that the node numbering order of the element is consistent with the local coordinate system direction and the element has not been flipped. A Jacobian determinant value less than zero indicates that the node translation caused geometric flipping of the element and mesh topology degradation. If the check finds that the Jacobian determinant value of any element associated with the translated node is less than zero, the node is reverted to its position before the translation, and the translation is abandoned to ensure that the mesh always maintains a valid topology.

[0038] For each h-shaped candidate element, the principal direction of the equivalent plastic strain gradient within that element is first determined. This principal direction reflects the direction of the most drastic spatial change in plastic deformation, typically perpendicular to the shear band. The h-shaped candidate element is then bisected once along the principal direction to generate two sub-elements. The bisecting surface passes through the geometric center of the h-shaped candidate element, and its normal direction is parallel to the principal direction. Unlike traditional isotropic subdivision (which uniformly divides a tetrahedron into eight sub-tetrahedrons), directional bisecting along the principal direction provides higher mesh density in the direction of the most drastic plastic strain change, without increasing the number of elements in other directions, thus achieving a better balance between computational accuracy and computational cost.

[0039] refer to Figure 3 , Figure 3 The triangular mesh in the diagram represents the local mesh region near the pile-soil interface. The area marked with a blue arrow on the right is the region of concentrated plastic strain, and the direction of the blue arrow indicates the principal direction of the equivalent plastic strain gradient. The red squares mark the nodes added due to subdivision, i.e., suspended nodes. Figure 3 As can be seen, the triangular elements near the blue arrows have been bisected along the main direction. Each subdivided element changes from one triangle to two sub-triangles. The bisecting surface passes through the geometric center of the original element, and its normal direction is parallel to the main direction. The upper left element, far from the plastic strain concentration area, retains its original size and is not subdivided. The trigger criterion for h-type element subdivision is the plastic strain index. After each load increment step converges, the equivalent plastic strain at each Gaussian integration point of each element is extracted, and the maximum value is taken as the element's plastic strain index. The arithmetic mean of the plastic strain indices of all elements in the global domain is calculated. Elements whose plastic strain index exceeds the threshold obtained by the product of the arithmetic mean and the preset h-type magnification factor are marked as h-type candidate elements. The preset h-type magnification factor is larger than the preset r-type magnification factor. The purpose is to restrict the subdivision operation to the local area where plastic deformation is most significant and avoid a sharp expansion of the total number of degrees of freedom. Directional bisectioning along the principal direction of the equivalent plastic strain gradient, rather than isotropic uniform subdivision, is used because the spatial variation of plastic strain is most drastic in the direction perpendicular to the shear band, and relatively gentler in the direction along the shear band. Directional subdivision can provide higher mesh density in the direction of most drastic variation, while not increasing the number of elements in other directions, achieving a better balance between computational accuracy and computational cost. For each suspended node generated by the subdivision, the translational degrees of freedom of the suspended node are constrained to be linear interpolations of the translational degrees of freedom of the two adjacent non-suspended nodes, and the micropolar rotational degrees of freedom of the suspended node are also constrained to be linear interpolations of the micropolar rotational degrees of freedom of the two adjacent non-suspended nodes. This constraint ensures that the displacement field and the micropolar rotational field remain continuous at the boundary between refined and unrefined elements, avoiding non-physical jumps in the stress field. Figure 3The h-type directional subdivision shown Figure 2 The r-type node redistribution shown is executed sequentially in each load increment step, and the two work together: the r-type operation improves the resolution of the rotation gradient concentration region by migrating existing nodes, and the h-type operation improves the resolution of the plastic strain concentration region by adding new nodes and elements. The two are adaptively adjusted for two different types of field variable gradients in the Cosserat micropolar constitutive model.

[0040] After the node redistribution and element subdivision operations are completed, the shape function derivative matrices of each element are recalculated based on the updated mesh, and the global stiffness matrix and load vectors are reassembled. The global stiffness matrix includes conventional stiffness sub-blocks corresponding to classical displacement degrees of freedom, couple stress stiffness sub-blocks corresponding to micropolar rotation degrees of freedom, and coupled stiffness sub-blocks between displacement and micropolar rotation. The load vectors include equivalent nodal forces corresponding to externally applied vertical uniformly distributed loads and body loads (such as self-weight). After reassembly, the process continues to the next load increment step, performing Newton-Raphson iterations until convergence, and then performing adaptive rh-type mesh repartitioning again. This cycle is repeated until the entire load loading process is completed.

[0041] After all loads are applied, three mechanical response quantities corresponding to each parameter combination are extracted from the calculation results. The pile top bearing capacity is obtained by summing the vertical reactions of all boundary nodes on the pile top surface, that is, summing the vertical constraint reactions or equivalent nodal forces of each node on the pile top surface in vertical components to obtain the total vertical force borne by the pile top surface. The total settlement of the composite foundation is the vertical displacement of the center node of the pile top surface, that is, the absolute value of the vertical displacement component of the node at the geometric center of the pile top surface under the final load state. The differential settlement is the maximum value of the difference in vertical displacement between the center node of the pile top surface and the preset observation node on the soil surface between the piles. The preset observation node is set at the midpoint of the soil surface between adjacent piles. The differential settlement reflects the degree of uneven deformation between the pile body and the soil between the piles.

[0042] Because different parameter combinations undergo varying degrees of r-type node redistribution and h-type element subdivision during incremental loading, the final mesh topology and node count for each parameter combination are not entirely identical. To enable comparison of the response results of all parameter combinations in the same vector space and subsequent order reduction analysis, it is necessary to uniformly map the global nodal displacement field and micropole rotation field obtained by each parameter combination on its respective adaptive mesh to the nodes of a unified reference mesh. The mapping method employs inverse distance interpolation. For each node on the unified reference mesh, several nodes closest to that node are found in the adaptive mesh (in this embodiment, the 8 closest nodes are selected). The reciprocal of the distance from these nodes to the nodes of the unified reference mesh is used as the interpolation coefficient. The displacement and micropole rotation values ​​of these nodes are then weighted and averaged according to the interpolation coefficient to obtain the mapped value on the nodes of the unified reference mesh. In optional embodiments, radial basis function interpolation or projection mapping based on element shape functions can also be used for interpolation. After mapping, each set of parameters yields a vector that is completely consistent with the number of nodes and degrees of freedom of the unified reference mesh. The three displacement components and three micropole rotation components of each node are arranged in the order of node number to form the unified dimension response snapshot vector corresponding to the set of parameters.

[0043] It should be noted that the role of the unified reference grid is solely to provide snapshot vectors with consistent dimensions for the intrinsic orthogonal decomposition, thereby satisfying the requirement of equal-dimensional column vectors in the snapshot matrix for subsequent singular value decomposition. The spatial resolution of the unified reference grid is configured based on its ability to faithfully represent the three engineering scalar responses extracted in this step: pile top bearing capacity, total settlement of the composite foundation, and differential settlement. Therefore, it can be determined once during the initial mesh generation stage. The three engineering scalar responses are essentially spatial integrals or displacement components of specific nodes, and the unified reference grid is not required to reproduce the high-resolution details achieved by the adaptive grid in the local shear zone point by point. The capture of the local high-resolution details of the shear zone is completed by the respective rh-type adaptive grids of each parameter combination in the original solution stage, and objectivity is ensured by the physical feature length determined by the Cosserat micropolar constitutive relation rather than the grid size. The intrinsic orthogonal basis matrix extracted from the snapshot matrix in the subsequent singular value decomposition in step 2 reflects the low-rank dominant mode of the response field among different parameter combinations, which has an averaging effect on the grid discretization error, further reducing the impact of the inverse distance interpolation on the final accuracy of the reduced-order model.

[0044] After step 1 is completed, each parameter combination corresponds to a unified-dimensional response snapshot vector, and all parameter combinations generate several vectors with identical dimensions. Taking 80 parameter combinations in this embodiment as an example, the unified reference mesh contains approximately 12,000 nodes, each node has 6 degrees of freedom (3 displacement components and 3 micropole rotational components), so the dimension of each unified-dimensional response snapshot vector is 72,000. Directly performing a complete Cosserat micropole finite element solution on each of the 80 parameter combinations would consume a large amount of computational resources, and in the subsequent sensitivity analysis and parameter optimization stages, thousands or even tens of thousands of parameter combinations need to be evaluated. To resolve this contradiction, a reduced-order fast evaluation model is constructed based on the intrinsic orthogonal decomposition method, transforming the prediction of the high-dimensional response field into an interpolation problem in the low-dimensional coefficient space, reducing the computation time of a single evaluation from several hours to the millisecond level.

[0045] Arrange the uniform-dimensional response snapshot vectors corresponding to all 80 parameter combinations column-wise to form a snapshot matrix. .matrix The number of rows is the total number of all degrees of freedom on the unified reference grid. (In this embodiment) The number of columns is the total number of parameter combinations. (In this embodiment) ),Right now It is The matrix, where Represents the total number of degrees of freedom. This represents the total number of snapshots. Each column of the snapshot matrix represents the complete response field of the composite foundation under the final load state under a specific set of parameter combinations, while each row represents the variation of a specific degree of freedom under different parameter combinations.

[0046] For snapshot matrix Perform singular value decomposition. Singular value decomposition decomposes any matrix into the product of three matrices: ;in for The left singular vector matrix, Each column is called a left singular vector, representing a fundamental mode of the response field in the space of degrees of freedom; for A diagonal matrix, where the elements on the diagonal are singular values. Arranged in descending order, each singular value reflects the energy percentage of the corresponding basic mode in all snapshots; for The transpose of the right singular vector matrix, Each column describes the coefficient distribution of the corresponding basic pattern across different parameter combinations. Superscript This represents the matrix transpose operation.

[0047] Starting from the first order, the squares of the singular values ​​are accumulated sequentially. The cumulative energy ratio is defined. For the front The proportion of the sum of squares of singular values ​​of order 1 to the total sum of squares of all singular values: ;in To retain the previous The cumulative energy ratio in the first-order mode, For the first Singular values ​​of order, This represents the current accumulated order. The total number of singular values ​​is the total number of snapshots. From Start increasing step by step The value, until The preset energy retention threshold is reached for the first time. In this embodiment, the preset energy retention threshold is set to 0.9999, meaning that 99.99% of the energy is retained. In actual calculations, for the geopolymer composite foundation problem in this embodiment, the first 8 to 15 orders of left singular vectors can usually reach this threshold. The specific order depends on the complexity of the influence of the parameterized design variables on the spatial distribution pattern of the response field.

[0048] In an optional implementation, the preset energy retention threshold can be between 0.999 and 0.99999. A higher threshold results in more retained orders and higher reconstruction accuracy, but a weaker reduction effect; a lower threshold results in fewer retained orders and faster computation speed, but may result in some loss of detail accuracy.

[0049] Projecting the uniform-dimensional response snapshot vector corresponding to each parameter combination onto the intrinsic orthogonal basis matrix yields the reduced-order coefficient vector. The projection operation is calculated as follows: ;in For the first The reduced-order coefficient vector corresponding to the combination of parameters. It is the transpose of the eigenorthogonal basis matrix. For the first The unified-dimensional response snapshot vector corresponding to the combination of parameters. Since the column vectors of the intrinsic orthogonal basis matrix are mutually orthogonal and have a magnitude of 1, the above projection operation essentially extracts the components of the high-dimensional response vector along the directions of each orthogonal basis. Reduced-order coefficient vector. The dimension is (In this embodiment, it is 12), the physical meaning of each component is the first... The projection amplitude of the response field of the parameter combination onto the corresponding order eigenorthogonal basis mode.

[0050] Radial basis functions By distance For the independent variable, a quadratic radial basis function is selected in this embodiment, and its expression is: ;in Let be the Euclidean distance between two parameter vectors. These are shape parameters that control the smoothness of the radial basis functions. Typically, the distance between known sample points in the parameter space is taken as 0.5 to 2 times the average distance. In this implementation, the optimal value is determined by leave-one-out cross-validation: one sample point is removed sequentially, and the remaining 79 sample points are used to construct an interpolation mapping to predict the reduced-order coefficient vector of the removed sample points. The prediction error is calculated, and several candidate points are traversed. After selecting the value, choose the one that minimizes the average prediction error. The value is used as the final shape parameter.

[0051] In an optional implementation, the radial basis function can also be a Gaussian radial basis function. or inverse quadratic radial basis functions .

[0052] In an alternative implementation, the interpolation mapping can also use a Kriging surrogate model or support vector regression instead of radial basis function interpolation. The Kriging surrogate model can simultaneously provide the predicted value and the estimate of the prediction uncertainty, making it suitable for scenarios requiring adaptive point supplementation.

[0053] Radial basis function interpolation mapping and intrinsic orthogonal basis matrices together constitute the order-reduced fast evaluation model. For any new combination of parameterized design variables, the evaluation process of the order-reduced fast evaluation model consists of the following steps: First, the values ​​of the five new parameterized design variables are combined to form a parameter vector. The predicted order reduction coefficient vector is obtained by radial basis function interpolation mapping. Then, the predicted reduced-order coefficient vectors are linearly combined using the eigenorthogonal basis matrix to reduce them to the global response field on the unified reference grid. ;in The reconstructed global response field vector, It is an eigenorthogonal basis matrix. The vector for predicting the order reduction coefficients is used. The reconstructed global response field vector corresponds one-to-one with the nodes on the unified reference grid, containing three displacement components and three micro-polar rotation components for each node. Finally, the predicted values ​​are obtained from the reconstructed response field using the same extraction method as in step 1 for pile top bearing capacity, total settlement of the composite foundation, and differential settlement: the predicted value of pile top bearing capacity is obtained by summing the vertical components of all boundary nodes on the pile top surface; the predicted value of total settlement of the composite foundation is obtained by reading the vertical displacement components of the center node on the pile top surface; and the predicted value of differential settlement is obtained by calculating the maximum value of the vertical displacement difference between the center node on the pile top surface and the preset observation node on the soil surface between the piles.

[0054] In an alternative implementation, step 2 can also directly establish a regression mapping from parametric design variables to pile top bearing capacity, total settlement of composite foundation, and differential settlement in the reduced-order coefficient vector space, skipping the reconstruction step of the global response field and further accelerating the evaluation process. However, this approach loses global response field information and is not suitable for application scenarios that require visualization of global stress or displacement distribution.

[0055] A reduced-order rapid evaluation model is used to sample and evaluate the value range of the parameterized design variables. Uniform sampling is performed within the value range of the five parameterized design variables at a preset sampling density. The preset sampling density refers to the number of sampling points evenly distributed within the value range of each parameterized design variable. In this embodiment, 20 evenly spaced sampling points are arranged on each parameterized design variable, and the total number of sampling points for the full combination of the five parameterized design variables is [number missing]. The total evaluation time for 3.2 million sampling points is still within an acceptable range, given that the single evaluation time of the reduced-order rapid evaluation model is on the order of milliseconds. The predicted values ​​of pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each sampling point are obtained using the reduced-order rapid evaluation model.

[0056] In an optional implementation, if the number of parameterized design variables is too large, resulting in an excessively large total number of sampling points for the entire combination, a quasi-random sequence can be used for sampling to reduce the total number of sampling points, or the number of sampling points on each parameterized design variable can be reduced to 10 to 15.

[0057] refer to Figure 4 , Figure 4 It contains 6 sub-images. Figure 4 (a) to Figure 4 (e) The effects of five parametric design variables—pile length, pile diameter, pile spacing, geopolymer elastic modulus, and soil moisture content—on the conditional distribution shift of pile top bearing capacity are shown. Figure 4 (f) The ranking results of the Borgonovo sensitivity index for the five parameterized design variables are summarized. Figure 4 (a) to Figure 4 In each subplot of (e), the horizontal axis represents the pile top bearing capacity in kN, and the vertical axis represents the probability density in kN. / kN. The thick black solid line in each sub-figure represents the unconditional probability density curve of the pile top bearing capacity, which is the overall distribution curve obtained by fitting the predicted pile top bearing capacity values ​​of all sampling points to a Gaussian kernel density estimate without limiting the values ​​of any parametric design variables. The four colored curves represent the conditional probability density curves of the predicted pile top bearing capacity values ​​within each interval after dividing the value range of the corresponding parametric design variables into several intervals according to a preset number of segments. The blue dashed line, red dotted line, green dotted line, and purple double dotted line correspond to intervals 1 to 4, respectively. Observation Figure 4(a) It can be found that the four conditional probability density curves of pile length have a significant horizontal shift compared to the unconditional probability density curve. The peak positions of the four curves are significantly dispersed along the horizontal axis, and the shape of each conditional distribution is more concentrated and sharper than that of the unconditional distribution. This indicates that the change in the value of pile length has a strong influence on the distribution of the bearing capacity at the top of the pile. Figure 4 (a) The Borgonovo sensitivity index value, marked in the upper right corner, indicates the pile length. This is the highest among the five parametric design variables. In comparison, Figure 4 (b) shows the pile diameter and Figure 4 (e) The conditional probability density curve of soil moisture content almost coincides with the unconditional probability density curve. The peak positions of the four colored curves are nearly identical to the black curve, with minimal offset. The corresponding Borgonovo sensitivity index values ​​are respectively... and This indicates that the pile diameter and soil moisture content have a relatively weak impact on the distribution of bearing capacity at the pile top. Figure 4 (c) shows the pile spacing and Figure 4 (d) shows that the elastic modulus of the geopolymer lies between the two, with a moderate degree of shift in its conditional distribution. The corresponding Borgonovo sensitivity index values ​​are respectively... and The Borgonovo sensitivity index measures the intensity of parameter influence by quantifying the deviation between conditional and unconditional distributions. Specifically, it is calculated by taking the smaller of the conditional probability density curve value and the unconditional probability density curve value at each point along the range of pile top bearing capacity within each interval, performing numerical integration, subtracting 1 from the integral result to obtain the local offset of the interval, and then taking the arithmetic mean of the local offsets of all intervals. Figure 4 (f) Arrange the five parametric design variables in descending order of their Borgonovo sensitivity index values ​​using a horizontal bar chart. The black bars correspond to the top three parametric design variables: pile length, pile spacing, and geopolymer elastic modulus. The gray bars correspond to the bottom two: pile diameter and soil moisture content. The red numbers to the right of the bars represent the corresponding Borgonovo sensitivity index values. The top three parametric design variables will be selected as key design parameters for determining the optimal design parameter range.

[0058] In an optional implementation, the preset number of segments can be between 5 and 20. A higher number of segments provides a more detailed characterization of the conditional distribution, but a reduction in the number of data points within each segment may decrease the accuracy of kernel density estimation. The kernel density estimation method can also be replaced by histogram estimation or adaptive kernel density estimation.

[0059] The process of calculating the Shapley interaction allocation value is as follows. The Shapley interaction allocation value originates from the concept of Shapley values ​​in cooperative game theory and is used to measure the contribution of each parameterized design variable to the variability of the mechanical response in its interaction with other parameterized design variables. The Borgonovo sensitivity index measures the influence of a single parameterized design variable on the shape of the response distribution, but does not consider the interaction effects between parameters. For example, the effects of pile diameter and pile spacing on bearing capacity may be strongly coupled—the sensitivity of pile spacing may change as the pile diameter increases—the Borgonovo index cannot distinguish such interaction contributions. The Shapley interaction allocation value, however, fairly allocates the interaction effects to each parameter by traversing all possible subsets of parameters.

[0060] Consider the five parameterized design variables as five participants in a cooperative game. Define the characteristic function: for any subset of parameters... ( (for any subset of the five parametric design variables), characteristic function values For a subset The change in variance of the predicted mechanical response when the parametric design variables are fixed to specific values.

[0061] For each parameter subset, the characteristic function value is calculated as follows: The parameterized design variables within the parameter subset are fixed to their discrete sampled values ​​in a preset sample set. The parameterized design variables outside the parameter subset are iterated through all discrete sampled values ​​in the preset sample set. The preset sample set consists of all sampling points generated by the aforementioned uniform sampling. After the parameterized design variables within the parameter subset are fixed to a specific set of values, the predicted mechanical response values ​​generated by the free variation of the parameterized design variables outside the parameter subset constitute a conditional sample set. The variance of the conditional sample set is calculated, and then the average variance of all possible fixed values ​​within the parameter subset is taken. The final average conditional variance is the characteristic function value of the parameter subset. Its physical meaning is: when the parameter subset When the values ​​of the parametric design variables are known, the remaining uncertainty of the mechanical response is considered.

[0062] The first After adding the parameterized design variables to the parameter subset, an extended subset is formed. The characteristic function values ​​of the extended subset are calculated in the same way. The difference between the eigenfunction values ​​of the expanded subset and the eigenfunction values ​​of the original parameter subset is taken as the first... A parameterized design variable in a subset of parameters Marginal contribution below ,in Indicates the first Add a parameterized design variable to the parameter subset. The change in conditional variance resulting from this. The physical meaning of the marginal contribution is: after considering a subset of parameters... After considering the impact of all parametric design variables, the third... The variance increment that can be further explained by each parameterized design variable.

[0063] The corresponding Shapley combination weights are determined based on the size of the parameter subset. Let the parameter subset be... Include There are 1 parameterized design variable, and the total number of all participants is 1. Then the parameter subset The corresponding Shapley portfolio weights are: ;in For scale The Shapley combination weights corresponding to the subset of parameters. for factorial, for factorial, for The factorial of 5 is equal to 120. This refers to the total number of participants, i.e., the number of parameterized design variables. The source of this weighting formula is: in all... In the order of participants, the parameter subset Just ranked in The number of permutations before each parameterized design variable is: The Shapley value requires that each permutation contributes equally, therefore the weighting is based on the permutation ratio. For example, .

[0064] The first The marginal contribution of each parameterized design variable under all 16 parameter subsets is multiplied by its corresponding Shapley portfolio weight, and then summed to obtain the result. Shapley interaction assignment values ​​for 5 parameterized design variables. The above process is performed sequentially on the 5 parameterized design variables to obtain 5 Shapley interaction assignment values.

[0065] After calculating the Borgonovo sensitivity index and Shapley interaction assignment value for the five parameterized design variables, the Borgonovo sensitivity index and Shapley interaction assignment value for each parameterized design variable are multiplied to obtain the overall sensitivity ranking value. The logic behind this multiplication and fusion is as follows: the Borgonovo sensitivity index captures the strength of a parameter's influence on the overall shape of the response distribution, while the Shapley interaction assignment value fairly measures the parameter's contribution to variance in an interactive environment. Using the Borgonovo index alone may overestimate a parameter that has a large impact on the distribution shape but a small contribution to the total variance; using the Shapley value alone may overlook a parameter with a moderate contribution to variance but a significant impact on the distribution shift. After multiplying the two, only parameters that perform outstandingly in both distribution influence and variance contribution can obtain a higher overall sensitivity ranking value, thus more accurately identifying the truly critical design parameters.

[0066] After calculating the comprehensive sensitivity ranking values ​​for the three mechanical response quantities (pile top bearing capacity, total settlement of composite foundation, and differential settlement), the sum of the comprehensive sensitivity ranking values ​​for each mechanical response quantity is taken as the final ranking criterion. The five parametric design variables are arranged from largest to smallest according to the comprehensive sensitivity ranking values, and the top three parametric design variables are selected as key design parameters.

[0067] In an optional implementation, the number of key design parameters can also be adaptively determined based on the distribution characteristics of the comprehensive sensitivity ranking values. For example, the top few parameterized design variables whose cumulative comprehensive sensitivity ranking values ​​account for more than 80% of the total ranking values ​​can be selected as key design parameters.

[0068] A three-dimensional mesh is formed by uniformly spacing sampling points within the value ranges of three key design parameters, while the other two non-key design parameters are fixed at the median of their respective ranges. Uniformly spacing meshing refers to evenly distributing sampling points within the value ranges of each of the three key design parameters to form a three-dimensional mesh. In this embodiment, 50 equally spaced points are distributed on each key design parameter, resulting in a three-dimensional mesh containing... A gridded traversal of points.

[0069] The pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each gridded traversal point are evaluated using a reduced-order rapid evaluation model. Gridded traversal points that simultaneously meet the following conditions are marked as feasible traversal points: the pile top bearing capacity is not lower than the preset lower limit, the total settlement of the composite foundation does not exceed the preset upper limit, and the differential settlement does not exceed the preset upper limit. The preset lower limit of bearing capacity is determined based on the design load and safety factor; in this embodiment, it is taken as 1.2 times the design load. The preset upper limit of settlement includes an upper limit of 200 mm for the total settlement of the composite foundation and an upper limit of 30 mm for the differential settlement; the specific values ​​are determined according to engineering design specifications and building deformation control levels.

[0070] In optional implementations, the preset bearing capacity constraints and preset settlement constraints can also be set in different grades according to different engineering levels, or additional constraints such as minimizing material usage can be introduced.

[0071] The closed interval formed by the minimum and maximum values ​​of each critical design parameter corresponding to all feasible traversal points is output as the optimal design parameter interval for the critical design parameter. Taking pile length as an example, if the minimum pile length among all feasible traversal points is 8.4 meters and the maximum is 13.2 meters, then the optimal design parameter interval for the pile length is the closed interval from 8.4 meters to 13.2 meters. Similarly, the optimal design parameter intervals for each of the three critical design parameters are determined. Engineering designers can then determine the specific design parameter values ​​within the optimal design parameter interval based on factors such as construction conditions, material supply, and cost control.

[0072] The reason for outputting an optimal design parameter range rather than a single optimal point is that the actual construction of geopolymer composite foundations inevitably involves parameter fluctuations, and the actual values ​​of pile length, pile diameter, and pile spacing during on-site construction will deviate from the design values. Providing a range rather than a single point allows designers to flexibly adjust within the range, while ensuring that the bearing capacity and settlement meet the design requirements under any combination of values ​​within the range, thereby improving the engineering robustness of the design scheme. The method described in this invention was used to conduct numerical verification on a geopolymer mixing pile composite foundation in a typical soft soil site. The computational domain was taken as a horizontal range extending outward from the pile center to three times the pile spacing and a vertical range five times the pile diameter below the pile tip. Tetrahedral elements were used for initial mesh generation, with approximately 12,000 nodes and 50,000 elements in the initial mesh, each node having 6 degrees of freedom, for a total of approximately 72,000 degrees of freedom. Latin hypercube sampling was used to generate 80 parameter combinations as snapshot samples within the range of five parametric design variables. Incremental load was applied to each parameter combination, and rh-type adaptive mesh repartitioning was performed after convergence of each load increment step. Singular value decomposition was performed on the uniform-dimensional response snapshot vectors of all 80 parameter combinations. The cumulative energy ratio of the first 12 left singular vectors reached over 99.99%, forming an intrinsic orthogonal basis matrix. A dense sampling network of 3,200,000 sampling points was constructed within the range of five parametric design variables, with 20 equally spaced sampling points per variable. A single evaluation was performed on each sampling point using a reduced-order fast evaluation model. Experimental data showed that the full-order Cosserat micro-finite element solution for a single parameter combination took approximately 2 to 4 hours, while the reduced-order fast evaluation model of this invention took approximately 5 to 15 milliseconds for a single evaluation, achieving an overall speedup of [missing data]. The evaluation of all 3,200,000 sampling points can be completed within hours. The key design parameters determined by the product ranking of the Borgonovo sensitivity index and the Shapley interaction assignment value are pile length, pile spacing, and geopolymer elastic modulus; the optimal design parameter range is pile length from 8.4 to 13.2 meters, pile spacing from 1.8 to 2.6 meters, and geopolymer elastic modulus from 400 to 700 MPa; any combination of values ​​within the range satisfies the engineering constraints that the pile top bearing capacity is not less than 1.2 times the design load, the total settlement does not exceed 200 mm, and the differential settlement does not exceed 30 mm.

[0073] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and systems in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.

Claims

1. A method for optimizing parameters of geopolymer composite foundations based on finite element analysis, characterized in that, Includes the following steps: Step 1: Initially mesh the computational domain of the geopolymer mixing pile-soil composite foundation and record it as a unified reference mesh. Set translational degrees of freedom and micro-polar rotational degrees of freedom on the mesh nodes and adopt the Cosserat micro-polar constitutive relation within the element. Set parameterized design variables and generate several sets of parameter combinations. Perform incremental load loading on each set of parameter combinations. After each load increment step converges, perform node redistribution operation based on the rotational gradient index and perform element subdivision operation based on the plastic strain index. Extract the pile top bearing capacity, total settlement of the composite foundation and differential settlement corresponding to each set of parameter combinations, and map the response field to the unified reference mesh to form a unified dimension response snapshot vector. Step 2: Perform singular value decomposition on the snapshot matrix composed of all uniform dimension response snapshot vectors to extract the intrinsic orthogonal basis matrix. Project the uniform dimension response snapshot vector of each parameter combination onto the intrinsic orthogonal basis matrix to obtain the reduced order coefficient vector. Construct an interpolation mapping from parameterized design variables to the reduced order coefficient vector. The reduced order fast evaluation model is composed of the interpolation mapping and the intrinsic orthogonal basis matrix. Step 3: Use the reduced-order rapid evaluation model to sample and evaluate the value range of the parameterized design variables, calculate the Borgonovo sensitivity index and Shapley interaction assignment value of each parameterized design variable, multiply them to obtain the comprehensive sensitivity ranking value, select key design parameters, traverse and evaluate within the value range of the key design parameters and screen feasible solutions that meet the preset bearing capacity constraints and preset settlement constraints, and output the optimal design parameter range of the key design parameters.

2. The method according to claim 1, characterized in that, In step 1, the initial mesh is divided into an initial tetrahedral mesh, which divides the geopolymer mixing pile, the surrounding soil, and the bearing layer at the pile tip into three material subdomains. Each mesh node is given three translational degrees of freedom and three independent micropolar rotational degrees of freedom. The parametric design variables are pile length, pile diameter, pile spacing, geopolymer elastic modulus, and soil moisture content, totaling five. Several sets of parameter combinations are generated by Latin hypercube sampling within the value range of the five parametric design variables. Each load increment step is solved until convergence through Newton-Raphson iteration.

3. The method according to claim 1 or 2, characterized in that, In step 1, the Cosserat micropole constitutive relation establishes a mapping between the couple stress tensor and the micropole rotational gradient tensor by introducing an intrinsic characteristic length parameter, so that the numerical solution of the shear zone at the pile-soil interface depends on the intrinsic characteristic length parameter rather than the mesh size.

4. The method according to claim 1, characterized in that, In step 1, the node redistribution operation includes: traversing each cell in the current mesh, extracting the micro-polar rotation gradient tensor at each Gaussian integration point of the cell, and taking the maximum value of the Frobenius norm as the rotation gradient index of the cell; calculating the arithmetic mean of the rotation gradient indices of all cells in the global domain; marking cells whose rotation gradient indices exceed the threshold obtained by the product of the arithmetic mean and the preset r-type amplification factor as r-type candidate cells; for each node of each r-type candidate cell, translating the node along the direction of the fastest decrease in the rotation gradient index value, with the translation step size being the product of the shortest cell side length connected to the node and the preset translation ratio factor; after each node translation is completed, checking the Jacobian determinant value of all cells associated with the node, and if the Jacobian determinant value of any cell is less than zero, the node is rolled back to its position before translation.

5. The method according to claim 1, characterized in that, In step 1, the element subdivision operation includes: traversing each element in the current mesh, extracting the equivalent plastic strain at each Gaussian integration point of the element and taking the maximum value as the element's plastic strain index; calculating the arithmetic mean of the plastic strain indices of all elements in the global domain; marking elements whose plastic strain indices exceed the threshold obtained by the product of the arithmetic mean and the preset h-type magnification factor as h-type candidate elements; for each h-type candidate element, performing one mid-surface bisection along the principal direction of the equivalent plastic strain gradient within the h-type candidate element to generate two sub-elements; for each suspended node generated by the subdivision, constraining the translational degrees of freedom of the suspended node to the linear interpolation of the translational degrees of freedom of two adjacent non-suspended nodes, and simultaneously constraining the micro-polar rotational degrees of freedom of the suspended node to the linear interpolation of the micro-polar rotational degrees of freedom of two adjacent non-suspended nodes.

6. The method according to claim 1, characterized in that, In step 1, the unified reference grid is the node coordinates and topological relationships of the initial grid; the pile top bearing capacity is obtained by summing the vertical reactions of all boundary nodes on the pile top surface; the total settlement of the composite foundation is the vertical displacement of the center node on the pile top surface; the differential settlement is the maximum value of the vertical displacement difference between the center node on the pile top surface and the preset observation node on the soil surface between the piles; the response field is mapped to the unified reference grid by mapping the global node displacement field and micro-polar rotation field obtained by combining each set of parameters on the adaptive grid to the nodes of the unified reference grid through inverse distance interpolation.

7. The method according to claim 1, characterized in that, In step 2, the singular value decomposition obtains the left singular vector sequence and the corresponding singular value sequence. Starting from the first order, the squares of the singular values ​​are accumulated step by step until the proportion of the accumulated result to the total sum of the squares of all singular values ​​reaches the preset energy retention threshold for the first time. The corresponding first few left singular vectors are then grouped into an eigenorthogonal basis matrix.

8. The method according to claim 1 or 7, characterized in that, In step 2, the interpolation mapping is a radial basis function interpolation mapping, which is constructed with the parameterized design variable values ​​corresponding to all parameter combinations as input nodes and the corresponding reduced-order coefficient vectors as output nodes. For a new combination of parameterized design variables, the predicted reduced-order coefficient vector is obtained through radial basis function interpolation mapping, and then the predicted reduced-order coefficient vector is linearly combined and restored to the global response field on the unified reference grid through the eigenorthogonal basis matrix. The predicted values ​​are obtained in the same way as the pile top bearing capacity, total settlement of composite foundation and differential settlement in step 1.

9. The method according to claim 1, characterized in that, In step 3, the sampling evaluation involves uniformly sampling within the range of values ​​of the five parametric design variables at a preset sampling density, and obtaining the predicted values ​​of the pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each sampling point through a reduced-order rapid evaluation model. The calculation of the Borgonovo sensitivity index includes: for the i-th parameterized design variable, the range of values ​​of the i-th parameterized design variable is divided into several intervals according to a preset number of segments; Within each interval, the predicted values ​​of the mechanical response corresponding to all sampling points falling within the interval are collected, and the conditional probability density curve of the mechanical response within the interval is fitted using the Gaussian kernel density estimation method. Simultaneously, the unconditional probability density curve is fitted to the predicted values ​​of the mechanical response corresponding to all sampling points using the same method. Within each interval, the smaller of the conditional probability density curve value and the unconditional probability density curve value is taken point-by-point along the value range of the mechanical response, and numerical integration is performed. The integral result is subtracted from 1 to obtain the local offset of the interval. The arithmetic mean of the local offsets of all intervals is taken to obtain the Borgonovo sensitivity index of the i-th parameterized design variable.

10. The method according to claim 9, characterized in that, In step 3, the calculation of the Shapley interaction allocation value includes: treating all parameterized design variables as participants in a cooperative game; for the i-th parameterized design variable, enumerating all parameter subsets that do not contain the i-th parameterized design variable; for each parameter subset, fixing the parameterized design variables within the parameter subset to their respective discrete sampled values ​​in a preset sample set, and iterating through all discrete sampled values ​​of the parameterized design variables outside the parameter subset in the preset sample set, calculating the variance of the predicted mechanical response value as the characteristic function value of the parameter subset; after adding the i-th parameterized design variable to the parameter subset, calculating the characteristic function value of the expanded subset in the same way, and using the difference between the characteristic function value of the expanded subset and the characteristic function value of the parameter subset as the marginal contribution under the parameter subset; determining the corresponding Shapley combination weight according to the size of the parameter subset, and all parameter subsets are assigned a weight based on the variance of the predicted mechanical response value. The sum of the marginal contributions after weighting the Shapley combination yields the Shapley interaction assignment value of the i-th parameterized design variable. The key design parameters are the top three parameterized design variables selected after sorting all parameterized design variables by comprehensive sensitivity from largest to smallest. The optimal design parameter interval is determined by performing an equally spaced grid traversal within the value range of the three key design parameters, evaluating the pile top bearing capacity, total settlement of the composite foundation, and differential settlement at each grid traversal point using a reduced-order rapid evaluation model, and selecting all feasible traversal points where the pile top bearing capacity is not lower than the preset lower limit and the total settlement and differential settlement of the composite foundation do not exceed the preset upper limit. The closed interval formed by the minimum and maximum values ​​of each key design parameter corresponding to all feasible traversal points is taken as the optimal design parameter interval for the key design parameters.

Images