## 1 Introduction  
Numerous practical projects widely use two-layered strip foundations. To reasonably determine the ultimate bearing capacity of two-layered strip footing, scholars have conducted many studies on the solution methods, such as the diffusion angle method [1], the Hanson weighted average method [2], shear failure theory [3, 4], and the limit analysis method [5, 6, 7]. These methods only take into account the velocity model and soil energy dissipation when applying the upper-bound limit analysis theorem. With this approach, the ultimate load can be solved immediately without the complex elastic-plastic body being deformed. The upper-bound limit analysis theory has nonetheless been adopted by numerous researchers due to its simplicity even if it frequently somewhat overestimates the ultimate load. Two-layered soils differ mechanically from homogeneous soil foundations in that they have more complex failure modes and bearing processes. The ultimate carrying capability of a layered soil foundation must therefore be determined differently. Determining the ultimate bearing capacity of a two-layered foundation is crucial in both theory and practice. Scholars have relatively limited research on the critical sliding surface of double-layer foundations. Therefore, studying the ultimate bearing capacity of double-layer foundations has great practical significance. How to determine the critical sliding surface of double-layer foundations and calculate the corresponding ultimate bearing capacity of foundations to reasonably evaluate the safety status of foundations is still a subject that needs to be continuously improved. With the development of numerical methods, the finite element method is more and more used to study the ultimate bearing capacity of foundation. [PERSON] [8] used FEM to calculate the numerical results of the undrained vertical bearing capacity of a rough ring foundation on two layers of clay, and compared the small displacement finite element prediction with existing experience, analysis, and numerical solutions to determine the coupling effects of the dimensionless parameters. [PERSON] [9] proposed a smooth spread method based on enrichment elements for numerical modeling of regions including weak and strong discontinuities. The triangular background grid is used to discretize the region. [PERSON] [10] mainly studied the influence of shear strength anisotropy on the bearing capacity of surface foundation near the slope. The lower limit analysis method is used with finite element and linear programming techniques to predict the limit pressure. Experimental methods can also be used to analyze the mechanical properties of layered soil. [PERSON] [11] analyzed the mechanical properties of backfill soil when the bedding plane is inclined, and as the dip angle of the interlayer increases, the crack distribution of the fine filler becomes more concentrated. Numerous studies have been conducted recently on various optimization techniques to locate the crucial slip surface and ascertain the ultimate bearing capacity of the foundation. Based on the limit analysis, [PERSON] et al. [12] created an optimization model of the two-layered clay foundation through the rigid block discrete system, then used the dichotomy method to find the critical slip surface and determine the ultimate bearing capacity. By using an updated Hoek-Brown failure criterion, [PERSON] [13] examined the bearing capacity of a shallow rigid foundation on rock mass under horizontal seepage. [PERSON] [14] used the upper-bound theorem in conjunction with the spatial discretization technique to calculate the ultimate bearing capacity of the foundation in nonhomogeneous and anisotropic clays. For the reinforced foundation under the strip footing, [PERSON] [15,16] combined the critical sliding field method with the limit equilibrium slice method to locate the critical slip surface and compute the related ultimate bearing capacity.