Finite element method algorithm for geotechnical applications based on Runge-Kutta scheme with automatic error control

Abstract This paper introduces a novel explicit algorithm to solve the finite element equation linking the nodal displacements of the elements with the external forces applied via means of non-linear global stiffness matrix. The proposed method solves the equation using Runge-Kutta scheme with automatic error control. The method allows any Runge-Kutta scheme, with the paper demonstrating the algorithm efficiency for Runge-Kutta schemes of second to fifth order of accuracy. The paper discusses the theoretical background, the implementation steps and validates the proposed algorithm. The numerical tests show that the proposed method is robust and stable. In comparison to the iterative implicit methods (e.g. Newton-Raphson method), the new algorithm overcomes the problem of occasional divergence. Furthermore, considering the computation time, the fifth-order accurate scheme proves to be competitive with the iterative method. It seems that the proposed algorithm could be a powerful alternative to the standard solution procedures for the cases with strong nonlinearity, where the typical algorithms may diverge.


Introduction
The Finite Element Method approximates the solution of the continuous mechanical balance equation by calculating the unknown displacement u at a set of discrete points only, typically at the element nodes, as a response to the external load f ext . This requires solving a large set of equations, with coefficients given in the global stiffness matrix K where: In most engineering applications, the stiffness matrix K is changing in a non-linear fashion with the increment of external load f ext . The stiffness matrix nonlinearity can be a consequence of both the material behaviour and the geometric nonlinearity. In the common geotechnical applications, the material behaviour nonlinearity is dominant. Therefore, this paper disregards geometrical nonlinearity.
To solve Eq. (1), the Finite Element Method needs an efficient, stable and robust numerical algorithm. Typically the method of choice belongs to the family of iterative Newton-Raphson scheme, for example: modified, initial stress or accelerated (Bathe, 2006), which usually offers a satisfactory rate of convergence. Yet, the Newton-Raphson scheme requires the initial iteration result close enough (i.e. within the convergence radius) to the correct solution to converge and an even closer result in order to get the quadratic rate of convergence. On the other hand, if the load increment is large compared to the nonlinearity of the problem, the initial iteration result may fall outside of the convergence radius and the Newton-Raphson scheme may fail to reach the solution. If this happens, the codes after significant number of iterations usually cut the load increment to half and repeat the process, hoping that the convergence will be reached. If that does not work, ultimately, the solution may be abandoned with an error (Section 8 gives an example of such situation).
The method also requires the evaluation of the inverse of the Jacobian (first order derivatives) matrix at each iteration which imposes additional numerical expense and difficulties in a large system of equations. To avoid the problems related to the repetitive calculation of the Jacobian, the quasi-Newton methods evolved (e.g. Broyden-Fletcher-Goldfarb-Shanno (BFGS)) (Matthies and Strang, 1979;Avriel, 2003). The main idea behind this class of methods is that the Jacobian matrix in Newton-Raphson scheme needs to be estimated only in the first iteration and then, in the subsequent iterations, it will be approximated based on the Jacobian from the previous iteration. In the original contribution, Brayden (Broyden, 1965) even proposed a formula to update the Jacobian inverse directly based on the inverse of the preceding iteration which would further save the computational effort. Although the quasi-Newton methods offer occasionally successful https://doi.org/10.1016/j.compgeo.2020.103841 Received 4 March 2020; Received in revised form 10 September 2020; Accepted 12 September 2020 alternative when the classical Newton-Raphson scheme fails, they still require the initial guess to be close to the correct solution to ensure convergence. Gens and Potts (1988) carried out an early study on the different methods used to solve nonlinear equations in geotechnical applications with a constitutive model based on the critical state concept. The study included, in addition to the iterative procedures, the explicit first-order forward Euler scheme (tangent stiffness method). They concluded that there is no generic statement on the recommended method to be used and it is necessary to equip the numerical code with different solution methods so that the most suitable one can be used depending on the solved problem. In a later study, however, Potts and Ganendra (1992) concluded that the modified Newton-Raphson method is the most robust and economical method, which is now the standard in most of the Finite Element codes. To improve the convergence of the scheme, several numerical techniques can be used, including the arc length control and the line search method (Bathe, 2006;Crisfield, 1983). Abbo and Sloan (1996) provided an improved second-order explicit scheme (Euler-modified Euler) with error control of the drift from the load path. Although, the scheme showed to be robust and stable even in applications that involve critical state soil models , the iterative methods were significantly faster, hence the scheme was not widely adopted. Recently, the numerical software started including the automatic (algorithmic) differentiation (AD) (Griewank and Walther, 2008), which means that the accurate derivation of the stiffness matrix is available without extra numerical burden. This opens new possibilities for novel numerical algorithms, more stable and robust than existing ones.

Importance and significance of the proposed method
This paper presents a novel alternative method which treats Eq. (1)  Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 as a differential equation and applies Runge-Kutta explicit schemes with automatic sub-incrementation (sub-stepping) and error control to solve it. The paper discusses first the theoretical background and the numerical implementation steps. Subsequently, the method is tested on a typical geotechnical problem of a shallow foundation resting on an elastoplastic soil. In the tests, the new method has been more stable than most common existing solutions. The proposed algorithm based on the 5th order Runge-Kutta scheme is also competitive with the iterative schemes regarding the time of calculations.

Roman
The new method of solution builds on a well-tested, robust scheme with mathematically proven convergence properties. Consequently, its adoption may be beneficial not only in the Finite Element software for geotechnical engineering, but in all Finite Element codes used for analysis of highly non-linear problems, for instance coupled thermohydro-mechanical-chemical analyses.

Solution algorithm
In the case of nonlinear material behaviour, the stiffness matrix in Equation (1) depends on stresses and strains. Those are related to displacements, hence the stiffness matrix is ultimately displacement de- , where f int is the vector of internal nodal forces. This allows us to rewrite Equation (1) as a differential equation: where u d i is the infinitesimal increment of displacements and f d i is the infinitesimal increment of externally applied forces over an algorithmic time infinitesimal increment t d . The dependency of stiffness matrix on displacements stems from the complex constitutive behaviour of soil (material) which, in most cases, involves nonlinear elasticity and plasticity. Typically, Eq. (2) is solved incrementally allowing for linearization by applying the total external force in increments f i per each time increment t. To enhance the mathematical clarity of Eq. (2), one introduces the dimensionless time = t t t T ( )/d 0 where t and t 0 are the algorithmic time at the end and at the start of the load increment resulting in 0 T 1. Then, following the chain rule, the incremental displacements per time increment will be (2) yields:

Substitution in Equation
Eq.
(3) is a first-order ordinary nonlinear differential equation which can be solved using iterative methods (implicit) such as the Newton-Raphson method. Other possibility is to use the explicit methods that solve the equation directly without iterations but assume load increments f i that satisfy certain accuracy condition. Fig. 1 and Fig. 2 show, for a one-dimensional case, how both methods approach the solution u i for the i th load increment f i . Based on Eq.
(3), for each load increment f i , one calculates the corresponding displacement increment u i and adds this value to the previously calculated total displacements u i 1 to determine the new u i . The solution starts with known initial displacement u 0 , thus Eq. (3) is an initial value problem. Despite the rapid convergence rate of the implicit methods, their convergence is conditional as the initial prediction needs to be sufficiently close to the correct solution. Therefore, the Newton-Raphson method divergence happens especially often in the case of material behaviour with strong nonlinearity. Another major disadvantage of the implicit methods is that they do not have any load path error control meaning that it could converge if big loading increments are adopted, but to a wrong answer (Abbo and Sloan, 1996). On the other hand, the explicit methods are stable; however, the only explicit method commonly used in the past, the Forward Euler method, tends to drift away from the correct solution with the advancement of the solution. To control the drift, small load increments should be adopted which slow down the solution tremendously and make the Forward Euler method unattractive.
Although the explicit Runge-Kutta methods received wide recognition for integrating constitutive laws on the Gauss point level (Sloan, 1987;Sloan et al., 2001;Sołowski and Gallipoli, 2010), they got little attention on their applicability for solving the global finite element equations. Apart from the contribution by Abbo and Sloan (1996) who proposed an explicit scheme with error control that can be considered as a second-order scheme when compared to the first-order Forward Euler scheme, the authors are not aware of any other major contribution.
When solving Eq. (3) with Runge-Kutta method, the estimation of u i is taken as a weighted sum of explicit evaluations of u i k ( ) at predefined positions k along the given algorithmic time interval T. The number of these evaluations is directly related to the required accuracy of the solution, yielding a scheme with different orders (e.g. first, second … fifth and higher). Employing a high order Runge-Kutta scheme leads to a fast convergence and an accurate solution with well controlled errors. It should be noted that the Runge-Kutta method is a standard numerical method that can be reviewed elsewhere, see Lee and Schiesser (2003) for example.

Runge-Kutta explicit scheme for load-deflection estimation
The first-order Runge-Kutta method is equivalent to the forward explicit Euler method where the displacement increment is the direct outcome of Eq. (3) with = dT 1, see Fig. 2 and Eq. (4). Obviously, the accuracy of the first-order scheme is heavily dependent on the load increment size. To achieve a sufficiently accurate solution, very small increments should be adopted leading to a numerically expensive solution. For Runge-Kutta second-order scheme, the estimation of the displacement increment u i needs two stages = k ( 1, 2) of secondary displacement increment evaluations u i k ( ) (see Fig. 2 for a graphical clarification): Stage = k 1: 1 and (1) and f int (2) represent the internal forces at the corresponding displacements u i 1 and + u u i i 1 (1) , respectively. Consequently, the second order estimation of the displacement increment u i is: As can be noticed, the weighing factors are identical here for the two estimations and equal ½. One important aspect to mention is that several Runge-Kutta methods have embedded a lower order solution, which can be used to estimate the error. This error is given as the difference between the current scheme displacement estimation and the estimation with the one order lower scheme. For example, for the second-order scheme presented above, the error Err is given as: Runge-Kutta schemes with higher order of accuracy follow a similar algorithm. By knowing the correct weights and positions to estimate each secondary stage evaluation, the displacement increment is determined as a weighted sum of the performed evaluations. For an arbitrary high order method, Eq. (4) and Eq. (5) read: Here u i k ( ) represents the updated displacements according to the adopted Runge-Kutta scheme. The stiffness matrix In this contribution, the evaluation of the stiffness matrix K u ( ) i k ( ) is carried out using the automatic differentiation (AD) during the finite elements assembly process as will be discussed later. The estimated displacements are then calculated as: where s is the total number of evaluation stages that is dependent on the scheme order n (Lee and Schiesser, 2003). By subtracting the lower order solution from the one order higher solution, the error can be estimated as: Tables 1-6 give the coefficients kj , k and k for Eqs. (9), (10) and (11) for Runge-Kutta schemes of 2nd, 3rd, 4th and 5th order. Note that any Runge-Kutta scheme, e.g. those given in (Lee and Schiesser, 2003), can be used -tables give only the coefficients for schemes tested in this paper.

Runge-Kutta explicit scheme with error control
The error in the Runge-Kutta scheme is dependent on the size of the load increment f i . For the same increment size, the higher-order schemes yield predictions that are more accurate. Additionally, the method assesses the error which can be used to automatically choose the next load increment size, so that the error is bound within a desirable range, similarly as in the case of stress integration (Abbo and Sloan, 1996;Sloan, 1987;Sołowski and Gallipoli, 2010). The proposed steps of the method to solve the global finite element equations are illustrated in Diagram 1 and listed in Table 7.
The scheme starts by assuming that the first load sub-increment where j = 1 is the number of the current sub-increment, equals the full load increment f i by using an initial sub-increment size of T j = 1. After calculating the corresponding displacement using f j by employing Eq. (9) and Eq. (10), the occurred relative error Error is estimated by applying Eq. (11). If this error is less or equal to the imposed tolerance D Tol by the user, the current sub-increment is accepted and the program moves to the next load increment (see Diagram 1 and Table 7). On the other hand, if the error is greater than D Tol , the current sub-increment is rejected and instead of T j = 1, a new smaller value is estimated using Eq. (12). At this point the calculations are repeated with the new sub-increment size. The algorithm will continue reducing the increment size until the error criterion is satisfied and the step is accepted.
For both accepted and rejected sub-increments, a new sub-increment size is estimated again using Eq. (12). The procedure continues until the full increment is applied yielding and n sub is the resulted total number of accepted sub-increments.
The size of the new sub-increment T new is estimated based on the old given sub-increment size T old , the scheme order of accuracy n and the estimated error Error in the current sub-increment. Thus, the new Table 1 Stages for the 2nd order scheme.

Table 2
Stages for the 3rd order scheme.

Item Estimate
Stage 1 displacement estimate: (2) 1 Third order displacement estimate: A.A. Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 sub-increment size is Here, is the relative change in the next sub-increment size, D Tol is the required accuracy imposed by the user and is a safety factor, as this sub-increment size estimation is approximate. Such a choice of the new sub-increment size should keep the error in the next step below the tolerated value (Gustafsson, 1992). To reduce the chance of spurious increase/decrease of the steps due to random correlation between the Table 3 Stages for the 4th order scheme.

Table 6
Coefficients for the 5th order scheme. two solutions used for estimation of the error and to avoid unnecessary failed sub-increment sizes, the resulted from Eq. (12) is limited by value of and 0.1. Previous experience (Abbo and Sloan, 1996;Sloan et al., 2001) showed that the recommended values are in the range 0.7 0.9 and 1.1 2.0. The authors used the values = 0.9 and = 1.2, leading to 0.1 1.2 in the calculations. This means that the increase of the next load sub-increment is capped at 120% of the current load sub-increment. As can be interpreted from Eq. (12), the new sub-increment size will be reduced in case of rejected step ( > Error D Tol ).

Note on the evaluation of the global stiffness matrix by automatic differentiation
The evaluation of the global stiffness matrix is an important step during the solution for displacement increment as given by Eq. (8). Typically, that would need the assembly of the global stiffness matrix from the individual finite element stiffness matrices respecting the relevant degrees of freedom (Smith et al., 2013). Thebes code (Abed andSołowski, 2017, 2019), based on Numerrin numerical solver (Laitinen, 2013), follows a different method to evaluate the stiffness matrix by the direct differentiation of the internal nodal forces with respect to the current nodal displacements. That is carried out using the automatic (algorithmic) differentiation (AD) (Griewank and Walther, 2008) where the stiffness matrix is evaluated during the assembly of the internal nodal forces. The main idea behind AD is that any function is executed as a sequence of basic operations (e.g. additions, subtractions, multiplications, etc.). By applying the chain rule repeatedly to these operations, the derivative of such a function can be computed with no truncation errors (i.e. to the machine precision). In order to achieve that, a parallel program translates the code of the respective function (that is not necessarily present in closed form but most of the time as a computer program) into a new code that contains derivatives. The application of this method is not well recognized in the field of computational geomechanics but has received full consideration in the general field of the finite element method (Tijskens et al., 2002; Rothe and Diagram 1. A flow chart illustrates the required steps by the proposed Runge-Kutta scheme.

A.A. Abed and W.T. Sołowski
Computers and Geotechnics 128 (2020) 103841 Hartmann, 2015; Zwicke et al., 2016;Mitusch, 2018). The theoretical details of algorithmic differentiation are out of the scope of this paper, however the interested reader is referred to (Griewank and Walther, 2008;Bischof et al., 1992Bischof et al., , 1996 for more details. In any case, the way the global stiffness matrix is determined does not affect the essence of the proposed method for solving the finite element equations.

Comparison of a single Newton-Raphson iteration versus one Runge-Kutta calculation stage
To assess the efficiency of the proposed scheme, we compare it to the Newton-Raphson scheme, typically used in Finite Element Method calculations. During a standard full Newton-Raphson iteration k within the load increment (step) i, the following operations are performed: 1. Form the tangent stiffness matrix K T (Jacobian matrix) at the estimated displacement from previous iteration: 3. Compute the new displacements + u i k 1 to be used in the next iteration where = r f f ext int is the force residual (out of balance) vector. The iterations continue until the imposed convergence criterion is satisfied.
On the other hand, the operations for one stage k of Runge-Kutta schemes follow the steps: 1. Evaluate the stiffness matrix for a stage k using the estimated displacements for that stage: Table 7 Required steps to integrate the global finite element equations using the Runge-Kutta scheme.
1. Enter with the estimated total displacements from the previous loading increment ui 1, the applied load fi 1, the new load increment fi. Provide values for the tolerated error D Tol , the minimum load increment size ΔT min , the safety factor for load increment size χ and the threshold value to control the maximum growth of the next load increment α. Choose the required order of the Runge-Kutta scheme (n). Initialize T = 0 and T 1 = 1. 2. Loop as long as < = T 1 and perform a to d, otherwise, go to Step 3 a. Based on the chosen order of Runge-Kutta scheme (n), evaluate the displacements for each stage k using formula (9)  c. Estimate the error vector Err n ( ) based on Eq. (11) depending on the chosen scheme order, then estimate the relative error: For the 2D case, the error norm Err n ( ) is calculated as: where Err is the displacement error component being estimated, according to each order, based on Tables 1-6 and m is the total number of nodes in the mesh. Note that Err (0) and Err (1)

Calculate the displacement increment using the estimated stiffness matrix in
Step 1 and the provided loading increment:   A.A. Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841   8 3. Calculate the displacements to be used in the next stage + k 1: By comparison, one notices that the number of computational operations in a single iteration of Newton-Raphson method is roughly comparable to that performed during one stage of a Runge-Kutta scheme (though the Runge-Kutta method does require more summation operations in Step 3 to estimate the displacements at the required Table 11 Time and number of sub-increments required by each scheme for the shallow footing problem with D Tol = 1.0E−1.    Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 position for evaluating the stiffness matrix). That gives us a way to approximately compare the computational time in both methods by comparing the number of required iterations and the required calculation stages in order to solve a given problem.

Numerical applications
The following examples discuss the performance of the proposed method based on simulations of the behaviour of a shallow foundation on an elastoplastic soil. The section provides simulations results and computation time with both the proposed schemes and the Newton-Raphson iterative scheme.
It is worth noting that the steps of the proposed method are independent of the adopted material behaviour model on Gauss point level. Thus, in principle, any suitable constitutive model for soil behaviour can be used including linear elasticity, nonlinear elasticity, elastoplastic etc. However, in this paper, in order to test the convergence of the method in case of a challenging material behaviour, the soil is modelled as Modified Cam Clay (MCC) which is an elastoplastic model based on the principles of critical state soil mechanics (Schofield and Wroth, 1968;Wood, 1990). Assuming a constant Poisson's ratio μ, the model requires the following soil parameters: 1) the initial specific volume v 0 at the reference isotropic effective pressure p′ = 1 kPa; 2) the elastic unloading-reloading index κ to reproduce the elastic behaviour; 3) the plastic compression index λ to capture the behaviour during plasticity (both κ and λ are estimated in the plane v-ln p′, see Fig. 3a); and 4) the slope of the critical state line M to predict failure (see Fig. 3b). The MCC is a volumetric hardening model where the hardening parameter is the isotropic preconsolidation pressure p c and the yield function F is given as an ellipse in p′-q plane (see Fig. 3b):  A.A. Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 where the isotropic pressure In the case of overconsolidated clay, the initial isotropic preconsolidation pressure is calculated based on the measured pre-overburden pressure POP (the maximum vertical effective stress ever experienced by the soil) and Eq. (13) yielding: where K o NC is the coefficient of at-rest earth pressure as predicted by MCC model for normally consolidated clay (Brinkgreve and Vermeer, 1992): The stress integration in case of MCC model can be carried out using serval standard schemes including the explicit schemes with error control (Sloan, 1987;Sołowski and Gallipoli, 2010) and implicit schemes (Borja and Lee, 1990;Abed, 2008). Thebes code offers both options to integrate MCC model, however in the following examples only the implicit scheme for the local stress integration is used.

Circular footing on Modified Cam Clay soil
The first example is an analysis of a circular shallow foundation. The 1.0 m radius footing rests on a dry homogeneous slightly over-consolidated clay with POP = 20 kPa. This POP value, based on Eqs. (14) and (15), leads to an initial isotropic overconsolidation pressure = p 18.25 c kPa. The unit weight of the clay is assumed to be 17.0 kN/m 3 and the stresses were initialized with at rest soil pressure coefficient Table 8 lists the adopted MCC parameters in this analysis. The footing applies a uniformly distributed stress of 100 kPa directly on the  A.A. Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 soil in axisymmetric conditions, representing the case of a perfectly flexible footing. The theoretical ultimate bearing capacity of the footing is around 110 kPa (Vesic, 1975). Fig. 4 shows the finite element model and the problem geometry. The mesh consists of 2100 4-noded quadrilateral finite elements with 4 stress integration points per element. Standard mechanical boundary conditions are applied to the problem boundaries with constrained horizontal deformations on the vertical boundaries and fully constrained deformations at the bottom boundary. The numerical analysis is carried out using the finite element code Thebes (Abed andSołowski, 2017, 2019). Originally, the code used full Newton-Raphson scheme to solve the finite element equations, which will be used here as a reference for comparison purposes. The 100 kPa external stress is applied as one loading increment, and then the automatic scheme divided it into the appropriate sub-increments to satisfy the required accuracy using the desired order of the Runge-Kutta scheme. Automatic sub-incrementation is also used for the Newton-Raphson method to constrain the maximum relative error to I TOL = 1.0E−3. It is worthwhile to mention that the Newton-Raphson sub-incrementation follows similar steps to that in Table 7 with three main differences: 1) The displacements in Step 2 are estimated using Newton-Raphson iterations; 2) The "Error" estimation in phase (c) of Step 2

Discussion of the numerical results
Tables 9-11 report the numerical results of the analyses obtained with a range of displacement error tolerance D Tol and a range of Runge-Kutta schemes of different orders. The Tables give the number of required sub-increments, the number of stages and the relative time with reference to t ref -the time required by the 5th order scheme. Figs. 5-7 graphically show the applied footing pressure versus the Euclidian norm of displacement in the studied domain for the different methods adopted in this study. Each graph includes the results using the Newton-Raphson method for the comparison purposes. For example Fig. 8, Fig. 9 and Fig. 10 illustrate the predictions of 3rd, 4th and 5th order schemes with D Tol = 1.0E−1, 1.0E−2 and 1.0E−3, respectively. The results show that for D Tol = 1.0E−3, all Runge-Kutta schemes approximated successfully the correct solution. However, the number of the required sub-increments to obtain the solution differs significantly from one scheme to another. It is worth noting that in the case of the 5th order scheme, the deviation from the correct solution is relatively small, no matter which tolerance is adopted. For this example, the full Newton-Raphson algorithm diverged around the applied footing stress    Abed and W.T. Sołowski Computers and Geotechnics 128 (2020) 103841 of 75 kPa illustrating one of the significant merits of the new method, which always succeeded to complete the calculations and provide the results without numerical difficulties. In addition, the 5th order scheme is considerably faster than the Newton-Raphson scheme as it succeeded to compute the results for the full 100 kPa loading in about 25% faster than the time needed by the iterative scheme to reach 75 kPa loading before diverging. That is confirmed by using the approximate assessment of the number of calculations based on count of the total number of calculation stages and iterations. The 5th order scheme completed the calculations with 108 calculation stages versus 131 iterations done by the Newton Raphson scheme before it diverges and triggers continuous unsuccessful attempts to converge by reducing consecutively the load increment size. Additionally, the Newton-Raphson scheme could easily reach 200 iterations if the iterative method would be able to converge and compute the final results (corresponding to full load increment applied). Table 9, Table 10, Table 11 and Fig. 8, Fig. 9, Fig. 10 illustrate the results of the different schemes but with a similar tolerance. The graphs clearly show the significant differences in the required number of subincrements (illustrated by dots) in each scheme to reach the prescribed accuracy. It is also clear that the 5th order scheme is the most accurate and requires the lowest number of sub-increments which verifies the theoretical expectations. Fig. 11 shows the sub-increment sizes with loading for the case of the 5th order scheme and D Tol = 1.0E−3. The graph illustrates the general tendency of the error control scheme to accelerate the increase of the load increment size as long as the load-displacement curve is smooth, and the slope is not changing dramatically (pure elasticity). Once the curve changes in a non-smooth manner (soil becomes plastic as stress reaches the preconsolidation pressure value, which is a nonsmooth transition), the load increment size reduces in order to control the error. Once plastic behaviour is present, the load-displacement curve becomes smooth again, the scheme starts accelerating and the step size increases quickly again. Fig. 12 shows that the occurred relative error is kept below 1.0E−3 throughout the loading sequence which verifies the capability of the procedure to constrain the error below the prescribed threshold.

Comparison to commercial codes results
Exactly the same problem was modelled using Plaxis2D (Brinkgreve et al., 2016) and OptumG2 (Krabbenhoft et al., 2018), two well-established commercial software for finite element calculations in geotechnical applications. Fig. 13 depicts the results of the calculations in terms of the predicted vertical displacements directly under the centre point of the footing (point A in Fig. 4). Plaxis could not advance further than 60 kPa where the calculations stopped indicating soil collapse. Even after activating the arc-length control option, Plaxis code still diverges at that point. It is worth mentioning that Plaxis uses the modified Newton Raphson iterative scheme to solve the balance equations. A similar observation applies to OptumG2 which shows even earlier divergence around 50 kPa. The full Newton-Raphson method as implemented in Thebes code diverged around 75 kPa while the new Runge-Kutta scheme succeeded to finish the calculations smoothly. These calculations demonstrate the potential of the new method to be competitive if compared to the implicit iterative method. Table 12 and Fig. 14 report the results of the code predictions for the same footing problem but with an overburden pressure POP = 1.0 kPa to model a normally consolidated clay. The analysis involves the Newton-Raphson method and the 5th order explicit scheme. Similar observations to that in the previously discussed case of higher POP of 20 kPa (overconsolidated soil) apply here as well. The analysis run with the 5th order Runge-Kutta scheme completed successfully while the calculations with the iterative Newton Raphson procedure failed under external stress of 71 kPa. In addition, the 5th order scheme required 50 sub-increments and 300 calculations stages to go through the complete loading. The iterative scheme needed 341 iterations to reach 71 kPa before the divergence. The actual relative time needed by the Newton-Raphson method is again about 25% in excess of the time needed by the 5th order scheme to carry out a similar analysis. This example illustrates that for plasticity-dominant-problems, where the iterative procedures need an increasing number of iterations to converge, the 5th order scheme possesses, besides the convergence merits, competency with respect to the analysis efficiency and speed.

Conclusions
This paper introduced a new method to solve the load-displacement finite element equations based on Runge-Kutta scheme. The method is implemented with different local truncation errors up to 5th order. Additionally, the procedure includes a routine that automatically controls the load increment size to restrict the load path error within a desired range. The performance of the procedure is verified by solving challenging shallow foundation problems supported by an elastoplastic material. The tackled numerical applications compared the proposed method with other common methods used for the solution, including the Newton-Raphson scheme and the arc-length method. The primary findings are: 1) in comparison to Newton Raphson iterative method, the new explicit method proved to be not only more stable but also competitive with respect to the speed of solution especially when the 5th order Runge-Kutta scheme is employed; 2) the convergence of the proposed scheme is assured even in highly nonlinear problems, while the other schemes fail; 3) the presented procedure is universal, and the proposed method can be extended to even higher-order Runge-Kutta schemes easily. Finally, more research is needed to check whether the benefits of the proposed method are also present in other highly nonlinear problems. For example, those may involve coupled multi-physics where numerical instability and divergence are relatively common, for instance in the case of coupled thermo-hydro-mechanical-chemical (THMC) problems.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.