Implementation of non-trivial boundary conditions in MPM for geotechnical applications

(1)

Delft University of Technology

Implementation of non-trivial boundary conditions in MPM for geotechnical applications

Remmerswaal, Guido; Bolognin, Marco; Vardon, Phil; Hicks, Michael; Rohe, Alexander

Publication date 2019

Document Version Final published version Published in

Proceedings of the Second International Conference on the Material Point Method for Modelling Soil–Water–Structure Interaction

Citation (APA)

Remmerswaal, G., Bolognin, M., Vardon, P., Hicks, M., & Rohe, A. (2019). Implementation of non-trivial boundary conditions in MPM for geotechnical applications. In Proceedings of the Second International Conference on the Material Point Method for Modelling Soil–Water–Structure Interaction: 8 – 10 January 2019, University of Cambridge, United Kingdom

Important note

To cite this publication, please use the final published version (if applicable). Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

(2)

2

Implem

1_F ABSTRACT This paper de applicability o treatments (no boundary for geotechnical has been app practical prob KEY WORD INTRODUC In continuu momentum Point Metho geotechnical Despite the often limite computation estimation o analyses and Even though al.,2016; W Element Me (usuallyfixed mechanicsor application o This paper method for conditions m (Shields, 19 which remai embankmen 2nd _Internationa

Anu

MPM Res

mentation

G.R Faculty of Civil T escribes curren of the Materia on-trivial BCs r solid materia engineering. T plied on top of blems. DS:Boundary c CTION um mechanic conservation od (MPM), l large deform abundance o ed to a stren n of factors o of the conse d displacemen h MPM has ang et al., 20 ethod (FEM) d) BCs locate r often flux-b of MPM to sl summarises n the applicati maybe used 36; Einstein, ins present du nt. al Conference o

ura3D

search Commun

n of non-t

Remmerswaa l Engineering a nt work on the al Point Metho s) in MPM for als. Both treat The flow cond f aslope leadin conditions; Geo cs, boundary n) with a set which decou mation proble of small defor ngth analysis of safety. App quence of fa nt design crit been alread 016; Coelho on an edge ed at the grid based in fluid lope stability novelbounda ion of a surf to simulate, 1950;etc.). A uring failure. n the Material P nity

trivial bou

a

al1,*, M.Bolog and Geoscienc 2_Deltares *_{E-mail: g} e implementati d (MPM) to g (i) flow condi tments are req dition BCs hav ng to failure. B otechnical eng y value prob of additional uples mass f ems(Sulsky et rmation mode s. Especially plying MPM ailure, which eria. dy applied to et al.,2018), of the backg d edges. Howe ds, i.e. they a problems, no ary treatments face traction for example A surface trac . This traction Point Method f - 61 -

undary co

pplication

gnin1, P.J. Va

ces, Delft Univ , Delft, The Ne g.remmerswaal

ion of non-triv geotechnical bo itions, as well quired, for exa ve been used to Both example gineering practi

blems usually l constraints, from the bac

t al., 1994). els, both anal

in slope st M to slope stab h, for examp o many slope these analys round grid. N ever, BCs are are decouple on-trivial BCs s in MPM fo on a moving , flow within ction can be u n could, for e for Modelling S

onditions

ns

rdon1, M.A. H versity of Tech etherlands l@tudelft.nl vial boundary oundary values

as (ii) the app ample, to estim o simulate a su s show that th ice; Material P y involve d namely bou ckground grid lytical as wel tability theor bility theory ple, is requir e stability pr ses often dire Nowadays, M e more comm d from the g s are required or flow cond g boundary ( n a slope or used to apply example, repr Soil-Water-Struc

in MPM

Hicks1, A. Ro hnology, Delft, conditions(BC s problems. It plication of a s mate the cons ubcritical flow hese new treat Point Method; S differential eq ndary condit d, has prove ll as numeric ry, current p enables a qu red for risk roblems(Ande ectly applied MPM applicat monly associa grid in MPM. d as well (Cor ditions (Zhao (Remmerswa erosional an y a load, for e resent a build Pa cture Interactio

for geote

ohe2 The Netherlan Cs) for improv summarises no surface traction sequence of sl

and the surfac tments are use Slope stability quations (e.g tions (BCs). T en to be use cal, geotechni practice is lim uantitative an assessment, ersen et al.,2 d the BCs fro

tions are ofte ated with the m

. Therefore, f rtis et al.,201 et al.,2018) aal,2017). Flo nd deposition example on to ding founded aper number: www. www. n

chnical

nds ingthe general ovel boundary n on a moving lope failure in ce traction BC eful in solving y. g. mass and The Material eful in many ical design is mited to the nd qualitative vulnerability 2010;Soga et om the Finite en limited by mass in solid for a general 8). and a novel ow boundary nal processes op of a slope, on top of an mpm2019.10 .anura3d.com .mpm2019.eu l y g n C g d l y s e e y t e y d l l y s , n

(3)

2nd International Conference on the Material Point Method for Modelling Soil-Water-Structure Interaction

- 62 -

IN/OUTFLOW BCS FOR UNIFORM AND NON-UNIFORM FLOWS

Recent advances in MPM theory have made this method a robust unified framework with efficient and versatile coupling techniques for solid-fluid-structure-gas interaction and large deformations, especially when compared to more established civil engineering methods such as Smoothed Particle Hydrodynamics (SPH) (Yan et al., 2018). To reduce the computational cost of MPM, large domains are often truncated and confined between artificial boundaries (ABCs)(Zhao & Liang, 2016; Bolognin et al., 2017; Martinelli et al., 2017). The locationof the required ABCs can be fabricated by intuition, experience, asymptotic behaviour and numerical experimentation (Papanastasiou et al., 1992). For truncation of the domain in fluid flows, ABCs are required which allow the fluid to flow in and out of the area of interest. This paper describes the implementation of an MPM BC, which improves the conditions under which material points are introduced or removed at the boundaries. The novel BC treatment thereby allows the description of well-posedproblems that are able to simulate truncated fluid flow mechanics.

To simulate a segment of a flow, the boundary must allow mass to enter andleave the domain ata predefined flow rate, and (simultaneously) control the pressure imposed from the water outside the boundary. For a solvable governing equation, all boundariesare required to only have one BC imposed. Moreover, as mass is entering and leaving the domain (in MPM), material points(MPs) must be introduced or removed (to ensure mass continuity). In flow conditions at least one BC should control the kinematics, i.e. fix the acceleration or velocity, otherwise the problem is not well-posed. However, the inflow and outflow conditions cannot simultaneously control kinematics, as physically impossible situations may then arise (e.g.not enough mass in the domain to satisfy the BCs).Many flow segments finishin either subcritical hydraulic conditions or zero pressure conditions, so a pressure boundary is more appropriate for the outflow boundary. Consequently, the inflow boundary has been selected to control the kinematics.

The novel in/outflow BCs are respectively very similar to the classic flux (velocity) and traction controlled (pressure) BCs of FEM. The difference is that due to the material discretisation into MPs, the MPs move in and out of the domain. This means specific methods for the Eularian component must be developed.

In order to impose a velocity and fixed (zero) acceleration, a zero acceleration is applied on the nodes of the boundary condition, see Figure 1a. In addition, additional material points must be added. To facilitate this an additional row of ghost elements (inflow elements) are added immediately next to the domain. Material points are added in these elements such that the water density is maintained. These points are given an initial velocity of the required flux and a fixed acceleration. Neither the velocity nor the acceleration are changed until the points enter the domain. Details of the implementation can be found in (Zhao et al., 2018).

Conversely, an outflow BC needs to impose a prescribed force (Neumann BCs), representing, for example, hydrostatic pressure (Figure 1b), in the same manner as FEM. Also in this case an additional set of elements was used in order to remove the points that leave the computational domain. As soon as material points enter these elements they are removed from the calculation.

Figure 1 Illustration of the (a) inflow BC with inflow elements (green), and (b) outflow BC with outflow elements (red) (Zhao et al., 2018).

(4)

- 63 -

In this paper a steady uniform inflow is considered for the inflow boundaries and a fixed pressure is consideredfor the outflow boundaries. However, it would be straightforward to update the procedure for a temporally varying inflow.

In order to validate the theoretical formulation and numerical implementation of in/out flow BCs, a subcriticalrectangular open channel flow problem has been simulated. An initiallyempty domain is considered of 0:15 x 0:15 x 0:01 m, in the x-, y-, z-directions, with frictionlesswalls, a prescribed horizontal inflow velocity of 1.0 m/s and a zero acceleration along the inflowboundary as indicated in Fig. 2a. The domain is discretised by line artetrahedral elements and it is modelled with only one element in thez-direction since it is an axisymmetric problem with only the x- and y-directions shown in Fig.2. A band of outflow elements is attached on the right side of the computational mesh with aprescribed hydrostatic pressure on its nodes. Outflow elements and their BCs are only activated when MPsenter the adjacent elements. The water is modelled by a simple Newtonian compressible constitutivemodel. It has an initial density ρ0 = 1 Mg/m3, dynamic viscosity µ = 1 GPa·s and bulkmodulus K = 20 MPa. The water bulk modulus was reduced by a factor of 100 from reality inorder to increase the time step size, as an explicit time integration scheme has been adopted.

Fig. 2 shows the simulation results in terms of (b) pressure and (c) velocity at different time steps. During the initial time interval of about 0.8s atransient flow is observed. The water front resembles a dam break simulation with constant water level during this period, until it reaches the wall, where in this case the material points are removed. After the initial interval, a steady uniform flow is achieved and maintained until the simulation time is complete. Note that the prescribed traction BC at the outflow boundary nodes causes unreasonably high pressures near the bottom when the flow is not yet fully developed, which is especially evident at 0.15 s and 0.20 s. This is caused by the prescribed hydrostaticpressure which is currently initiated based on the final water depth. As the calculation continues, and the outlet water level keeps rising, this inconsistency vanishes. This phenomenon can be avoided if the hydrostatic pressure at the outflow elements would be updated during the calculation based on the water depth. However, the purpose of the present simulation is to indicate the capabilities of the proposed in/outflow BCs by developing and maintaining a steady uniform flow. The BCs are capable of fulfilling this purpose as, after 0.8 s, a constant velocity, a hydrostatic pressure and a horizontal free surface are achieved.

Figure 2 (a) Computational mesh of a rectangular channel and application of in/outflow BCs for the simulation of a subcritical flow. Representation of the fluid flow converging to a steady solution enforced by the downstream BC; (b) pressure field; (c)

(5)

2 TRACTION In FEM, the always know FEM is triv body bounda the body bo background trivial to us completely application boundary co nature are u boundaries h Even though close to the MPs during approach tak nodes. This To apply the boundary on is used in FE at the bound boundaries t therefore use The edge de Method(Seth These kerne accurate rep threshold. T and the edg boundary cr points on th from polyno Once the loc the boundar compute the defined norm of the surfac (a) Figure 3I segmentation 2nd _Internationa N BOUNDA e nodes of the wn(Cortis et ial. However ary is genera oundaries as t grid, free bo se the MPs fixed BCs a of BCs at th onditions whi unaligned bo have been add h the MPs are boundary in g the entire d ken in this pa removes unn e traction con nto the nodes EM (and for a dary, wherea the location o ed to locate th etection meth hian, 1996). I el functions c presentation o This threshold ge detected b rosses an edg he backgroun omial sections cation of the e ry to the node e nodal force mal to the sur ce. The tractio

Illustration of a n of the contou al Conference o ARY CONDIT e mesh coinc al.,2018). Th r, application ally unknown they represen oundaries are or the backg at the edge o hese location ich are not al oundaries (at dressed by Co e not located the initial co deformation p aper is to app necessary inac ndition at the s. This mappi aligned BCs as for unalign of the bound he boundary hod, called th In PFM a pro can deform a of the materia d is defined i by PFM, as s ge of the back nd grid with c s of each orde edge has been es surroundin es, which ar rface, tangent on in Figures

a surface tracti ur with the bac

n the Material P TION FOR M cide with the herefore, appl n of BCs to m a priori in M nt centres of m e not guarant ground grid of a structur n is generally ligned with t least during ortis et al. (20 at the bound ondition (e.g. process, rega ply the BC ex ccuracies and e nodes, FEM ing is the sam

in MPM). Ho ned BCs all t dary must be based on the he Proximity oximity field according to al domain. T in the initial shown in Fig kground grid composite Bé er. n found, FEM ng the edge ( re directly ap

tial to the sur 3b and 3c ha (b) ion BC on a fre ckground grid, 50 kPa along Point Method f - 64 - MOVING B boundary, an lication of m moving bound MPM. Indeed mass. Moreo teed to coinc nodes for th re, can be ea y similar to the backgroun g many steps 018); here tra dary, BCs hav Wang et al., ardless of the xactly at the b d means that m M shape func me integration owever, for a the nodes su known in or information y Field Meth is created by the deforma The edge of th condition to gure 3a. Usin

are found. F ézier curves M shape funct (Figure 3c). G pplied in the rface or along s been define ee boundary: ( and (c) extern ga part of the fr for Modelling S BOUNDARIE nd the locatio most BCs to th daries in MP d, in MPM th over, due to th

cide with the he application asy to align FEM. There nd grid edges s of a simula action bounda ve generally b 2018). For s e distance to boundary, by models will n ctions are use n along the bo aligned BCs th urrounding th der to apply stored at the od (PFM), is y a summatio ations compu he domain is achieve a go ng this thresh Finally, the ed (Figure 3b). tions are used Gaussian inte MPM gove g a predefined ed to be in the (a) contour of p al loads applie free boundary. Soil-Water-Struc ES on of the edg he boundarie PM is a challe he MPs are by he fact that t background n of BCs. C along a bac efore, the foc s. Indeed, mo ation). Fixed aries are addr been applied implicity, the o the body b y distributing not have to be ed to map the oundary of th his integratio he boundary a the BC. An MPs(Remme s an impleme on of kernel f uted at the M s found by co ood fit betwe

hold the poin dge is constru The edge ca d to integrate egration along erning equatio d direction in e vertical dire (c) proximity field ed to the nodes cture Interactio ge of a contin s of a contin enge, as the y definition n the MPs flow d grid nodes. Certain BCs, ckground gri cus in this s oving bounda d displacemen ressed. to material p e BCs are app boundary. In it to the back e unnecessaril e condition fr he continuum on involves on are required. edge detectio erswaal, 2017 entation of th functions aro MPs, which a omparing the een the mater nts at which ucted by conn an thereby be the surface t g the edge ca ons. The tra ndependent fr ection. d, (b) detected due to a vertic n nuum body is nuum body in location of a not located at w through the So, it is not for example id edge. The section is on aries by their nt(unaligned) points located plied at these contrast, the kground grid ly refined. rom the body m body which nly the nodes For moving on method is 7). he Level Set und all MPs. allows for an e field with a rial boundary the material necting these e constructed traction from an be used to ction can be rom the angle

edge after cal traction of s n a t e t e e n r ) d e e d y h s g s t . n a y l e d m o e e

(6)

- 65 -

In order to demonstrate the capabilities of PFM, the effect of a surface traction on top of a soil slope has been simulated. A 10 m high, 45° clay slope, founded on a 3 m thick clay layer, is loaded under gravity. The slope and foundation layer are placed on top of a fixed boundary representing a stiff stratum. After the initial stresses have been generated an additional surface tractionq is applied over a length of 5 m. This traction can represent, for example, a small building founded on top of the slope. Two cases have been studied, in which q is applied between 5-10 m (called case 1) and 10-15 m (called case 2) from the crest. During deformation the horizontal location of q remains fixed, i.e. it is always applied between either 5-10 m or 10-15 m from the crest. However, q follows the vertical deformation of the slopeas it is applied on top of the edge computed with PFM. The traction has been fixed to be in the vertical direction. Furthermore, to investigate if PFM correctly computes the effect of different surface tractions, multiple simulations have been performed in which q has been varied, i.e.q = 0, 20, 30, 50 or 100 kPa. The soil has a Young’s modulus of 1000 kPa and a Poisson’s ratio of 0.45.An undrained cohesion ( ) softening Von Mises constitutive model has been used. The strength of the soil has been varied, with an initial cohesion of = 35, 40, 45 or 50 kPa, to study the effect of a surface traction on slopes which are either stable or unstable under their own weight. The residual cohesion is fixed to = 0.5 .

All models were setup to run a simulation of 60 s. However, in some analyses of case 1, a single material point detached from the continuum body and the simulation ended prematurely. This was due to the geometry and application of the surface traction, where a large stress concentration occurs due to the load being applied at the location where the gravitational slip surface cuts through the top of the slope. This occurred in four simulations ( = 40 or 45 kPa and q = 50 or 100 kPa). However, most of the deformation occurred before the premature end in all these simulations and the effect on the results shown in this paper is therefore small.

Figure 4 shows the MPs at the end of a simulation of case 1,with an initial cohesion of 50 kPa and an applied traction q = 100 kPa. Furthermore, the edge detected with PFM for this simulation is given in red. In general, the detected edge well representsthe boundary of the domain. However, sharp corners are rounded and small gaps can occur within the material. Smoothing techniques, such as Gaussian smoothing, can be used to reduce the occurrence of these gaps. Unfortunately, sharp corners are rounded further by Gaussian smoothing.The blacklineshows the edge detected for the same slope loaded only by gravity, i.e. q = 0 kPa, which remains stable. Finally, the green line represents the same slope loaded by the same load (q = 100 kPa), but atthe location of case 2.A slope failure is still observed for case 2, but the deformation is much smaller when compared to case 1.

Figure 4End of a simulation with = 50 kPa;MPs and their colour represent the displacement for q = 100 kPa (case 1); PFM edge for q = 100 kPa (case 1) coloured in red;PFM edge for q = 100 kPa (case 2) coloured in green;PFM edge forq = 0 kPa

coloured in black.

The maximum vertical displacement is presented in Figure 5. For q = 0 kPa, large deformations, i.e. slope failures, are observed for the two weakest materials. For increasing q the maximum deformation tends to increase as shown in Figure 5. Figure 5a show that for case 1, the maximum displacements increase with an increase ofq, because the surface traction is applied on top of the gravitational slip surface. However, for case 2 the maximum displacement for q = 20 kPaslightly decreases compared to a slope without a surface traction (see Figure 5b). The surface traction causes heaving of the soil on either side of the traction. Due to the fact that q is applied to the side of the location of the gravitational slip surface, the maximum deformation is reduced by the heaving caused by the surface traction.An increase of q above 30 kPa causes deformations beneath the tractionwhich exceed the maximum deformation caused by the gravitational failure for q = 0 kPa. Therefore, for q larger than 30 kPathe maximum deformation once again increases compared to q = 0 kPa as shown in Figure 5b.

(7)

- 66 -

For the stronger soils ( = 45 or 50 kPa) limited deformations are observed for low values of q. After qhas reach a certain value, depending on the strength of the soil, a sharp increase in the deformation is observed, i.e. a slope failure occurred. PFM can therefore be used to compute a safety limit for q,while also providing the behaviour when this limit has been reached. Decreasing the step size inq will improve the accuracy of the limit load. Finally, as expected the maximum deformations in case 1 are larger than in case 2, because the load is applied on top of the gravitational slip surface in the former case whereas it is to the side of the slip surface in the latter. Due to the adopted step size inqa difference in the safety limit cannot clearly be observed.

(a) (b)

Figure 5Maximum vertical deformation for a slope loaded by a vertical surface traction for various (kPa) and q, (a) case 1, (b) case 2.

CONCLUSIONS

Novel BCs have been outlined, which allow MPM to be more generally applied in hydro-geotechnical applications. Two numerical applications have been presented, in which these BCs are shown to give reasonable results. However, further evaluation is required to determine the accuracy of the BCs. Being able to prescribe in/outflow BCs (i) allows true steady-state conditions,(ii) reduces computational cost, (iii) simplifies the geometry, and (iv) improves the general applicability of the method.PFM has been shown to be capable of computing a limit load in MPM, while also providing the behaviour of the structure once this limit load has been reached. Both the in/outflow BCs and PFM can be further developed to increase the general applicability of MPM.

ACKNOWLEDGEMENTS

The Authors are thankful for the technical support from the MPM research group at Deltares, The Netherlands. The Anura3D software is made available by the Anura3D MPM Research Community. Marco Bolognin is grateful for the financial support by the NWO Project “MPM-FLOW: Understanding flow slides in flood defences” (Grant No. 13889). Guido Remmerswaal is part of the “Perspectief research programme All-Risk” with project number P15-21 4, which is (partly) financed by NWO Domain Applied and Engineering Sciences.

REFERENCES

Andersen S & Andersen L. (2010). Modelling of landslides with the material-point method.Computational Geosciences, 14(1), 137–147.

Bolognin M, Martinelli M, Bakker KJ & Jonkman SN.(2017). Validation of material point method for soil fluidisation analysis. Journal of Hydrodynamics, Ser. B, 29(3), 431–437.

Coelho BZ, Bolognin M, Nuttall JD, Rohe A &Aboufirass A. (2018). Assessment of dike safety within the framework of large deformation analysis with the material point method. Proceedings of Numerical Methods in Geotechnical Engineering IX, Porto, Portugal, 657–665.

Cortis M, Coombs WM, Augarde CE, Brown MJZ, Brennan A & Robinson S. (2018). Imposition of essential boundary conditions in the material point method. International Journal for Numerical Methods in Engineering, 113, 130–152.

(8)

- 67 -

States Department of Agriculture Soil Conservation Service, Washington, D.C., USA.

Martinelli M, Rohe A & Soga K. (2017).Modelingdike failure using the material point method. Procedia Engineering, 175,341–348.

Papanastasiou TC, Malamataris N&Ellwood K. (1992). A new outflow boundary condition. International Journal for

Numerical Methods in Fluids, 14(5), 587–608.

Remmerswaal G.(2017). Development and implementation of moving boundary conditions in the material point method. MSc Thesis, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, The Netherlands.

Sethian JA. (1996). Theory, algorithms, and applications of level set methods for propagating interfaces. ActaNumerica, 5, 309–395.

Shields A. (1936). Application of similarity principles and turbulence research to bed-load movement. Translation of the Ph.D. Thesis of A. Shields, Publication No. 167 Hydrodynamics Laboratory, CalTech Library, Pasadena, USA.

Soga K, Alonso E, Yerro A, Kumar K&Bandara S. (2016). Trends in large-deformation analysis of landslide mass movements with particular emphasis on the material point method. Geotechnique, 66(3), 248–273.

Sulsky D, Chen Z & Schreyer HL.(1994). A particle method for history-dependent materials.Computer Methods in Applied

Mechanics and Engineering, 118(1–2),179–196.

Wang B, Vardon PJ & Hicks MA. (2016). Investigation of retrogressive and progressive slope failure mechanisms using the material point method.Computers and Geotechnics, 78, 88–98.

Wang B, Vardon PJ & Hicks MA. (2018). Rainfall-induced slope collapse with coupled material point method. Engineering

Geology, 239,1–12.

Yan X, Li C, Chen X & Hu S. (2018). MPM simulation of interacting fluids and solids. Computer Graphics Forum, 37(8), 183–193.

Zhao X, Bolognin M, Liang D, Rohe A, & Vardon PJ (2018). Development of in/outflow boundary conditions for MPM simulation of uniform and non-uniform open channel flows. Computers & Fluids, 179, 27–33.

Zhao X & Liang D. (2016). MPM Modelling of seepage flow through embankments. In International Ocean and Polar