A numerical model for settlement analysis of circular plates on multilayered soil

A numerical model for settlement analysis of circular plates on multilayered soil Enrique Justo * , Isabel Gonz´ alez-de-Le´ on, Manuel V´ azquez-Boza Department of Building Structures and Ground Engineering, University of Seville, Avda Reina Mercedes 2, Seville 41012, Spain ARTICLE INFO Keywords: Circular plate Elastic foundation Multilayered soil Settlement Soil-structure interaction Elastic continuum method ABSTRACT This paper presents a numerical model for calculating settlements and contact stresses of a circular plate resting on an elastic subgrade. The method is an extension of the elastic continuum method developed by Poulos and Davis for piles. Soil settlements are calculated with Mindlin’s equations. Plate settlements are calculated through a finite difference approximation of Kirchhoff’s equations for thin plate bending. The method, originally devised for homogenous soils, has been extended for multilayered soils using the Steinbrenner approximation. Model validation was performed by comparing results with a finite element solution and with previously published methods. The results prove that the method provides a very good approximation for homogenous soils and also for multilayered soils in which soil stiffness increases with depth, while for layered soils with stiffness decreasing with depth the Steinbrenner approximation was found not to be sufficiently accurate. Compared to alternative numerical methods, such as those using variational calculus, the proposed method has the advantage of its greater simplicity. 1. Introduction The problem of calculating the settlements of a circular plate resting on an elastic foundation has important applications in structural and geotechnical engineering. Relevant examples include the foundation of oil and water tanks. It is a soil-structure interaction problem, in which the soil and plate behaviours are coupled, so that each one influences the other. 1.1. Plate resting on homogeneous elastic soil A number of methods have been applied to the settlement analysis of a plate resting on a homogenous elastic soil. The simplest representation of the soil response is the Winkler model [1], in which the soil medium is substituted by a number of independent linear springs. The Winkler model has the well-known shortcoming of neglecting the shear interaction between adjacent soil elements. Thus, a loaded region only produces deformation in the soil lying directly below. Subsequently, improved Winkler models introducing various kinds of interaction among the springs have been presented, such as the two-parameter Pasternak foundation [2], but they require new parameters that are difficult to obtain in practical applications. Some strictly analytical solutions have been published for the problem of a plate resting on a Winkler soil, but in all cases the results are provided in terms of special functions or singular integrals, whose evaluation is not straightforward. Pavlou and Vlachakis [3] and Koreneva and Grossman [4] presented analytical solutions for the problem of a plate resting on a Winkler foundation based on Hankel integral transforms and Bessel functions. Foyouzat and Mofid [5] and Krutii et al. [6] published solutions for a plate on a Winkler soil with a variable modulus, obtained in the form of infinite power series. The problem of a rectangular plate resting on a Pasternak foundation was addressed analytically by Shi et al. [7] and a solution in terms of infinite Fourier series was provided. Peng et al. [8] calculated the settlements of a rectangular plate on a Pasternak foundation using neural networks to solve the differential equations, but their method required a high number of iterations in the training stage. Vlasov and Leont’ev [9] modelled the soil as a simplified continuum with null horizontal displacements and the vertical displacements approximated by a product of separated functions. Although the soil response obtained in this way is more realistic than that calculated with Winkler-based models, it has the limitation of neglecting the soil horizontal displacements, which produces an artificially stiff foundation response [10]. Subsequently, some authors have applied the Vlasov foundation model to calculate the settlements of a plate using energy * Corresponding author. E-mail addresses: [email protected] (E. Justo), [email protected] (I. Gonz´ alez-de-Le´ on), [email protected] (M. V´ azquez-Boza). Contents lists available at ScienceDirect Engineering Structures journal homepage: www.elsevier.com/locate/engstruct https://doi.org/10.1016/j.engstruct.2025.120005 Received 11 September 2024; Received in revised form 30 January 2025; Accepted 25 February 2025 Engineering Structures 331 (2025) 120005 Available online 5 March 2025 0141-0296/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). methods with a variational approach. Vallabhan and Das [11], and more recently, Gunerathne et al. [12] used variational calculus and finite differences to calculate the settlements of a circular plate on a homogenous soil described with a Vlasov model. Gunerathne et al. compared their results with those obtained with the Finite Element Method (FEM) and reported a reasonably good agreement, but they had to introduce a correction in Young’s modulus and Poisson’s ratio of the soil to compensate for the stiffer response obtained with their model compared to the FEM due to the assumption of null radial displacements [12]. Besides, variational methods have the disadvantage of their complex mathematical formulation. Typically, an iterative procedure is required to obtain the solution [12]. The Boundary Element Method (BEM) has also been applied to address the problem of a plate on a semi-infinite elastic foundation (e.g., [13,14]), but it has the disadvantage of its mathematical complexity. Finally, FEM models (e.g., [15]) provide the most accurate approximation to the problem, but they are too expensive for daily engineering applications. In this context, the development of simpler, flexible and straightforward numerical methods is still an area of research interest for both academic and practical engineering applications [12]. A simpler and yet accurate representation of the problem is that provided by the elastic continuum method [16]. This method uses a closed-form solution of the elasticity equations to model the soil behaviour, while the plate settlements are obtained through a finite-difference approximation of Kirchhoff’s differential equations. Although elastic continuum models have been widely applied to pile problems [17], their use for plate foundations has been much more restricted. Chakravorty and Ghosh [18] solved the problem of a circular plate resting on an elastic half-space using the Boussinesq solution to compute the soil settlements. They derived a general solution for any kind of loading and obtained results for three particular loading cases, but their validation, carried out only against experimental data, was incomplete. Milovic [19] solved the axisymmetric problem using the Schleicher equations [20] to calculate soil settlements. A major drawback of his model is the mathematical complexity introduced by the Schleicher equations, which require evaluating elliptic integrals to obtain the soil settlements. Besides, Milovic provided results for the case of uniform loading but did not validate them with reference methods (e. g., FEM). Melerski [21] solved the problem of a circular plate on a homogenous anisotropic semi-infinite soil using energy methods for the plate and Konig’s equation for the soil [22], but again he did not validate the method against a more accurate solution. In this paper, we present a model for settlement and stress calculation of a circular plate on an elastic foundation based on the elastic continuum approach, which represents the soil behaviour better than Winkler-based models. The method considers the soil beneath the plate as a continuum but, unlike Vlasov models, it does not require neglecting horizontal displacements. Compared to previous studies using the elastic-continuum method [18,19], the main advantages of the proposed model presented are: (1) the use of Mindlin’s equations to calculate soil settlements at any depth, allowing application of the method to multilayered soils with the Steinbrenner approximation; and (2) a comprehensive validation of results against alternative solutions and FEM models. 1.2. Plate resting on multilayered soil The methods described in the previous section consider the soil as a semi-infinite homogeneous, isotropic and linear-elastic continuum. However, most soils are not homogeneous. In real engineering situations, the foundation soil is normally composed of multiple layers with different elastic properties. The consideration of multilayered soil significantly complicates the mathematical formulation of the problem. As a consequence, most studies use FEM (e.g., [23,24]) or BEM models to calculate the settlement of plates resting on soils with multiple layers (e.g., [25]). FEM models provide a good approximation to the problem, but the time and resources required for model preparation and analysis are too great for routine calculations in engineering practice. On the other hand, BEM models have the drawback of their cumbersome mathematical formulation. Very few research studies have presented analytical or semianalytical models for the settlement analysis of a plate on a multilayered soil. Xiao & Yue [26] published an analytical method for a perfectly rigid plate using Fourier integral transforms and Fredholm integral equations, that had to be solved using a numerical technique. Wang et al. [27,28] used Hankel transforms to solve the soil elasticity equations, but their solution includes semi-infinite integrals of the Bessel functions, for which there is no standard method of calculation [27]. Gunerathne et al. [29] presented a variational method for multilayered soils based on their previous study for homogenous soils. Their results showed reasonable agreement with those of the FEM but, due to the assumption of zero radial displacements, their model showed a stiffer response than the FEM solution. El-Garthy and Galel [30] used a Vlasov model to study the behaviour of a rectangular plate on multilayered soil using finite elements and finite differences to solve the equations, again with the assumption of zero radial displacements. Elhuni et al. [10] applied an energy approach using finite elements and finite differences to calculate the settlements of a circular plate on a multilayered soil subject to horizontal and vertical loads. The soil was considered as a continuum Vlasov model with the improvement that horizontal soil displacements were explicitly taken into account. Their results were more accurate than those of previous Vlasov-based models, at the expense of a greater mathematical and numerical complexity. In this paper, we present a numerical method based on the elastic continuum approach to calculate the settlements and contact stresses of a circular plate on a multilayered soil. The method uses Mindlin’s equations [31] to calculate the soil settlements. Unlike Boussinesq or Schleicher solutions, used in previous studies with the elastic continuum method [18,19], Mindlin’s equations give soil settlements at any depth, allowing the model to be extended to the case of multilayered soils with the Steinbrenner approximation [32]. Compared to studies using energy methods and Vlasov type foundation models for the soil [29], the method has the advantage of less mathematical complexity while providing more accurate results in most practical cases. FEM and BEM models, while accurate, may not be a cost-effective option for everyday engineering calculations. The Steinbrenner approximation has been combined with elastic continuum methods in previous studies aimed to calculate the settlement of piles, but it has never been used for plates. Poulos and Davis [17] used this method to calculate the settlements of a pile in a finite layer underlain by a rigid base. In a previous study [33], the authors of this paper extended the elastic continuum method to piles in multilayered soil applying the Steinbrenner approximation, and showed that it provided good agreement with FEM calculations and field tests. The method proposed in this paper has a simpler mathematical formulation than previously published models, and yet provides a more accurate approximation of the soil response for uniform soils and the most part of multilayered soils due to the fact that it does not require simplifying assumptions about the soil displacements or stresses. The method has been validated by comparing results with an FEM solution, showing that the elastic continuum method can be successfully applied to plates on multilayered soils in which the soil stiffness increases with depth, which is the most common situation in real engineering practice. 2. Materials and methods 2.1. Definition of the problem and basic assumptions The problem under study is the settlement analysis of a circular plate of radius R p resting on a multilayered soil foundation. For Mindlin’s equations [31] to be applicable, the soil is initially considered to be a E. Justo et al. Engineering Structures 331 (2025) 120005 2 homogeneous isotropic elastic half-space (Fig. 1). Modifications of the basic analysis for multilayered soil will be introduced later based on the Steinbrenner approximation. The plate and the subgrade soil are divided into elements subject to a constant contact stress. The problem solution is obtained by imposing compatibility of soil and plate displacements. The plate’s material is linear-elastic, homogenous and isotropic, and its thickness is assumed to be small enough to apply Kirchhoff’s theory of thin plate bending. The plate deflections are assumed to be small in comparison with its thickness. The circular plate is supposed to be subject to axisymmetric loading. As the problem is axisymmetric, it is sufficient to consider deflections in one diametrical section of the plate. The plate–soil interface is considered perfectly smooth, which means that a free relative movement is allowed between the plate and the underlying soil mass without the development of any interface shear stress [34]. This assumption is common in plate-soil interaction problems (e.g., [19,26]) and in most cases it represents approximately the real behaviour of the plate-soil contact. 2.2. Plate and soil discretisation The plate and the subgrade soil are discretised as shown in Fig. 2. Each radial sector of the plate is divided into n elements. The circular area of the soil in contact with the plate is also divided into n elements per radial sector. After testing the model discretisation with an increasing number of elements, it was determined that a mesh size value (n) of 36 elements provided sufficient accuracy for the research objectives. A comparative analysis of results using different mesh sizes is shown in Section 3.1.2. The problem is formulated in terms of two sets of variables: the contact stresses p i at the plate-soil interface and the vertical displacements w i at the centre of each element. 2.3. Plate displacement equations The plate displacements were calculated using Kirchhoff’s theory of thin plates. According to this theory, the differential equation for axisymmetric bending of circular plates resting on an elastic foundation has the form [35]: d4w dr4+2 r d3w dr3−1 r2 d2w dr2+1 r3 dw dr=1 D• (q−p)(1) where: q is the uniform load. p is the contact stress at the plate-soil interface. w is the vertical displacement. r is the radial coordinate. E p and ν p are the elastic constants of the plate. D is the flexural rigidity of the plate, defined as follows, where h is the thickness of the plate: D=Eph3 12(1− ν p2) The internal forces acting on a differential element of the plate are (Fig. 3): Mr= − D(d2w dr2+ ν r dw dr)(2) Mt= − D(1 r dw dr+ ν d2w dr2)(3) Q=D(d3w dr3+1 r d2w dr2−1 r2 dw dr)(4) where: M r is the radial bending moment. M t is the tangential bending moment. Q is the shearing force. As shown in Fig. 4, the boundary conditions at the plate’s centre (r=0) are obtained using the equations of axial symmetry of the plate deflections: w−1=w1(5) w−2=w2(6) At the outside edge (r=R p ), free-edge conditions must be applied: Mr=0 (7) Q=0 (8) Eqs. (1)-(4) can be reformulated in finite difference form. Using a Fig. 1. Problem definition. Fig. 2. Plate and soil discretisation with contact stresses acting on the plate and adjacent soil. Fig. 3. Internal forces on a differential element of the plate. E. Justo et al. Engineering Structures 331 (2025) 120005 3 symmetric central difference approximation, the derivatives appearing in Eqs. (1)-(4) can be transformed into algebraic expressions: wIV i=wi−2−4wi−1+6wi−4wi+1+wi+2 δ4(9) wIII i=−1 2wi−2+wi−1−wi+1+1 2wi+2 δ3(10) wII i=wi−1−2wi+wi+1 δ2(11) wI i=−wi−1+wi+1 2δ(12) where δ is the radial dimension of the plate element and i is the element number (from 1 to n). Substituting Eqs. (9)-(12) into Eq. (1), the plate displacement equation for element i is obtained in algebraic form: Awi−2+Bwi−1+Cwi+Dwi+1+Dwi+2=qi−pi(13) where: A=1 δ4−1 rδ3(14) B=−4 δ4+2 rδ3−1 r2δ2−1 2r3δ(15) C=6 δ4+2 r2δ2(16) D=−4 δ4−2 rδ3−1 r2δ2+1 2r3δ(17) E=1 δ4+1 rδ3(18) Eq. (1) can be applied as it is to elements 3 to n-2, whereas for elements 1,2, n-1 and n the displacements at the fictitious elements -1, -2, n+1 and n+2 (Fig. 4) must be eliminated from Eq. (13) using the boundary conditions at the centre and edge of the plate. Using Eqs. (5)-(6) to eliminate the fictitious vertical displacements w - 1 and w -2 we obtained the plate equations for elements 1 and 2. For i=1: (B+C)w1+ (A+D)w2+Ew3=q1−p1(19) For i=2: (A+B)w1+Cw2+Dw3+Ew4=q2−p2(20) To eliminate w n+1 and w n+2 we applied the free-edge boundary conditions (M r =0, Q =0) at the edge of the plate, which required evaluating the first, second and third derivatives [Eqs. (2) and (4)] at r=R p . For that purpose, we used an ad-hoc finite difference formulation to ensure that w(R P ) does not appear in the equations, since the edge of the plate does not belong to the mesh [Eqs. (21)-(23)]. wIII r=Rp =−wn−1+3wn−3wn+1+wn+2 δ3(21) wII r=Rp =wn−1−wn−wn+1+wn+2 2δ2(22) wI r=Rp =−wn+wn+1 δ(23) Using Eqs. (2),(4),(7),(8) and (21)-(23) to eliminate w n+1 and w n+2 , we obtained the plate equations for elements n-1 and n. For i=n-1: Awn−3+Bwn−2+(C−2 k3−k2 E)wn−1+(D+k3+k1 k3−k2 E)wn =qn−1−pn−1(24) For i=n: Awn−2+(B−2 k3−k2 D−k3+k2 k3−k2 E)wn−1 +(C+k3+k1 k3−k2 D+k1E+k2 k3+k1 k3−k2 E)wn=qn−pn (25) where: k1=1+2 ν δ Rp (26) k2=1−2 ν δ Rp (27) k3=3+(1+ ν )δ2 R2 p (28) The system of algebraic Eqs. (13),(19)-(20) and (24)-(25) can be expressed in matrix form: [Ip]⋅[w] = [q] − [p](29) where [I p ] is the n by n matrix of coefficients for plate settlement, [q] is the vector of distributed loading at the plate and [p] is the vector of contact stresses. 2.4. Soil displacement equations The soil displacements are calculated using Mindlin’s equation for the displacement caused by a point load within a homogeneous halfspace [31] for the particular case of a force acting on the soil surface: w=P(1+ ν s) 8 π Es(1− ν s)[8(1− ν s) R1 +4z2 i R3 1](30) where: E s and ν s are the elastic constants of the soil. R 1 , r and z i are defined in Fig. 5 and have the following relationship: Fig. 4. Boundary conditions at the plate. E. Justo et al. Engineering Structures 331 (2025) 120005 4 R2 1=r2+ (zi)2 The soil displacement I s ij at point i due to a uniform unit stress acting on element j is obtained by integrating Eq. (30) along the surface of the jth element. For each element j, the integration domain is an annulus (Fig. 6), where: r2=x2 i+ ρ 2−2 ρ xicosθ(31) Is ij =∫2 π 0∫Ro Ri w ρ dθd ρ (32) and R i and R o are the inner and outer radius of the annulus All integrals were evaluated numerically using software MATLAB. The singularities that appear in the vertical displacement function w [Eq. (30)] when i=j did not cause difficulties with the numerical integration routine in software MATLAB because they occur at the boundaries of the integration domain. By adding up the contributions of the stresses p j of all elements in the plate-soil interface (j=1…n), the vertical displacement of the soil at the centre of element i can be expressed as: wi=∑n j=1Is ij⋅pj(33) where Is ij is the vertical displacement factor for element i due to the contact stress p j , calculated with Eq. (32). Eq. (33) can be written in vector form: [w] = [Is]⋅[p](34) 2.5. Compatibility of displacements Equalling the plate and soil displacements and combining Eqs. (29) and (34), a compatibility equation in terms of displacements is obtained: [[Ip] + [Is]−1][w] = [q](35) By solving the system of equations, a solution for w can be obtained. Subsequently, the contact stresses p can be derived from Eq. (34): [p] = [Is]−1⋅[w](36) Using Mindlin’s equation for σ z [31], the vertical stress at any depth within the soil can be obtained integrating Eq. (37) in a similar way to that used to calculate the vertical displacements [Eq. (32)]: σ z=3 2 π Pz3 R5 1 (37) 2.6. Consideration of layered soil The method can be extended to include multilayered soil problems by means of the Steinbrenner approximation [32]. The soil displacement at an element i in layer k caused by a unit stress acting on element j can be calculated as the sum of the settlements of the layers below. For example, in the case depicted in Fig. 7: Is ij =Is ij(Ek) − Is kj(Ek) + Is kj(Ek+1) − Is k+1,j(Ek+1) + Is k+1,j(Ek+2) − Is k+2,j(Ek+2) where [Is hj(Eh) − Is h+1,j(Eh)]is the settlement of any layer h calculated with the Eq. (32) using the modulus of elasticity E h of that particular layer. 3. Method validation 3.1. Circular plate in homogenous soil The proposed method has been verified by comparing the results with those recently published by Gunerathne at al. [12] for two cases of circular plates in homogenous soil (Fig. 8), and also with an FEM analysis carried out with PLAXIS 2D software. Case A consists of a perfectly rigid concrete plate with R p =0.5 m, h p =150 mm, E p Fig. 5. Effect of a vertical force acting on the surface of a semi-infinite solid. Fig. 6. Integration scheme for element j at the soil surface. Fig. 7. Model for multilayered soil. E. Justo et al. Engineering Structures 331 (2025) 120005 5 =30 GPa, and ν p =0.15. Case B features a perfectly flexible steel plate with R p =5 m, h p =10 mm, E p =200 GPa, and ν p =0.3. Both plates are subjected to a uniform load q of 60 kPa. The soil parameters are in both cases: E s =20 MPa, and ν s =0.2. The PLAXIS model for soil and plate is shown in Fig. 9 for the rigid plate case. The axisymmetric model was built using 15-node triangular soil elements and 5-node plate elements [36]. The depth and width of the soil model were adjusted to 200R p and 50R p , respectively, to avoid disturbance in the solution due to edge constraints. The relative increase in settlements obtained when we raised the soil’s depth from 50 to 200R p was almost 5 % in some cases, which is a significant difference in the context of this research. Horizontal displacements were constrained along the axis and the outer edge, while fixed-ended boundary conditions were defined at the bottom of the soil. The type of mesh used was “very fine”, with a subsequent mesh refinement in the plate area. The numerical model comprised 9238 elements and 74770 nodes. We used a Linear Elastic material model for the plate and the soil. The settlement results for the rigid and flexible cases are compared in Fig. 10, showing that in both cases our solution closely matches the one obtained with PLAXIS. Besides, the settlements calculated with this method show better agreement with the FEM results than those obtained by Gunerathne et al. [12] using a more complex variational solution. An additional advantage of the elastic method is that it does not require the correction of the soil elastic parameters, since radial displacements are taken into account in the derivation of Mindlin’s equations from Elasticity theory. Fig. 11 shows that the contact stresses calculated with the elastic method for the (a) rigid and (b) flexible plates are both in very good accordance with the FEM results. Only local discrepancies can be observed near the plate edge, where the stress gradient is very high. Comparison with the variational solution published by Gunerathne et al. is not possible as their model did not explicitly calculate stresses. In real engineering situations, the foundation plates are not necessarily perfectly rigid or perfectly flexible, so we also tested the method for one plate of intermediate rigidity. To compute the plate-soil relative stiffness we have used the relative raft stiffness defined by Brown [37]: KR=Ep Es(1− ν 2 s)(h Rp)3 (38) The settlement results and contact stresses of a plate of intermediate rigidity (K R =0.1) are portrayed in Fig. 12. In this case too, the results obtained with the elastic continuum method are in close agreement with the FEM solution. 3.1.1. Sensitivity analysis We conducted a sensitivity analysis to investigate how the stiffness of the plate influences the method’s accuracy. For that purpose, we calculated settlements and contact stresses of plates with values of relative stiffness (K R ) of 10, 1, 0.1 and 0.01. For values of K R greater than 10 and smaller than 0.01, the plate behaviour can be assimilated to the Fig. 8. Cases for model validation: (A) Perfectly rigid plate. (B) Perfectly flexible plate. Fig. 9. Characteristics of the PLAXIS model. E. Justo et al. Engineering Structures 331 (2025) 120005 6 perfectly rigid or perfectly flexible cases, which have already been discussed. The results of the analysis are plotted in Fig. 13. From Fig. 13(a) it can be observed that the settlements under the plate approach the FEM solution for stiffer plates, whereas the discrepancy is higher for more flexible plates. Meanwhile, the settlements at radial distances outside the plate present a constant deviation from the FEM solution, irrespective of the plate stiffness. In fact, the results obtained with both methods are independent of K R for radial distances greater than 1.5 R P approximately. As for the contact stresses (Fig. 13b), the results were in close agreement with the FEM solution in all the cases. Only small variations were observed at the area of high stress gradient near the edge of the plate, where a finer mesh would be required to achieve a perfect match. 3.1.2. Discretization analysis A mesh convergence study with an increasing number of elements has been carried out to find the minimum number of elements for which the discretisation has no influence on the solution. We tested meshes with 9, 18, 36, and 72 elements. The results in terms of settlements and contact stresses are plotted in Fig. 14 for the perfectly rigid plate and in Fig. 15 for the perfectly flexible plate. It can be observed that there is practically no difference in the settlements for meshes of 9 elements and above. As for the contact stresses, the results only diverge at the high stress gradient area near the edge of the plate but even in that region the stresses converge acceptably from 36 elements on. As a consequence, we used a 36-elements mesh for the purpose of this study, although a mesh of only 9 elements would be perfectly adequate for practical engineering calculations where high computational accuracy is not required. The low discretisation requirement is an additional advantage of the proposed numerical method. 3.2. Circular plate in multilayered soil The performance of the model in a multilayered soil was verified by comparing the results obtained with those published by Gunerathne et al. [29] using the variational method and also with an FEM solution obtained with the PLAXIS software. Fig. 16 shows the multilayered soil model used for validation. Fig. 10. Settlements of (a) a perfectly rigid and (b) a perfectly flexible circular plates on homogeneous soil. Fig. 11. Contact stress σ z distribution of (a) a perfectly rigid and (b) a perfectly flexible circular plates in homogeneous soil. E. Justo et al. Engineering Structures 331 (2025) 120005 7 Fig. 17(a) shows that settlements calculated with this method are in better agreement with the FEM solution than those obtained with the variational method. Likewise, as shown in Fig. 17(b), the contact stresses are in very close accordance with the FEM results. 4. Results and discussion 4.1. Circular plate on three layered soil We tested the performance of the proposed numerical model for three multilayered soil profiles (Fig. 18) whose relative stiffnesses reproduce those used by Poulos in a classic study of piles in multilayered soil [38]. The soil stiffness increases with depth in soil 1 and decreases with depth in soil 3. Soil 2 represents an intermediate situation. The results are compared with an FEM solution calculated with PLAXIS. The multilayered PLAXIS model is similar to the one described in a previous section used for homogenous soil. For each soil profile, we have calculated the settlements and contact pressures of three plates of different stiffnesses: (a) rigid plate, (b) flexible plate and (c) plate of intermediate rigidity (K R =0.1). The dimensions and elastic parameters for the rigid and flexible plates (a) and (b) are identical to those used for homogeneous soil (see Fig. 8). For the plate of intermediate rigidity (c), the plate parameters are E p =30 GPa, ν p =0.15, R p =10 m, h =0.4 m. The settlement results for the perfectly rigid plate are plotted in Fig. 19. For soil 1 and soil 2 the settlements obtained with the elastic method are in acceptable agreement with those calculated with FEM, with relative errors below 10 % (Table 1). For soil 3, where soil stiffness decreases with depth, our results differ significantly from those produced by FEM, with the relative error reaching 16.7 %. For the perfectly flexible plate (Fig. 20) and the plate of intermediate stiffness (K R =1) (Fig. 21) the trend is similar: the calculated settlements are reasonably accurate (error <10 %) for soils 1 and 2, whereas for soil 3 the errors reach 25 % (Table 1). Additional calculations carried out with plates of K R values of 1 and 0.01 provided similar outcomes. These results are coherent with previous findings by Poulos [38] who in his study of piles in multilayered soil also reported greater errors when approximate solutions based on the elastic method were applied to soils with softer underlying layers. Fig. 22 to Fig. 24 display the contact stresses for the rigid, flexible and intermediate plates in the three multilayered soil profiles analysed. It can be seen that the results in terms of contact stresses show a good agreement with the FEM solution in all cases. Only small discrepancies Fig. 12. (a) Settlements and (b) contact stresses of a circular plate of intermediate rigidity in homogeneous soil. K R =0.1 (Ep =30 GPa, ν p =0.15, Rp =10 m, h = 0.4 m, Es =20 MPa, ν s =0.2). Fig. 13. (a) Settlements and (b) contact stresses of circular plates of different relative stiffness K R =10, 1, 0.1 and 0.01 in homogeneous soil. E. Justo et al. Engineering Structures 331 (2025) 120005 8 can be observed for soil 3 due to the inaccuracy of the Steinbrenner approximation for soils with stiffness decreasing with depth, as we mentioned above. {{{Fig. 23, Fig. 24}}} Based on these results we can conclude that the elastic continuum method with Steinbrenner’s approximation yields satisfactory outcomes for cases where the soil stiffness increases with depth, which is the most common scenario in real engineering applications. However, the method is not suitable for soils with denser layers overlying softer layers (i.e., road pavement foundations). For these type of soils, Steinbrenner’s approximation cannot properly reproduce the multilayered soil behaviour and a more elaborated method should be used. For example, models using a variational approach have been shown to be effective in previous studies [29]. 4.2. 4.2 Circular plate on two layered soil To investigate why Steinbrenner’s approximation cannot replicate the behaviour of soils with denser upper layers, we analysed the problem of a plate resting on a soil with only two layers that have markedly different moduli of elasticity (Fig. 25). In soil A, the softer layer is above, whereas in soil B the denser layer is above. On this occasion the focus is on the soil response, for that reason we have used a perfectly flexible plate for the analysis. The influence of the plate is therefore minimised, and the stress distribution at the plate-soil interface approaches the Fig. 14. (a) Settlements and (b) contact stresses of a perfectly rigid circular plate on a multilayered soil with different values of mesh sizes (n). Fig. 15. (a) Settlements and (b) contact stresses of a perfectly flexible circular plate on a multilayered soil with different values of mesh sizes (n). Fig. 16. Multilayered soil model used for validation. E. Justo et al. Engineering Structures 331 (2025) 120005 9