Co-simulation framework of discrete element method and multibody dynamics models

Abstract

Purpose

Bulk material-handling equipment development can be accelerated and is less expensive when testing of virtual prototypes can be adopted. However, often the complexity of the interaction between particulate material and handling equipment cannot be handled by a single computational solver. This paper aims to establish a framework for the development, verification and application of a co-simulation of discrete element method (DEM) and multibody dynamics (MBD).

Design/methodology/approach

The two methods have been coupled in two directions, which consists of coupling the load data on the geometry from DEM to MBD and the position data from MBD to DEM. The coupling has been validated thoroughly in several scenarios, and the stability and robustness have been investigated.

Findings

All tests clearly demonstrated that the co-simulation is successful in predicting particle–equipment interaction. Examples are provided describing the effects of a coupling that is too tight, as well as a coupling that is too loose. A guideline has been developed for achieving stable and efficient co-simulations.

Originality/value

This framework shows how to achieve realistic co-simulations of particulate material and equipment interaction of a dynamic nature.

Citation

Lommen, S., Lodewijks, G. and Schott, D.L. (2018), "Co-simulation framework of discrete element method and multibody dynamics models", Engineering Computations, Vol. 35 No. 3, pp. 1481-1499. https://doi.org/10.1108/EC-07-2017-0246

1. Introduction

In the handling of bulk materials, such as iron ore, coal and woodchips, equipment and particles interact with each other. Their complex interaction has discouraged the advancement of analytical models for dynamic interaction problems, and development of handling equipment is mostly done by trial and error of physical prototypes. Replacing this part of the design process with inexpensive virtual prototypes would reduce the costs and time required to complete a design iteration. However, the required single computational method that can capture both particle motions and equipment dynamics does not exist; therefore, separate solvers have to be used instead.

One method to simulate the behaviour of handling equipment is the multibody dynamics (MBD) simulation. This method numerically simulates systems composed of multiple bodies each having mass, inertia and degrees of freedom (Whittaker, 1970; Wittenburg, 2007; Meijers, 1997). The bodies are connected with each other by means of joints, cables, contacts or other kinematic or force constraints. The bodies and constraints lead to the equations of motion of the system which can then be solved. Overall, the MBD method has proven to be a useful tool for motion analysis of multibody systems (Langerholc et al., 2012).

The discrete element method (DEM) is a particle-based method that can be used to simulate the behaviour of a bulk material. The method computes the individual behaviour of each particle, studying its interactions with neighbouring particles and walls (Cundall and Strack, 1979). By calculating the interaction forces, the resulting motion can be computed with the help of the equations of motion. After collecting all the information of the particles, it is possible to study the behaviour and flow of a bulk material. This is a computationally intensive method, which has only recently become applicable for large scale problems with the recent increase in computational power.

A combination of both methods could capture both material behaviour and equipment behaviour and predict the equipment performance as presented in Figure 1 (Coetzee et al., 2010). The DEM could compute the loads from the bulk material on the equipment geometry and feed these values to the MBD. This method takes these loads and calculates the corresponding movements of the geometries. These movements are sent back to the DEM program which then can start computing behaviour of the discrete elements and the loads on the geometries again.

By repeating those steps repeatedly over time, both equipment and bulk material behaviour are captured by the co-simulation. This approach would allow for investigating and improving equipment designs quickly and inexpensively. However, a comprehensible approach for coupling these two computational methods was not found.

This paper establishes a framework for the development, verification and application of a co-simulation of DEM and MBD. To that, first, we investigate how the DEM and MBD can be coupled, which is a requisite for coupling the equipment model with the material models. Then we establish a method for verifying the coupling through series of simple tests with known analytical solutions, proving the correctness of the coupling program. Finally, we examine the robustness of the coupling, presenting a generic applicable guideline for obtaining stable results.

Section 2 documents the coupling for both directions, from DEM to MBD and from MBD to DEM, Section 3 demonstrates that the proposed method of coupling is capable of achieving accurate results by examining a selection of scenarios of coupled models where the co-simulation is compared to an analytical solution. Section 4 examines the stability of the co-simulation and offers an approach to selecting a feasible timestep, as well as a suitable communication interval for robustness.

2. Coupling of an equipment model and a material model

For a co-simulation of both particulate material and equipment, three software components are required as displayed in Figure 1. The first component is the DEM for simulation of particle behaviour and forces acting on the equipment. These loads are sent to the MBD component which consists of a mechanism of multiple bodies connected to each other through cables, joints and contacts. The load on the mechanism affects the position of its bodies and these new positions are sent back to the DEM program. The coupling consists of sending the position of the selected bodies from the dynamic model to the particle model and sending the load data of the selected bodies from the particle model to the dynamic model. A coupling server, the third software component, is used to communicate with the two simulation packages, sending input to them and asking for output from them.

The DEM software component used in this research is EDEM®, a package developed by DEM Solutions (2014a). The program can be accessed through the graphical user interface or through the extended application programming interface (API). As the data exchange needs to occur at regular intervals, a coupling without user intervention through the API is highly preferable, as this will eliminate the need for the user to manually couple the programs. By writing C++ user defined libraries, users can add their own programming to the DEM software. The API offers possibilities for user-defined contact models, particle body forces, particle generators and MBD coupling.

The MBD software component used in this research is Adams®, a package developed by MSC Software (2014). The program can collaborate with other programs through the Adams/Controls package. Previously, the only programs able to connect with the Adams®/Controls package were MSC Software’s EASY5® and MATLAB® from the MathWorks Inc., but recently this list has expanded to the Adams® external interface and the Functional Mock-up Language. The Adams® external interface offers programmers access to Adams®/Solver through a collection of functions in C and exchange information when desired.

Published research on coupling DEM and MBD programs has been limited so far. The studies by Yoon et al. (2011, 2012) and Park et al. (2012) investigated the inclusion of particles in a MBD package, focussing on the use of the computer’s graphical processing unit to perform calculations. Balzer et al. (2016) and Hess et al. (2016) used an open source coupling strategy through Functional Mock-up Interface. Edwards et al. (2013) and (Barrios and Tavares, 2016) used a coupling specific for high pressure roll grinding. A general purpose approach for coupling with Adams® has been written by Elliott (2000), and offers several handles how to accomplish a robust and accurate coupling with an external interface. Coetzee et al. (2010) modelled a dragline bucket where the DEM code PFC uses an external component that calculates the motion of the dragline bucket. A coupling written for EDEM® and Adams® by ESTEQ (2008) has been used in the earlier stages of this research. This coupling uses MSC. Software’s EASY5® to communicate with Adams® and is restricted to outdated versions of EDEM. This made it impossible to use recent versions of the EDEM® program and to simulate batches of co-simulations. Therefore, a light coupling server has been developed that could interact with the recent versions of both packages without the need to configure EASY5® for each co-simulation.

This coupling server software component is a self-written program which interacts with the DEM and MBD components. The coupling server couples the selected MBD bodies and the selected DEM geometries, sticking the two models together and behaving as glue between the two models. As both DEM and MBD do not simulate at the same speed and do not advance with the same timestep, the coupling server has to coordinate the exchange of data and the schedule required. For example, a stand-alone DEM simulation advances with a timestep of 1e-6 s, while a stand-alone MBD simulation prefers to simulate with a timestep of 1e-2 s or more. Running both the codes at the same timestep is undesirable, as this will increase the computational costs for the faster simulation and in addition increase computational costs because of the overhead of the coupling. Instead, the two simulations communicate at a communication interval nΔt. This is specified in the coupling and is for practical reasons a multiple n of the smallest timestep Δt specified in DEM. The coupling server works for both the directions, coupling the load data from DEM to MBD and the position data from MBD to DEM, and this is explained in the following subsections.

2.1 From DEM to MBD

The main purpose of linking the bulk material simulation to the equipment simulation is to perform the equipment simulation with loads caused by the bulk material. Starting point is a vector with the position and trajectory x(t) of a coupled body during a communication interval nΔt. This vector contains the position, orientation and velocities of the geometry. With this data of the geometry, the DEM simulation checks whether particles are in contact with the body during the communication interval and calculates the corresponding load data F(t) in equation (1):

(1) x ( t ) → DEM → F ( t )

The load data F of each coupled body consists of six components: both forces F and torques T in three directions. The forces and moments are oriented in the directions of the global coordinate system and around the body’s centre of mass. At each communication interval nΔ t, the load data F is requested from the DEM simulation for all the bodies to be coupled.

However, the load data F(t) cannot be applied directly in the MBD simulation, as this value is not necessarily constant during the communication interval nΔt. If the load calculated in DEM increases every timestep Δ t as shown in Figure 2, the correct value to be communicated with MBD is denoted by the grey area. Unfortunately, the dynamics calculation uses the communicated value without regard for its history during the communication interval, which in this case leads to an overestimation displayed by the Error area. The communicated value should therefore be somewhere between F(t − nΔt) and F(t).

This can be resolved by averaging the force during the last n timesteps that took place since the previous communication time in equation (2):

(2) F¯(t)=1n∑F(t)+F(t−Δt)+F(t−2Δt)+…+F(t−nΔt)

where F¯(t) is the average force at time t during the communication interval nΔt. This average force is computed by starting at zero and at each timestep adding the particular force data. This operation is performed for all coupled bodies and all of their components according to equation (3):

(3) F(t)→Force averaging→F¯(t)

The Adams® solver requests the inputs to the solver to be given at the same time as the outputs. This means that if the load data F¯ of time t is fed to the solver, it would return the accompanying positions and velocities of the bodies at time t. Unfortunately, the positions and velocities at time t were the starting point of equation (1). In fact, at the next communication between MBD and DEM the position data at time t + nΔt is required. In order to have the MBD solver deliver this, it requires load data at t + nΔ t.

To achieve the load data F¯(t+nΔt), the available load data F¯(t) have to be extrapolated (Elliott, 2000) in equation (4):

(4) F¯(t)→Force extrapolation→F¯(t+nΔt)

The load data are extrapolated to t + nΔt by finding a quadratic solution based on the last three available load data points as displayed in Figure 3. First step is finding a quadratic solution displayed in equation (5) which fits the three data points:

(5) F¯(t)=c1t2+c2t+c3

The solution is found by applying Cramer’s Rule (Cramer, 1750) to the system of linear equations in equation (6).

(6) TCc=F¯C=[t2t1(t−nΔt)2t−1nΔt1(t−2nΔt)2t−2nΔt1][c1c2c3]=[F¯(t)F¯(t−nΔt)F¯(t−2nΔt)]

According to Cramer, the three coefficients c₁, c₂ and c₃ can be defined [equation (7)] in terms of determinants defined in equations (8), (9), (10) and (11).

(7) c1=det(C1)det(TC),c2=det(C2)det(TC),c3=det(C3)det(TC)

(8) det(C1)=|f¯(t)t1f¯(t−nΔt)t−nΔt1f¯(t−2nΔt)t−2nΔt1|

(9) det(C2)=|t2f¯(t)1(t−nΔt)2f¯(t−nΔt)1(t−2nΔt)2f¯(t−2nΔt)1|

(10) det(C3)=|t2tf¯(t)(t−nΔt)2t−nΔtf¯(t−nΔt)(t−2nΔt)2t−2nΔtf¯(t−2nΔt)|

(11) det(TC)=|t2t1(t−nΔt)2t−nΔt1(t−2nΔt)2t−2nΔt1|

The three coefficients of equation (7) can then be applied in equation (5) to extrapolate the load data F¯ to t + nΔ t as shown in equation (12):

(12) F¯(t+nΔt)=c1(t+nΔt)2+c2(t+nΔt)+c3

This operation is performed for all force and moment components for each body at each communication interval.

Another benefit of the extrapolation of load data is that it works similar to a smoothing operation. This will prevent the DEM simulation of feeding erratic forces to the MBD solver which could jeopardize the stability of the co-simulation. The stability of a co-simulation will be discussed in more detail in Section 4.

2.2 From MBD to DEM

The MBD solver calculates the movement of the mechanism based on the external loads from DEM and the mechanism modelled in Adams® [equation (13)]. The movements of the coupled bodies from the MBD solver are required for the next iteration of the coupling in in EDEM®.

(13) F(t+nΔt)→MBD→x(t+nΔt)

In a way, the MBD solver is always ahead of the discrete element simulation, as it uses extrapolated loads to compute the position data. Before the outputs from the MBD solver can be applied in the DEM package, coordinate transformation is necessary [equation (14)]. Two differences between the codes have to be resolved for a successful co-simulation: position versus translation and Euler angles versus orientation matrix:

(14) x(t+nΔt)→Coordinate  transformation→x(t+nΔt)

Adams® uses the absolute position of a body x, while EDEM® uses the translation Δ x of the body with the original position x₀ of the body as origin. These differences can simply be overcome by applying equation (15):

(15) Δx=x−x0

Adams® uses a 3-1-3 body rotation sequence (Kane et al., 1983), and the transformation from Euler angles ψ, θ and ϕ to a rotation matrix O can be achieved by equation (16):

(16) O=[cos(ψ)cos(ϕ)−sin(ψ)cos(θ)sin(ϕ)−sin(ψ)cos(ϕ)−cos(ψ)cos(θ)sin(ϕ)sin(θ)sin(ϕ)cos(ψ)sin(ϕ)+sin(ψ)cos(θ)cos(ϕ)−sin(ψ)sin(ϕ)+cos(ψ)cos(θ)cos(ϕ)−sin(θ)cos(ϕ)sin(ψ)sin(θ)sin(ψ)sin(θ)cos(ψ)sin(θ)cos(θ)]

By default, both programs define the orientation around the current centre of mass.

At this point, the geometry data at t + nΔt is sent to DEM program. This includes the position, velocity, orientation and angular velocity of each coupled body. While the position and orientation are needed for correct positioning of the coupled body, this is not true for the velocities. These are used for the calculation of contact forces according to the definition of the Hertz Mindlin contact model (Mindlin and Deresiewicz, 1953).

The EDEM® software automatically interpolates the geometry data for all the DEM timesteps between t and t + nΔt according to the linear interpolation of the SLERP algorithm (Shoemake, 1985). The difference between the non-linear movement and the linear interpolation of the SLERP algorithm should be considered when selecting a suitable communication time nΔt. For example, a circular motion during the time nΔt will be interpreted as a straight line; therefore, n should be chosen in such a way that this effect can be ignored.

3. Verification of the coupling

The coupling of the two codes discussed in the previous section is verified in this section. The verification process is performed with four different tests, starting with a simple test and gradually increasing the complexity. At each simulation experiment, the simulation results are compared to the analytical solution. Without such a verification, it will remain unclear whether the implemented coupling can produce accurate results.

The four tests have been selected to test the implemented coupling step by step, with each test focussing on additional aspects:

1. The first test examines whether the implemented coupling can accurately predict the contact of a single particle colliding with a wall of a coupled body. Also, this test examines if the coupling server is interfacing correctly with Adams®. Adams® is given a known load of particle collisions and has to return the corresponding behaviour of the body. This test is limited to translations only.

2. Aim of the second test is to verify that the coupling server interfaces correctly with EDEM®. EDEM® will be given a known motion for a body and has to return the corresponding load data to the coupling server. The given motion consists of both translations and rotations.

3. The third test verifies that the two-way coupling works for a scenario limited to translational forces and movements.

4. The final test of the coupling of DEM and MBD method verifies that the two-way coupling works for a scenario including rotations and moments.

The agreement between the co-simulation results and the analytical results is evaluated according to the coefficient of determination R² described by Weisberg (2005). Because both the results are based on the same mathematical problem, a very high coefficient of determination is to be expected. Correlation should exceed 0.99, demanding high accuracy while allowing for numerical scatter in the co-simulation.

3.1 Particle–wall collision

The first test of the verification process is a simple test where particles are shot at and collide with a geometry. Instead of using two particles in DEM, one of them has been replaced with a body whose behaviour is calculated with MBD. Aim of this test is to investigate whether a simple collision between particle and geometry can be computed correctly. In this first scenario, displayed in Figure 4, a single particle was generated and collided with the cube, which only has translational degrees of freedom.

The particle has a radius of 10 mm, a coefficient of restitution of 0.5 and an initial velocity of u_p = 2 m/s. The cube was at rest at the start of the scenario and has a mass m_c of 10 kg. Different particle masses m_p were used in this test, ranging from 0.01 kg to 10 kg. During the simulation, where gravitational forces are absent, the particles collided with the cube, transferring their energy according to equations (17) and (18):

(17) vc=mpup+mcuc+mpCR(up−uc)mp+mc

(18) vp=mpup+mcuc+mcCR(uc−up)mp+mc

where v_p is the velocity after collision, u_p the initial velocity and m_p the mass of the particle and v_c, u_c and m_c for the cube. Simulations have been performed with a timestep of 5e-6 s.

Table I shows the velocities of the colliding bodies after the collision, as well as the error according to equation (19).

(19) Error=100(vSimulationvTheory−1)

It can be observed that the test particles of 100 times lighter showed acceptable agreement. However, simulations using heavier particles with a particle/cube ratio of 5 and 1 did not show good agreement with theory. This is presumably caused by the calculation of the particle/wall contact in EDEM®. The built-in contact models all assume that the mass of the wall or body is equal to 1e8 kg (Arumugan, 2014), while in fact it should be significantly less. This affects the value of the equivalent mass which in turn affects the calculation of the damping force. To reduce this error to an acceptable level, users should limit the geometry to particle mass ratios to at least 100 to 1 or use a customized contact model that does not assume the mass of the geometry. Simulations with high coefficients of restitution C_R will be affected less as the damping component of the contact force is smaller.

Next, the single particle is replaced by a stream of particles, each with a mass of 0.01 kg and hitting the cube with their own initial speed of 2 m/s minus the current speed of the cube. Particles are created at a rate of 200 particles per minute. Because of the incoming momentum of the particles, the cube starts accelerating, although the rate of acceleration decreases because of the increasing distance between the starting point of the particles and the position of the cube, requiring the particles to travel longer before they collide with the cube. Figure 5 shows the velocity and displacement of the cube calculated by the co-simulation and theoretical calculations based on equation (17). It can be observed that the agreement between co-simulation and theory is very good as the coefficients of determination R² exceed 0.999. This proves that the collision between particle and body is computed correctly.

3.2 Motorized rotating pendulum

The second test of the verification process used a box filled with particles connected to a pendulum, rotating at a constant angular velocity. Aim of this test is to assure that the coupling server correctly updates the position of the body in EDEM® and whether the load calculated in EDEM® corresponds with the theoretical approach. This will be achieved by comparing the forces acting on the joint of the pendulum predicted by the simulation to the gravitational and centrifugal forces of the box with particles. As the interaction forces do not affect the motion of pendulum, the forces on the joint and the torque required for the rotation can be calculated by equations (20), (21) and (22):

(20) Fx=−mω2rcosθ

(21) Fy=mω2rsinθ+mg

(22) Tz=Fxrsin θ−Fyrcosθ

(23) m=Σmp+mc

In this test, the pendulum has a mass m according to equation (23), consisting of Σm_p of 236.6 kg and m_c of 10 kg. The pendulum has a horizontal starting position and rotates counter-clockwise at a constant angular velocity ω of 72 deg/s. The radius r of the pendulum is 1.5 m, and gravity is acting at 9.81 m/s². Figure 6(a) shows the simulation of the motorized pendulum.

Figure 6(b) shows the torque required for rotating the particle box at the prescribed angular velocity. First, the box in the simulation needs to be accelerated to the prescribed angular velocity which results in a mismatch during the first 0.25 seconds as equations (20) and (21) assume constant velocity. For the remainder of test, it can be seen that the predictions of the co-simulation match very well with the theoretical approach based on equation (22), resulting in a coefficient of determination of 0.999. This confirms that the coupling server and EDEM® calculate the expected load output well while given a known input, both for translational movements as well as rotational movements.

3.3 Translating spring damper system

The third test of the verification process uses a box of particles connected to a translational spring damper system (Figure 7). Aim of this test is to investigate whether the coupling server and MBD process the load input correctly and if it results in the desired response of the system. This test is again limited to translational forces and movements. The response of the spring damper system can be described according to equations (23) and (24):

(24) md2xdt2+cddxdt+kx=mg

The box has a mass of 100 kg and the particles have a mass of 261.8 kg, while the spring has a constant k of 5,000 N/m, preloaded at 100 N and the damper has a damping coefficient c_d of 500 Ns/m. The particles have a coefficient of restitution C_R of 1 to exclude damping from the particles.

Figure 8 shows the velocity and position of the box of particles during the simulation. Because of weight of the particles and the box, the box drops down in search for a new equilibrium. Because of the damper present in the system, the velocity is reduced, damping the system and gradually slowing down the box. It can also be observed from Figure 8 that the response predicted by the co-simulation correlates well with the response following from equation (24), resulting in a determination coefficient of 0.998. This proves that the coupling server and MBD software Adams® properly transform the translational forces of the particles to translational movements of the mechanism.

3.4 Torsional spring damper system

The final test of the coupling of DEM and MBD method is performed by confirming the response of a torsional spring and damper system. Aim of this test is to verify that the MBD software Adams® and DEM software EDEM® are coupled correctly to the coupling server, both for translations and rotations. The motor of the pendulum example in Section 3.2 has been replaced by a torsional spring with a stiffness κ of 40 Nm/deg and a torsional damper C_d of 20 Nms/deg. The system can then by described by equations (23) and (25):

(25) d2θdt2I+dθdtCd+κθ=mgrcosθ

where I is the moment of inertia and the load on the system is determined by the mass of the box m_c = 50 kg, the total particle mass Σm_p = 266.4 kg and the angle of the pendulum θ. The initial position of the pendulum is horizontal.

In Figure 9, the angular velocity and rotation of the pendulum are shown. Due to the weight of the box with particles, the pendulum started accelerating. The torsional spring prevents the pendulum from reaching vertical position, whereas the torsional damper reduces the angular velocity. When the outcome of the simulation is compared to the equation of motion it becomes clear that the coupling of DEM and MBD produces very accurate results with a determination coefficient of 0.999.

This final test concludes the verification process of the coupling method described in Section 2. This proves that the developed coupling of DEM and MBD works as expected and is capable of accurately simulating systems where particles and mechanisms interact.

4. Coupling stability

Besides an accurate co-simulation of DEM and MBD, a robust coupling is also desired. Without robustness, simulations will fail or produce considerable errors. Both the solvers of the DEM and MBD have their own preferences for achieving a stable simulation with the lowest possible computational costs. However, when the two solvers are coupled in a co-simulation, these preferences sometimes conflict and need to be resolved. This section focuses on the stability of a co-simulation and aims to provide users with a guideline for robust and stable co-simulations while minimizing computational costs.

4.1 Stability of DEM

The stability of a DEM simulation is determined by the size of the timestep Δ t or Δ t_DEM, which is the stepsize the simulation uses to advance through time. It is essential that these timesteps are not too large, as contacts need to be detected in time in order to calculate the interaction forces correctly. If the selected timestep is too large, overlap between the particles might be aggravated, causing a disproportional response from the contact model. This response consists of interaction forces which are too large and leads to a large acceleration of the particle, therefore resulting in a large displacement during a timestep. It is likely that the next contacts of the particle will also be disproportional, initiating a chain reaction of unstable contacts which will result in the explosions of particles, a scenario well known by DEM users.

Choosing a stable timestep size is a topic of interest for many DEM users, as a very conservative estimate of the timestep will raise computational costs significantly. One of the available guidelines of choosing a suitable timestep is the definition of the Rayleigh timestep shown in equation (26) (Ning, 1995). This is the amount of time required for the propagation of a surface wave through a particle and has proven to be a useful tool for estimating a suitable timestep.

(26) ΔtR=πRρG1(0.1631ν+0.8766)

Users are recommended to take a fraction γ_R of the Rayleigh timestep Δt_R of 20per cent for systems with high coordination numbers (=¿4) and up to 40per cent for lower coordination numbers (DEM Solutions, 2014b).

Another guideline for determining the critical timestep is based on the eigenfrequency ω of a single particle (O’Sullivan and Bray, 2004) shown in equation (27). Here, a safety factor γ_ω of 0.8 is commonly considered to ensure a stable integration.

(27) Δtω=2mkn

The value for Δt_ω is usually smaller than Rayleigh critical timestep Δ t_R, as pointed out by Yade-DEM (2015). For example, the critical timestep Δt_ω for the particles of the material model is 20 per cent of Δ t_R, while for the large-scale model this ratio increases to 31 per cent. It can be concluded that the critical timestep based on the eigenfrequency Δt_ω is slightly more conservative than the Rayleigh critical timestep Δt_R when the safety factors γ_R and γ_ω are considered. In case of simulation scenarios where particles collide at high velocities, smaller timesteps should be taken to avoid disproportional response of the contact model.

4.2 Stability of MBD

The solver of a MBD simulation works completely different from the DEM simulation. As opposed to solving the governing equations in discrete timesteps in the DEM, the Adams®/Solver solves governing equations of the model in continuous time. Adams® default solver GSTIFF uses backward differentiation formulas and fixed coefficients for prediction and correction, based on the work of Gear (1971). It is a variable step size integrator, which means that the internal MBD simulation time does not advance at constant intervals, as it actually can slow down and reverse if the corrector has trouble converging.

Stability of the results is guaranteed by the corrector of the solver that monitors the error in the solution and checks if this is smaller than the specified corrector error tolerance. The corrector can force the solver to reduce its stepsize or repeat the previous steps at a smaller interval. It can be helpful to use the modified corrector instead of the original in models that have discontinuities in their force functions, for example when interacting with DEM. The modified correct only applies error control on a limited set of variables such as the displacements and is less strict on the prediction errors of the load data from DEM. The complete implementation of the corrector can be found in the Solver’s manual (MSC Software, 2013).

Other variants of the GSTIFF solver are the WSTIFF solver (van Bokhoven, 1975), which uses variable coefficients instead of the GSTIFF’s constant coefficients based on the assumption that the timestep does not change. Each time the timestep changes, the GSTIFF solver introduces a small error, while the WSTIFF solver prevents this. This makes WSTIFF is a more suitable solver for simulations with discontinuous forces such as contacts or interaction with other discrete functions such as the DEM.

4.3 Stability of co-simulation

The two different solvers’ approaches cause some challenges when it comes to combining MBD and DEM. Elliott (2000) demonstrates several examples where the co-simulation results are not computed accurately, because of combination of a continuous MBD solver and an external discrete solver. The MBD solver selects its own timestep, and only requires the error tolerance to be configured. However, each time a communication takes place, a timestep is forced in the MBD solver. The maximum MBD timestep Δ t_MBD therefore depends on the DEM timestep Δ t_DEM and the communication interval n [equation (28)]:

(28) ΔtMBD≤nΔtDEM

The stability of the coupling is investigated by examining co-simulations of the spring-damper system from Section 3.3. Table II presents the accuracy and costs of co-simulations using different DEM timesteps Δ t_DEM and communication intervals n. It shows three coefficients of determination: for the position x, velocity v and force F acting on the geometry. Most important is the accuracy of the position of the geometry, as an inaccurate position can lead to missed contacts and false contact forces. Accurate simulations have been marked grey when the correlation between the simulation and the analytic solution is above 0.99.

When the interval is chosen too large, the computed velocity starts to deviate from the analytical solutions. This is shown in Figure 10 when Δt_DEM = 1e−4 and n = 25. This is undesirable, as the geometry velocity is part of the contact interaction force calculation and therefore an error can propagate through the system. The determination coefficient of the load data on the geometry is less high compared to the position and the velocity of the geometry. This is mainly caused by the discrete nature of the particle contacts and often compensated at the next communication interval. For example, when the interaction force is exaggerated this triggers larger than accurate acceleration of the geometry, while at the next communication the contact force is likely to be underestimated due to the aggravated displacement of the geometry during the interval. However, this behaviour tends to escalate into instability the more iterative steps are taken, as the coupling is in fact too loose.

When the interval is chosen too small and the MBD solver is forced to have a very small timestep Δt_MBD, the accuracy of the solution is affected as well [equation (28)]. Figure 11 shows an example where a small DEM timestep Δt_DEM = 1e–5s is chosen in combination with a communication interval of n = 1, here the computed solution starts to differ from the analytical solution. By forcing the MBD simulations to use a very small timestep, numerical scatter is amplified, resulting in the so-called pinging of the mechanism. The pinging can be observed from the erratic curve of the velocity and position of the geometry in Figure 11. The proposed force extrapolation of Section 2.1 reduces this problem, as it also acts as a smoothing operator, although it obviously does not eliminate this problem.

The computational costs of a co-simulation also depend on the selection of the communication interval n. The critical timestep of the MBD solver is usually much larger than the DEM. Simply connecting both software components with n = 1 can be undesirable, as this requires both the components to run at the same speed and may result in additional computational costs while accuracy does not improve. It can be observed from Table II that co-simulations with a communication interval n = 1 dramatically increases computational costs with a factor of up to eight, because at each DEM timestep the two software components communicate with each other. As the communication interval increases, the computational costs decrease because the MBD solver is not forced to use very small increments. Simulations with a DEM timestep of 1e-4 s and an interval of 100 or more have even lower computational costs; however, their low costs are because of instability of the co-simulation, causing all the particles to leave the computational domain which affects the computational costs dramatically.

A guideline for achieving a stable and efficient co-simulation is shown in Figure 12 and consists of the following steps:

• For the DEM, a stable timestep has to be determined, based on the Rayleigh critical timestep Δt_R or the eigenfrequency critical timestep Δt_ω. This needs to consider particle masses, stiffness and impact velocities characteristic for the simulation. A safety factor γ needs to be considered; however, a very conservative factor will lead to an unnecessary increase in computational costs.

• The MBD solver needs to be chosen based on the presence of contacts in the simulation and the nature of the interaction between particles and equipment. If the interaction is intermittent or contacts are causing abrupt events in the simulation, the WSTIFF solver is recommended. The corrector and its error tolerance need to reflect the nature of the simulation. If the simulated system contains discontinuities such as contacts that are high impact collisions between geometry and particles, the modified corrector can be used, as it is less strict on force prediction errors.

• The communication interval n needs to be chosen in such a way that the DEM code is provided with sufficient updates of the geometry’s position in order to prevent a disproportional response from the contact model. A value of n = 5 can be selected as initial value. For simulations with a small Δ t_DEM and a large geometry mass a higher value for n can be helpful to prevent forcing the MBD solver to use small timesteps, reducing the risk of pinging.

• The stability of the co-simulation can be tested by experimenting with the communication interval n. By examining co-simulations with different intervals, the quality of the results can be assessed. A communication interval that is too large can be recognized by comparing the velocity profiles of the geometries to the derivative of the position data. An interval that is too small can be identified by examining the forces acting on the geometry. If these forces are extremely erratic, it is likely that the MBD will have difficulties in successfully computing the accompanying velocities and displacements.

• If the co-simulation is found to lack stability, its settings need to be reconfigured. A first step in achieving a more robust solution is to alter the communication interval, depending on the evaluation. If changing the communication interval does not produce the desired effect, the MBD solver needs to be reconfigured. As a last resort, as this has the largest impact on the computational costs, the DEM timestep can be lowered.

5. Conclusions

This paper shows how a DEM model can be successfully coupled with a MBD model into a co-simulation. By exchanging load and position data, both models can cooperate and together compute the interaction between bulk material and handling equipment. The two computational methods have been coupled in two ways, which consists of coupling the load data on the geometry from DEM to MBD and the position data from MBD to DEM.

The coupling has been tested thoroughly in several scenarios, starting with a simple scenario of a single collision between particle and geometry and concluding with a complex scenario combining translation and rotation. All tests clearly demonstrated that the coupling is successful in predicting particle-equipment interaction.

A guideline has been developed for achieving a stable and efficient co-simulation. The robustness of the coupling has been assessed, demonstrating cases where the coupling is too tight, as well as the effects of a coupling that is too loose. When the proposed guideline for a stable co-simulation is adopted, this coupling technique is ready for simulating large scale material–equipment interactions.

Future work will combine this developed coupling framework with a MBD model of real scale industrial equipment and large-scale DEM material models into a co-simulation. How accurate this co-simulation performs in predicting equipment performance will be investigated through on-site validation of a grab with iron ore pellets.

Figures

Figure 1.

A Co-simulation using MBD and DEM

Figure 2.

Force averaging in DEM

Figure 3.

Force extrapolation required for MBD solver

Figure 4.

Particle cube collision test

Figure 5.

A stream of particles hitting the cube, causing it to accelerate

Figure 6.

Motorized pendulum

Figure 7.

Schematic of translational spring damper system

Figure 8.

System response predicted by the co-simulation compared to the system response according to equation (24)

Figure 9.

Response of system calculated by a coupled simulation compared to the equation of motion in equation (25)

Figure 10.

Example of unstable coupling: the communication interval is too large Δ t_DEM = 1e–4 s

Figure 11.

Another example of unstable coupling: pinging of the MBD solver because the communication timestep is too small Δt_DEM = 1e-5s

Figure 12.

Achieving a stable co-simulation

Table I.

Particle–geometry collisions with u_p = 2 m/s and m_c = 10 kg

mp	mcmp	Theory		Simulation		Error (%)
		vp	vc	vp	vc	p	c
0.01 kg	1,000	−9.97e-1	3.00e-3	−9.98e-1	3.02e-3	0.141	0.672
0.02 kg	500	−9.94e-1	5.99e-3	−9.95e-1	6.03e-3	0.130	0.642
0.1 kg	100	−9.70e-1	2.97e-2	−9.74e-2	2.94e-2	0.365	0.783
1 kg	10	−7.27e-1	2.72e-1	−7.56e-1	2.75e-1	4.31	0.943
2 kg	5	−5.00e-1	5.00e-1	−5.54e-1	−5.10e-1	10.8	2.09
10 kg	1	5.00e-1	1.50e0	3.87e-1	1.62e0	22.7	7.84

Table II.

Effect of different timestep and communication interval Δ t_R = 5.21e–4s

ts = ΔtDEM(s)	ΔtMBD(s)	n=ΔtMBDΔtDEM(–)	ttotal (h: min: s)	R2 of x	R2 of v	R2 of F
1e-5	1e-5	1	4:03:43	0.09660	0.58219	0.00185
	5e-5	5	0:53:54	0.85019	0.96361	0.02700
	1e-4	10	0:36:29	0.99985	0.99995	0.68607
	2.5e-4	25	0:28:06	0.99999	0.99998	0.99117
	5e-4	50	0:25:29	0.99999	0.99998	0.99559
	1e-3	100	0:23:53	0.99998	0.99985	0.19334
	1e-2	1000	0:02:23	0.35532	0.02900	0.00065
1e-4	1e-4	1	0:15:08	0.99973	0.99992	0.71259
	5e-4	5	0:05:10	0.99999	0.99998	0.99687
	1e-3	10	0:03:57	0.99998	0.99992	0.39737
	2.5e-3	25	0:03:07	0.99992	0.72059	0.02915
	5e-3	50	0:02:35	0.97173	0.47166	0.00438
	1e-2	100	0:00:22	0.35368	0.03011	0.00070
	1e-1	1000	0:00:22	0.02381	0.00002	0.00004