Finite element simulation of the effect of loading rate on the stress-strain behaviour of Albany sand

Abu Hena Muntakim, Mohammed Saiful Alam Siddiquee and Kalum Priyanath Udagepola* Faculty of Engineering and Applied Science, Memorial University of Newfoundland, Canada. Department of Civil Engineering, College of Engineering, King Abdul Aziz University, Jeddah, Saudi Arabia. Department of Computer Science and Software Engineering, School of Information Technology and Computing, American University of Nigeria, Nigeria.


INTRODUCTION
In recent times, problems related with long-term creep deformation of sand deposits loaded with a heavy superstructure or secondary consolidation of saturated soft clay including a number of full-scale field cases, have attracted the attention of geotechnical engineers, for correctly understanding and accurately evaluating the viscous properties of geomaterials.Highly non-linear relationships of soil were the main obstacles in soil mechanics.With the development of different experimental and analytical methods, various constitutive models for defining soil behaviour have been published.The effect of strain rate on other materials has also been studied (Khan et al., 2011).
In order to simulate the effects of material viscosity on the stress-strain behaviour of geomaterial (i.e.clay, sand, gravel, and sedimentary soft rock), a set of stressstrain models within the framework of the general non-linear three-component model has been proposed by researchers (Di Benedetto et al., 2002).Three basic viscosity types have been published which are (i) Isotach, (ii) TESRA (temporary effects of strain rate and acceleration) or Viscous Evanescent and (iii) P&N (positive and negative viscosity).In this research, P&N viscosity type has been used to simulate the stressstrain behaviour of Albany sand, fine silica sand from Australia.This type of viscosity model is devised to simulate a kind of peculiar material behaviour under variable speed loading.

METHODOLOGY
The modelling of stress-strain behaviour of geomaterial is very challenging as stress-strain behaviour is highly non-linear in nature.The development of finite element analysis (FEA) has found a way to solve the boundary value problem with highly non-linear material property.
There are many commercially available FEM software today.Among them ABAQUS is the robust software that allows the user to model their own material model using user subroutine.But the challenge arises when the user wants to write their material model in another language than FORTRAN.In this study, this challenge has been successfully handled as the user subroutine for material model has been written in C++ and used.

Computational setup
ABAQUS/CAE, or 'complete abaqus environment' (a recursive acronym and backronym with an obvious root in computer-aided engineering) is used for both modelling and analysis of mechanical components and assemblies (pre-processing) and visualising the finite element analysis result.The full computational setup scheme is described here.At first, installation of the finite element software ABAQUS, FORTRAN compiler and C++ compiler is carried out.Then a change is done in windows environment variable to make FORTRAN and C++ compiler available to CMD.After that, a change in ABAQUS environment file is done to make *.lib and *.dll files available for ABAQUS.Finally, the verification procedure of ABAQUS is carried out to check whether all components are compatible with each other.
The material model code was written in C++ and it was compiled to .dll using C++ compiler.From .dll, using CMD and C++ compiler the .libfile was created.The finite element model was created using ABAQUS/ CAE.Using a FORTRAN interface the material model was called and analysis was performed.

Material model description
Di Benedetto et al. (2002) successfully simulated the rate dependent stress-strain behaviour of geomaterial observed in a number of laboratory stress-strain tests by the non-linear three-component model (Figure 1) using one-dimensional case.
Figure 2 illustrates the P&N viscosity type (Tatsuoka et al., 2008b) In this type of viscosity the viscous stress increment that developed at a given moment during subsequent loading, decays with an increase in instantaneous irreversible strain towards different residual values.The strength during monotonic loading (ML) at a constant strain rate decreases with an increase in strain rate.
In the framework (Figure 2) of the three-component model, the measured stress σ consists of two parts, which are the inviscid stress component σ f , and the viscous stress component σ v at the same Ɛ ir .Negative isotech type is a feature of σ v .Both positive TESRA type component and the negative Isotach type component at the other strain rate are the components of σ v .This can be observed when a step increase in έ at point B during ML at a constant strain rate.The stress-strain behaviour should be like A→B→D if there are only negative Isotach type components.But a behaviour like A→B→C→D, instead of A→B→D is observed in poorly graded relatively round and stiff-particle granular material.A step increase in έ (B→C) results in the same amount of immediate positive stress increase when the viscosity type is The specimens were loaded automatically.To control the cell pressure, a high precision gear-type axial loading system driven by a servo-motor together with an electric pneumatic pressure transducer was used.By increasing the effective stress from 20 kPa towards 400 kPa at an axials strain rate of 0.0625 %/min, the isotropic compression was performed.During the isotropic compression process to evaluate the vertical quasi-elastic Young's modulus, eight cycles of axial strain (double amplitude) of 0.001 − 0.003 % were applied at p = 50, 100, 200 and 300 kPa. Figure 5 shows the results from CD TC tests at different vertical strain rates on air dried Albany sand.

Pseudo-algorithm
Siddiquee et al. (2006) have developed the pseudoalgorithm, which was the revised form of the original solution technique of the DR method.Viscous effects were not included.
Isotach or combined or TESRA.After that, subsequent ML at a constant έ results in the decrease of σ v from a temporarily increased value (C→D) like the TESRA type.This feature was also found in the stress-strain behaviour of Albany sand.For this reason the P&N model is the appropriate viscosity type for simulating viscosity of Albany sand.A 0.3 mm thick latex rubber disc smeared with a 0.05 mm thick silicon grease layer (Tatsuoka et al., 1984) was used at the top and bottom ends of each specimen.In the return mapping algorithm (Ortiz & Simo, 1986), incremental elasto-plastic equations are solved at the first level of integration.Satisfying the consistency condition (abiding by the flow rule), the stress is returned to the growing yield surface.When calculating the viscous stress based on the P&N model, the stress is returned to the inviscid yield surface with an incremental integration during the second level of integration when it is necessary at each step of return mapping iteration.This scheme is presented in Figure 6.

Overall model calculation functions
In the user subroutine UMAT, the stress and hardening softening parameters are calculated from the strain and elastic modulus provided by ABAQUS.With the updated stress and hardening softening parameters ABAQUS carry out the non-linear boundary solution and provide strain and elastic modulus to UMAT.In this process the whole analysis is completed.The main function of user subroutine UMAT is UMAT_CPP.In UMAT_CPP the strain calculated by the ABAQUS solver is taken as the input and it calculates the stress in that given moment.At the end of this function, the stress is updated to ABAQUS.In this function the elastic modulus is calculated from Young's modulus and Poisson ratio.Failure surface is calculated in the function ReternMapping_.PlsticModel_ function calculates the reference curve.The function yldchk_ calculates the yield function.The calculation of invariants is done in the function invar_.The potential function and yield function is calculated in the function yieldf_.

Parameters used
Tatsuoka et al. (2008b) described various aspects of the simulation and represented simulation parameters for the stress-strain behaviour exhibiting the P&N viscosity.

framework
The present study is done using the generalised elastoplastic isotropic strain-hardening and softening model, which takes into account strain localisation associated with shear banding by introducing a characteristic width of shear band in the additive elasto-plastic decomposition of strain (Tatsuoka et al., 1993).The yield function is used as follows: Equation ( 1) is used as the growth function of the yield surface of the generalised Mohr-Coulmb type.I 1 is the stress invariant (i.e.hydrostatic stress component, positive in compression) and 2 is the second deviatoric invariant (i.e. the deviatoric stress).Siddiquee et al. (1999;2001) and Siddiquee and Tatsuoka (2001) have explained in detail about the growth function.
The plastic potential function, ѱ, is defined as; This plastic potential function of the Drucker-Prager type is similar to the yield function except that g(ϴ) in equation (1).Here in the analysis, stress dependent elastic parameters are used.

RESULTS AND DISCUSSION
In this study, experimental results of four TC tests of different strain rates have been simulated successfully.
All the stress-strain relationships depended on the rate of straining.The most innovative idea of this research was the modelling of the effect of strain rate on the shear strength of Albany sand.Usually most geomaterials increase in strength with the increase in rate of straining, but here in this case, the strength of Albany sand is reduced with the increase in rate of straining.The opposite effect of the rate of straining on Albany sand might have effects on the post-peak behaviour of the sand.It has been found that the FEM simulation matched the experimental data quite well up to the peak, and then it started to deviate.This tendency may be attributed to the non-unique nature of solution at the post peak range of the stress-strain behaviour.Further studies on the behaviour on the post-peak of Albany sand is underway.
In Figure 7 the simulated curve has been compared with the experimental data of TC test at a vertical strain rate 5.0 % / min.In this simulation the peak effective principal stress was 4.22 at irreversible shear strain 6.6 %.The simulated curve is largely deviated from the experimental curve after an irreversible shear strain of 13.4 %.

June 2016
Journal of the National Science Foundation of Sri Lanka 44(2) In Figure 8 the simulated curve has compared with the experimental data of TC test at a vertical strain rate 0.5 %/min.In this simulation the peak effective principal stress was 4.4 at an irreversible shear strain of 8.2 %.The simulated curve is largely deviated from the experimental curve after an irreversible shear strain of 5.51 %.
In Figure 9 the simulated curve has been compared with the experimental data of TC test at a vertical strain rate 0.05 %/min.In this simulation the peak effective principal In Figure 10 the simulated curve has been compared with the experimental data of TC test at a vertical strain rate 0.005 %/min.In this simulation the peak effective principal stress was 4.7 at an irreversible shear strain of 7.15 %.As this curve is accounted as base, the simulated and experimental curve is nearly the same.

CONCLUSION
In this research a visco-elasto-plastic model is used within the framework of three-component material model.The model is implemented as UMAT of commercially available software called ABAQUS.The effect of strain rate on Albany sand has been studied by the combination of elastic visco-plastic constitutive law and three component framework.The following conclusions can be drawn from this study.
TC test results of Albany sand at different strain rates have been modelled successfully into a commercially available package called ABAQUS using an unconventional approach.The P&N model was implemented into a generalised elasto-plastic isotropic strain-hardening nonlinear model in C++.The model is then embedded in the finite element computer programme ABAQUS.
The experimental data was successfully simulated up to the peak strain but the problem was that simulation deviated significantly in large deformation range.The simulated curve was deviated more or less from 7.5 % irreversible shear strain.The deviation was higher for vertical strain rate at 0.3 %/min.In this study user subroutine was written in C++ rather than FORTRAN.With the help of C++ compiler .dlland .libfiles were created and they were placed in an appropriated destination and the environmental file was updated.This made possible to call the user subroutine from FORTRAN and simulation of stress-strain behaviour of Albany sand.

From
the laboratory experiments(Tatsuoka et al., 2008a) it was found that four poorly graded granular materials, named, a) corundum A (aluminium oxide, Al 2 O 3 ), an artificial material (e max = 1.066 and e min = 0.865); b) Albany sand, a fine silica sand from Australia (e max = 0.804 and e min = 0.505); c) hime gravel, a natural fine gravel from a river bed in the Yamanashi Prefecture, Japan (e max = 0.759 and e min = 0.515); and d) Monterey No. 0 sand, a natural fine beach sand from the USA (e max = 0.860 and e min = 0.550), exhibited the P&N viscosity in the drained TC tests.In this paper, the experimental results of Albany sand are the main focus.The particle shape and size is shown in Figure3.The specific gravity was 2.64.Its maximum void ratio and minimum void ratio were e max = 0.804 and e min = 0.505, respectively.

Figure 3 :
Figure 3: Particle shape and size of Albany sand

Figure 4 :
Figure 4: Automated triaxial apparatus used in the present study

Figure 5 :Figure 6 :
Figure 5: Results from CD TC tests at different vertical strain rates on air dried dense Albany silica sand

Figure 7 :Figure 8 :
Figure 7: Experimental and simulated curve of effective principal stress, R vs irreversible shear strain at a vertical strain rate 5.0 %/min

Figure 9 :Figure 10 :
Figure 9: Experimental and simulated curve of effective principal stress, R vs irreversible shear strain at a vertical strain rate 0.05 %/min

Table 1 :
Viscosity parameters used to analyse the CD triaxial tests