Modelling of the Torsional Simple Shear Test with Randomized Tresca Model Properties in Midas GTS NX

The dynamic properties of soil obtained from the torsional simple shear test (TOSS) are assumed to be uniform throughout the specimen. For some exceptional soils, this may hold true, but for the most majority of soils that we examine, it is obviously not the case, and this level of non-uniformity depends on the conditions in which the soil was formed. In this paper, we discuss a method of modelling inherently non-uniform soil specimens by representing them with elements that have an elasto-plastic simple Tresca material model with different properties (Elastic young modulus and yield stresses). The combination of properties that can simulate the nonlinear behaviour of the soil is found and calibrated using a model of the TOSS test built in the finite element software Midas GTS NX. Furthermore, the influence of rigid inclusions in the soil is studied and the results show an increase in stiffness with the increasing percentage of inclusion in the soil.

1 Introduction (Masing 1926) first examined the concept of nonuniform failure of a brass specimen subjected to loading and unloading in tension. He suggested that the specimen does not have uniform stiffness along its cross-section, and elements failed and reached their yield strength at different times of the test. Yield levels are mostly related to local variations in dislocation density. This type of heterogeneous dislocation is found commonly in deformed materials especially in soil. Based on that notion, a soil specimen's nonlinear behavior can be predicted and modeled by a set of elasto-plastic elements, which yield at different times of the test. A previous study by (Ahmad and Ray 2021b) on a simulation of the Torsional Simple Shear Test (TOSS) with the 3D Finite Element (FEM) software Midas GTS NX showed that a simple material model such as the Tresca model can generate complex soil behavior by varying the properties of the elements. The study examined normal, lognormal, and binomial distributions of the elastic modulus and yield stresses in the specimen, and the results showed that the nonlinearity becomes more evident when the variation of properties is on a wider range and more variable (lognormal distribution). The same FEM model mentioned before was used in a previous study (Ahmad and Ray 2021a) to prove the ability of the Ramberg-Osgood (RO) model to accurately simulate the shear stress-shear strain curves obtained from cyclic and irregular TOSS tests. This paper aims to find the distribution of Tresca elements within the soil sample, which have different properties in Midas which will generate a nonlinear behavior that matches the RO model; hence the TOSS test results.

Testing Device and FEM Model
The dynamic soil testing was conducted using a combined Resonant Column-Torsional Simple Shear device (RC-TOSS) on a hollow cylindrical soil sample with an outer diameter of 6 cm and an inner diameter of 4 cm and a height of 14 cm. The device drives the top of the specimen in torsion using a set of two magnets surrounded by four coils, whereas the bottom of the specimen is fixed. The rotation is measured by the means of two proximitors which measure the gap between the fixed sensors and the targets that are mounted on the specimen as shown in Fig. 1. Torque is calculated from the known current passing through the coils during the test. The device is able Fig. 1 Cross-section of RC-TOSS testing device to run both cyclic and irregular loading histories and plot the shear stress-shear strain curves. A detailed description of the device and the methods of calibration can be found in (Ray 1983) and (Szilvágyi 2017).
The model built in the FEM software Midas GTS NX is also a hollow cylinder with the same dimensions as in the TOSS test. It consists of 4032 high order hexahedral elements. It is pinned at the bottom while the top nodes of the top elements are rigidly linked to a node in the center where a prescribed rotation is applied to simulate the loading conditions in the laboratory tests as can be seen in Fig. 2. A time varying function is assigned to the prescribed rotation which could be a sine wave in case of cyclic loading or an irregular time history like the time histories used in the laboratory TOSS tests.
The soil used in this study is dry sand extracted from the vicinity of the Danube River near Paks in Hungary. The soil is a fluvial sediment of the river, which has a wide variation of in-situ densities. It has a very low percentage of fines and a complex grain structure with heterogenous grain shapes varying between sub-angular and rounded shapes. The properties of the sample are shown in Table 1.

Ramberg-Osgood Model
This model was developed to predict the shear modulus degradation curves for soil by Faccioli et al. (1973) for site response analysis. However, originally it was introduced by Ramberg and Osgood (1943) for aluminum-alloy stainless-steel and carbon-steel sheets. The current form that is used in this study to describe the nonlinear shear stress-shear strain behavior is as shown in Eq. (1) for the one-way curve (backbone curve): where γ and τ are the shear strain and shear stress respectively. G max is the maximum shear modulus at very small strains. τ max is the maximum shear stress, which depends on the confining stress and the angle of friction and cohesion obtained from the triaxial test. α, C and R, are constants used to fit the curve with the test data.
The original Masing criteria were later integrated to the Ramberg-Osgood model in order to simulate the two-way loading. The criteria suggest that G max is constant for all the unloading and reloading curves, and the shape of these curves is similar to the backbone curve but scaled by a factor of 2. The extended Masing criteria that were first introduced by Jennings (1964) determine the path of the curve in case of irregular loading patterns, these rules indicate that whenever the maximum shear strain is exceeded, the unloading and reloading curves should follow the backbone curve, and in case of an intersection with a previous curve, the current loading or unloading  where τ i and γ i are the shear stress and the shear strain at the turning point respectively. The parameters of the RO model that fit the curve from our TOSS test on dry sand are: G max = 95,500 kPa, τ max = 44.168 kPa, α = 1, C = 1.55, R = 1.9.

Masing Model
In the Masing model, the system consists of successively yielding elements. Each element has an elastio-plastic yield behavior, which gives departure from elastic line as more and more elements yield as shown in Fig. 3. The plastic term of the stress-strain equation in the relative coordinates is divided by two upon unloading and reloading since the yield stresses have effectively doubled. The distribution of the yield levels depends on the microstructure of the material and the dislocation arrangement in the cross section. This distribution can be described by a probability density function. The study conducted by Ahmad and Ray (2021b) on the same FEM model in this study, used the elasto-plastic Tresca model where 512 different properties were distributed on 4032 elements. The lognormal distribution of the elastic modulus and yield stresses showed the most nonlinearity in the behavior (torque-rotation curves), and the higher the standard deviation that is assigned to the distribution, the higher the nonlinearity we obtain as can be seen in Fig. 4 and Table 2. Skelton et al (1997) used a discrete distribution by approximating the continuous distribution by an isosceles triangle, where the number of yielding elements increases in linear steps as the yielding stress also increases linearly until it reaches a maximum stress  and then decrease linearly as demonstrated in Fig. 5 (before the arrows). If n is the number of yield stress levels until the apex of the triangle, then N = 2n is the total number of yield stress levels, and the total number of elements Z involved in deforming up to σ max is given by: Calculating the total stress for the entire system is different from the Iwan model, because the elements are not connected in series or in parallel, but a combination of both, thus we calculate the total shear stress for the system as follows: where the summation from 1 to j includes all those elements that remain elastic after a loading of deflection , and the summation from j + 1 to N includes all of those elements that have slipped or yielded. In case we have 100 yield stress levels (N = 100) and a total of Z = 2550 elements we chose the yield stress levels to match a maximum stress of σ max = 100 kPa. We can see from Fig. 6 that using this distribution does not give control over the nonlinearity of the curvature to obtain a good fit with the Ramberg-Osgood curve which fits the testing data. Therefore, another method is suggested to fit the two curves, and to find the combination of yield stresses for the elements that gives the fit between the Tresca model and the test data.

Calibration and Curve Fitting
In order to obtain a fit between this model of discretely distributed Tresca elements (Eq. 4) with the RO model, a method where several iterations using the solver in Excel and Midas analysis are needed. For the calibration process, only one-way static loading is applied in Midas to save time. First, as the previous example we have a constant stiffness for all the elements, knowing that the young modulus (E) is related to the shear modulus (G) and Poisson's ratio ( ) as follows: We start with 100 stress levels as the previous example, knowing that the yield stress in pure shear y is equal to one-half the Tresca yield s y . Then the solver is used in Excel (the least square method) to find the number of elements at each stress level to minimize the sum of the squares of the differences between the lab test stress values and the calculated average stress obtained from (Eq. 4) for all of the elements in the sample at each strain amplitude, with the condition that the total number of elements (Z) is 4032 equals the number of elements in the model in Midas. The number of elements for some of the stress levels is shown in Fig. 5 (after the arrows). The process of assigning the random material properties to the elements in Midas is explained in (Ahmad and Ray 2021b) and only the number of the elements and the materials are different. Briefly, after creating the model shown in Fig. 2, Midas allows the export of the mesh and its materials to a text file (.fpn) that can be converted to a (.csv) file and imported to Excel. The properties of the elements in the text file can be replaced by the properties calculated earlier by the solver in Excel after being randomized using subroutines written in VBA in Excel. Subsequently, the sheet is converted back to a text (.fpn) file and imported to Midas as a mesh with its assigned material and properties. After importing the mesh, the boundary conditions (pinned bottom) and the rotation should be assigned again to the relevant elements.
For the initial analysis, we assumed that the strain is uniform across the specimen for all the elements i , which is not the case as observed after running the analysis on Midas with these properties, the nonuniformity of the strains at the last step is shown in Fig. 7a. This led to an average shear stress-shear strain curve lower than the calculations as shown in Fig. 7b. However, this run was used to obtain an initial strain for each element i in Midas to assign in the excel sheet for the calculations in the next step.
Based on the previous yield stress level distribution, in the next step we will have 4032 yield stress levels, so each of the elements can have a different yield stress, and for the Tresca model calculation, each element will have different strain amplitudes at every step which are extracted from Midas from the prior run (Fig. 7a). As a result, every element will have different strain steps depending on its location in the sample. However, the average of the strains for all the elements will always be the same and identical to the lab test results. The solver is used again to find a new distribution of the yield stresses in the sample to match the Tresca calculated curve with the RO curve ( Eqs. 1 and 4). These properties are imported into Midas again as an (.fpn) file, and the loading conditions are defined, and the analysis is run. After every iteration, the strains and stresses of all the elements are imported to Excel to find the average strain and average stress for the entire sample at each step and compare it again with the lab test results. The results in Midas will not match exactly with the calculations in Excel but it will get closer (Fig. 8), due to the change of strain by repeating this process shown in Fig. 9 each iteration will provide a distribution of the yield stresses in the sample that will give a curve closer to the RO curve after every iteration, and when we find exactly the correct distribution of strains in the sample, we will obtain a match between Midas, RO and excel calculations after a few iterations, in this case it took 7 iterations to obtain a match.
As a result, we found the combination of Tresca properties that matches the shear stress-shear strain curve obtained from RO model hence the test data. Figure 10 shows the fit for the cyclic loading, where a time variant sine function was used to change the analysis to dynamic. This model inherently applies the two Masing criteria. In case a comparison is made by plotting the torque-rotation curves at the center point of the rigid links, it also matches perfectly with torque-rotation curve of the lab test. Figure 11 shows the distribution of the yield stress levels among the elements. This distribution can be used to regenerate Tresca properties for elements in a different Midas model, merely by knowing the number of elements in the soil layer and by using the "Random number generation" tool in Excel. Furthermore, different aspects of nonuniformity can be studied, such as voids and dense localities. Fig. 8 Iterations to reach a match between Tresca and RO models Fig. 9 Process followed to reach a match between Tresca and RO models It should be mentioned that this analysis can be conducted and compared directly with the lab test results without the need for a prior Ramberg-Osgood analysis or curve fitting. It is only demonstrated here that both models can have perfect fit with the shear stress-shear strain curves obtained from the lab TOSS test but they can be completly independent.

Effect of Randomness and Rigid Inclusions
To study the effect of the distribution of the same properties in the sample, these properties were shuffled and assigned randomly to the elements and imported to Midas where the analysis was conducted. The results of 9 different random distributions of the same properties show a maximum deviation of 4.27% from the original RO shear stress-shear strain curve.
Sample preparation method for the TOSS test is chosen depending on the required initial density of the sample. The dense state is achieved with the dry pluviation method while the loose state can be achieved by pouring dry soil into the sample mold with a glass funnel. However, uniformity in the sample is not assured, and dense localities can appear in random locations in the sample due to an error in forming the sample, the existence of particles of higher diameter and strength, or the bonding between the particles. To study this effect, rigid inclusions in the sample in the FEM model were represented by a block of elements that has a linear elastic material property (it does not yield) with a higher elastic modulus of 5 GPa which can change according to the estimated cause of the inclusion. The number of elements in this block was increased gradually starting from one block of 3 × 3 × 3 elements (0.7% of the sample are inclusions) and ending with two blocks of 13 × 13 × 3 elements (25% of the sample are inclusions) as demonstrated in Fig. 12. The shear stressshear strain curves of the models with inclusions were compared to the original curves that matched the RO model (without inclusions). The results show an increase in stiffness (G max ) with the increasing amount of inclusion in the sample as shown in Fig. 13 and Table 3, and this increase is almost linear for our tested soil, and it is given by: where y is the increase in stiffness in percent, and x is the percentage of inclusions in the sample.
(6) y = 1.205 * x + 4.23 Based on the Torsional Simple Shear tests on dry sand, and after we confirmed that the Ramberg-Osgood model is able to simulate the dynamic behavior of soils by matching the hysteretic shear stress-shear strain curves with the lab testing results. Masing's observation that the sample is not uniform, and it does not yield at the same time along its cross section is investigated. The nonlinear dynamic behavior of the sample could be generated in Midas by using simple elasto-plastic elements (Tresca model) with varied properties (elastic modulus and yield stress). A match was achieved between the lab testing data and the Masing model using Tresca model by assigning a constant elastic modulus and finding the best combination of yield stresses for the elements that has a discrete distribution. This distribution was found using an iterative method in Excel's solver and Midas. This model allows to consider the non-uniformity of soil samples and layers when running different types of analysis such as site response analysis and study the effect of dense localities where it was found that the stiffness would increase with the increasing number of rigid inclusions in the sample. However, this model has a limitation that it does not apply the two extended Masing criteria, as a result, it cannot be used for irregular loading patterns, and it only can be applied by approximating earthquake records to cyclic loading histories. Further studies can be conducted in order to surpass this limitation as well as studying the effect of voids and deformities in the sample.
Funding Open access funding provided by Széchenyi István University (SZE). This study was funded by the Stipendium Hungaricum Scholarship. All authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript. The authors have no financial or proprietary interests in any material discussed in this article.

Data Availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose. The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Fig. 13
The increase in stiffness with the increasing number of inclusion elements