Two-scale modeling of granular materials: A FEM-FEM approach

Abstract

In the present paper, a homogenization-based two-scale FEM-FEM model is developed to simulate compactions of visco-plastic granular assemblies. The granular structure consisting of two-dimensional grains is modeled by the microscopic finite element method at the small-scale level, and the homogenized viscous assembly is analyzed by the macroscopic finite element method at large-scale level. The link between scales is made using a computational homogenization method. The two-scale FEM-FEM model is developed in which each particle is treated individually with the appropriate constitutive relations obtained from a representative volume element, kinematic conditions, contact constraints, and elimination of overlap satisfied for every particle. The method could be used in a variety of problems that can be represented using granular media.

Keywords

homogenization / two-scale / representative volume element / compaction / granular assembly / finite element method

Introduction

In recent years, there has been a dramatic proliferation of research concerned with the compactions of granular matter [1-5]. And the two main numerical methods used to simulate the compactions of aggregates are the discrete element method [6-8] and the lattice model [9]. The algorithm of the discrete element method is based on the finite difference formulation of the equation of motion and avoids the inversion process for stiffness matrices. Such an algorithm enables to handle quite a number of particles in simulations. The lattice model is a stiffness method for discrete element systems. Initial configuration of the model is similar to the lattice finite-element mesh. Except for the solving scheme, both methods are essentially the same. The discrete element approach has been used to analyze the problem of density compaction in a dry granular system under tapping, the compaction dynamics, the jamming phenomena, or even to investigate the total stress of viscous particles subjected to twisting moments.

Although the discrete element method has been widely used to study the mechanical relationships between the forces involved and resulting material deformations, it is difficult to track all the individual units in the system coupling of rigid particle motion and the deformation of deformable particles. However, the effect of deformation or strain-rate of individual deformable granules is considerable on powder compaction processes of several industries, including the pharmaceutical, ceramic, food, and household product manufacturing. This point has been presented by a wide variety of studies as follows. Green bodies prepared by compaction of alumina granules have been made transparent by an immersion liquid technique, and internal structure has been characterized with an optical microscope to study the effect of forming pressure on the internal structure [10,11]. Well-dispersed ceria-gadolinia oxide powders have been obtained from thoroughly isopropanol-washed coprecipitated oxalates, followed by calcination at 800°C. The characteristics of calcined powders and the microstructure of the green compacts have been found to be of great importance in the sintering behavior [12]. To study the effect of binder distribution in a powder granule on the internal structure of the compacts, the internal structure of green bodies with various PVA contents has been examined [13]. The effect of deformation properties of individual granules on the consolidation process of granule beds has been investigated [14]. The microscopic behavior of viscous materials under compaction has been studied [15], with focus on the evolution of the pore structure with increasing pressure, at different strain rates.

In the point of view of tracking the individual deformable particles, the finite element method is superior to the discrete element method. Wu et al. [16-19] have introduced a model whose initial configuration is similar to the rectangular finite-element mesh and showed its application to the problem of failures in asphalt pavements. In the finite element model, a single element represents an individual granule. To simulate the whole compaction process, an automatic mesh upgrade scheme has been generated, consisting of contact detection algorithms, mesh update algorithms, concave four-node element modification, and locking prevention subroutine. However, the problem of the finite element model is that if a single element is used to represent an individual particle, then the solution is not accurate enough. If we increase the number of elements for each granule, then the computation cost is too high to analyze global behavior of a general granular system. Standard finite element analysis codes running on modern large capacity workstations cannot analyze each structural entity from the global system at the microscopic level.

In the past two decades, two novel numerical methods have been developed to deal with multi-scale phenomena ubiquitous in science and engineering. One is called the bridging scale method [20-23]; the other the homogenization method [24-26]. The anatomy of the first method is that the solution of an abstract Dirichlet problem is decomposed as u=\overline{u}+u^{\prime }, where \overline{u} denotes the solution of the coarse-scale finite element model, and u^{\prime } represents the solution of the fine-scale finite element model. The main reason to add u^{\prime } is that it is convenient to track behavior of microstructures for localization problems. However, for a porous medium, u^{\prime }can be used to generate constitutive equations for the global system. It is called the homogenization technique.

Lately, a new numerical method combining the discrete element method and the finite element method is generated to simulate the compaction of aggregates of rigid particles [27,28]. At the microscopic level, the granular structure consists of 2D round rigid grains, modeled by the discrete element method. At the global scale level, the homogenized standard continuum is analyzed using the finite element method. However, constitutive models for granular viscous materials based on the two-scale FEM-FEM approach remain scarce.

In the present paper, a homogenization-based two-scale FEM-FEM model is developed to simulate compactions of visco-plastic granular assemblies. The granular structure consisting of two-dimensional grains is modeled by the microscopic finite element method at the small-scale level, and the homogenized viscous assembly is analyzed by the macroscopic finite element method at large-scale level. The link between scales is made using a computational homogenization method. The two-scale FEM-FEM model is developed in which each particle is treated individually with the appropriate constitutive relations obtained from a representative volume element, kinematic conditions, contact constraints, and elimination of overlap satisfied for every particle. The method could be used in a variety of problems that can be represented using granular media.

Problem statement

The proposed two-scale finite element method in the paper aims at simulating the compaction of granular assemblies. The significance of the method is the conception of two scales. Figure 1 presents the experimental study for granular polyethylene glycol made by Cuitino et al. and a schematic representation of the compaction process [15]. Using this case as an example, we illustrate how to accomplish finite element simulation at the two scales. As shown in Fig. 1, the specimen of granular aggregates is cylinder. We take the bisection profile of the cylinder specimen to make a plane strain analysis. The detail description about the two-scale FEM-FEM model is presented as follows.

The granular assembly is expanded to two levels, the microscopic level and the macroscopic level. The microscopic cell, named the representative volume element (RVE), consists of irregular solid particles. The distribution of particles in the representative volume element is obtained from an arbitrary digital image of a certain kind of granular assembly. Then we attempt to trace out its geometric figure for the sake of finite element simulation about the representative volume element as shown in Fig. 2. According to the homogenization theory, we assume that the macroscopic level is locally formed by the spatial repetition of very small microstructures, i.e., representative volume elements, as shown in Fig. 3. However, in the finite element simulation for the macroscopic level, the macro structure is considered as continuous. The relationship between the two levels is shown in Fig. 4.

To fulfill the whole compaction simulation for granular assemblies, it is necessary to make a connection between the macroscopic and the microscopic levels. The key point is the generation of an equivalent constitutive model of the macroscopic level from the representative volume element. The whole numerical scheme is called representative volume element method (RVEM). The procedure of RVEM is shown in Fig. 5. At both levels, the finite element method is adopted.

Two-scale formulations

In the numerical simulations shown in this paper, homogenization modeling is coupled with RVEM code to simulate the compaction process of granular viscous solids, in which the microstructure is implemented in a direct way from an image of a real piece of material. In the homogenization theory, it is usually assumed that a composite material is locally formed by the spatial repetition of very small microstructures [20], i.e., RVE. These elements, consisting of aggregates at the microscopic level, are used as cells in the global finite element model. The granular system is simulated simultaneously at two different levels. The macro constitutive relationship is generated in light of microscopic analysis. In the proposed numerical model, the material is considered to be granular at the small-scale level and to be a continuum at the large-scale level. In this section, finite element formulations at the micro and macro levels are investigated.

Microscopic formulations

At the microscopic level, particles of the granular assembly are considered as viscous incompressible material. In the point of view of viscosity, we investigate the constitutive equation, the relation of strain-rate to velocity, and the principle of conservation of linear momentum for incompressible flow. On the basis of the governing equations, we develop the numerical model and the corresponding iteration form for the representative volume element.

Governing equations

The principle of conservation of linear momentum of an incompressible viscous material is written as

(1){\sigma }_{ij,j}+\rho X_i=\mathrm{0}inV,

where {\sigma }_{ij} is the stress rate tensor, \rho is the density of the material, and X_i are the external forces. The velocity of the viscous incompressible material is allowed to vary with time. The incompressibility is achieved when the density changes following a fluid particle are negligible.

The constitutive equation for the incompressible fluid is given as

(2){\sigma }_{ij}=-p{\delta }_{ij}+\mathrm{2}G{ϵ}_{micij},

where {ϵ}_{ij} is strain rate tensor, G is modulus of viscosity, and p is the negative mean stress. The stress-strain rate relation can be written as

(3){ϵ}_{micij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(u_{\mathrm{1}i,j}+u_{\mathrm{1}j,i}),

where u_{\mathrm{1}i} is velocity vector at the microscopic level. The boundary conditions are as follows.

(4){\sigma }_{ij}{n}_j={\overline{t}}_{i}on{S}_{\sigma },

(5)u_{\mathrm{1}i}={\overline{u}}_{\mathrm{1}i}on{S}_u,

where n_{\mathrm{1}j} is a unit normal vector to an element on surface S, {\overline{t}}_i and {\overline{u}}_{\mathrm{1}i} are specified surface tractions and velocities, respectively. The Lagrange multiplier p for the constant condition can be incorporated to include compressible and incompressible fields.

The corresponding weighted-integral formulation is written as

(6)\delta \varphi ={\displaystyle {\int }_V\delta {ϵ}_{micij}\mathrm{2}G{ϵ}_{micij}dV}-{\displaystyle {\int }_V\delta {ϵ}_{micij}p{\delta }_{ij}dV}-{\displaystyle {\int }_V\delta p{ϵ}_{micii}dV}-{\displaystyle {\int }_{S_{\sigma }}{\overline{T}}_i\delta u_{\mathrm{1}i}dS}-{\displaystyle {\int }_V\rho X_i\delta u_{\mathrm{1}i}dV}.

The variation of f with respect to p yields

(7)\delta \varphi =-{\displaystyle {\int }_V\delta p{ϵ}_{micii}dV}.

The substitution of the finite element shape functions and the differentiation of Eqs. (7) and (3) yields the virtual work expression in the matrix form as follows.

(8)\delta {\varphi }_e={\{\delta u_{\mathrm{1}}\}}_e^T{[K]}_e{\{u_{\mathrm{1}}\}}_e-{\{\delta u_{\mathrm{1}}\}}_e^T{[Q]}_e{\{p\}}_e-{\{\delta p\}}_e^T{[Q]}_e{\{u_{\mathrm{1}}\}}_e-{\{\delta u_{\mathrm{1}}\}}_e^T{\{F\}}_e.

Numerical method

The assembly of all the element matrices yields the finite element matrix equation, written as

(9)\left[\begin{array}{ll}[K]& {[Q]}^T\\ [Q]& [\mathrm{0}]\end{array}\right]\left\{\begin{array}{l}\{u_{\mathrm{1}}\}\\ \{p\}\end{array}\right\}=\left\{\begin{array}{l}\{F\}\\ \{\mathrm{0}\}\end{array}\right\}.

Clearly, some diagonal terms in the matrix equation are zeros. The system can be solved iterative as follows:

(10)\left[\begin{array}{ll}[K]& {[Q]}^T\\ [Q]& -{\displaystyle \frac{\mathrm{1}}{\gamma }}[I]\end{array}\right]\left\{\begin{array}{l}{\{u_{\mathrm{1}}\}}^{n+\mathrm{1}}\\ {\{p\}}^{n+\mathrm{1}}\end{array}\right\}=\left\{\begin{array}{l}\{F\}\\ -{\displaystyle \frac{\mathrm{1}}{\gamma }}{\{p\}}^n\end{array}\right\},

where n is the current iterative step, and g is a large positive scalar.

From the second equation,

(11){\{p\}}^{n+\mathrm{1}}={\{p\}}^n-\gamma [Q]{\{u_{\mathrm{1}}\}}^{n+\mathrm{1}},

and inserting this in the first equation we have

(12)([K]-{[Q]}^T\gamma [Q]){\{u_{\mathrm{1}}\}}^{n+\mathrm{1}}=\{F\}-{[Q]}^T{\{p\}}^n,

We define the equivalent viscosity matrix of a compressible flow as follows:

(13){[K]}_{equ}=[K]+{[Q]}^T\gamma [Q].

The link between scales—homogenization

In homogenization theory, the total micro structural velocity is expanded asymptotically:

(14)\{u(x,x\text{'})=\{u_{\mathrm{0}}(x)\}+\eta \{u_{\mathrm{1}}(x,x\text{'})\}+{\eta }^{\mathrm{2}}\{u_{\mathrm{2}}(x,x\text{'})\}+\mathrm{...},

where u is the total velocity, u₀ is the macroscopic level velocity, u_i are microscopic velocity field, x is the macroscopic level coordinate system, x′ is the microscopic level coordinate system. The coordinates x and x′ are related by

(15)x\text{'}={\displaystyle \frac{x}{\eta }},

Derivatives with respect to x and x’ are written using the chain rule as

(16){\displaystyle \frac{d}{dx}}={\displaystyle \frac{\partial }{\partial x}}+{\displaystyle \frac{\mathrm{1}}{\eta }}{\displaystyle \frac{\partial }{\partial x\text{'}}},

The small velocity strain-rate vector is written as

(17)\{ϵ(u)\}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\left({\displaystyle \frac{\partial u_{\mathrm{0}i}}{\partial x_j}}+{\displaystyle \frac{\partial u_{\mathrm{0}j}}{\partial x_i}}\right)+\left({\displaystyle \frac{\partial u_{\mathrm{1}i}}{\partial x{\text{'}}_j}}+{\displaystyle \frac{\partial u_{\mathrm{1}j}}{\partial x{\text{'}}_i}}\right)+\eta \left({\displaystyle \frac{\partial u_{\mathrm{1}i}}{\partial x_j}}+{\displaystyle \frac{\partial u_{\mathrm{1}j}}{\partial x_i}}\right)+\mathrm{...}\right],

Eliminating the higher order terms for the sake of simplification yields

(18)\{{ϵ}_{total}\}=\{{ϵ}_{mac}\}+\{{ϵ}_{mic}\},

where

(19)\{{ϵ}_{mac}\}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_{\mathrm{0}i}}{\partial x_j}}+{\displaystyle \frac{\partial u_{\mathrm{0}j}}{\partial x_i}}\right)and\{{ϵ}_{mic}\}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_{\mathrm{1}i}}{\partial x{\text{'}}_j}}+{\displaystyle \frac{\partial u_{\mathrm{1}j}}{\partial x{\text{'}}_i}}\right),

The relation between the macroscopic and the microscopic strain-rates is given as

(20)\{{ϵ}_{mic}\}=[M]\{{ϵ}_{mac}\},

(21){\displaystyle {\int }_{\Omega }{\displaystyle \frac{\mathrm{1}}{V_{RVE}}}}{{\int }_{V_{RVE}}\{ϵ(\delta u)\}}^T[C][M]\{{ϵ}_{mac}\}d{V}_{RVE}d\Omega ={{\int }_{\Gamma }\{\delta u\}}^T\{t\}d\Gamma ,

which is the variational equilibrium equation at the continuum level. The interior integral is then the macroscopic viscosity:

(22)[\overline{G}]={\displaystyle \frac{\mathrm{1}}{V_{RVE}}}{\displaystyle {\int }_{V_{RVE}}[\overline{C}]}d{V}_{RVE},

where V_RVE is volume of a representative volume element, and [\overline{C}]=[C][M].

The equivalent macroscopic strain-rate tensor E is defined as

(23)E_{ij}(u)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{V}_{RVE}}}{\displaystyle {\int }_{{\Gamma }_e}(u_{\mathrm{1}i}{n}_j+u_{\mathrm{1}j}{n}_i)}d\Gamma ,

which can also be written

(24)E_{ij}(u)={\displaystyle \frac{\mathrm{1}}{V_{RVE}}}\left({\displaystyle {\int }_{V*}{ϵ}_{micij}(u)dV+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle {\int }_{{\Gamma }_i}(u_{\mathrm{1}i}{n}_j+u_{\mathrm{1}j}{n}_i)d\Gamma }}\right),

where V_RVE = volume or area of elementary cell, {\Gamma }_e = external boundary of V_RVE. These two terms correspond to the average of the microscopic strain rate of the matrix and of the hole respectively (V* = V_RVE ‒ V_c, V_c standing for the porous area).

The constitutive equation of the homogeneous equivalent material, relating the macroscopic stress \Sigma and strain-rate E is given as

(25)\{\Sigma \}=[\overline{G}]\{E\},

Generation of the stiffness matrix [\overline{G}] is the main point of this study. The scheme to calculate the stiffness matrix is introduced as follows. First of all, the microscopic simulation is made to get configuration in each time step. Secondly, the simple mechanical tests are made to generate the constitutive model for each time increment based on these configurations. Then, the macroscopic simulation is made using these stiffness matrices.

Macroscopic formulations

At macroscopic level, the granular assembly is regarded as linear viscous compressible continuum. Therefore, the constitutive law is different from the one at microscopic level. The different numerical model is used here. Because the viscous modulus matrix at macroscopic level is unknown, we assume the viscosity matrix is equivalent to the equivalent stiffness matrix obtained from the representative volume element.

Governing equations

The governing equations for a compressible material in a region V bounded by a surface S are

(26){\Sigma }_{ij,j}+\overline{\rho }{\overline{X}}_i=\mathrm{0}\text{in}V,

(27)E_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(U_{i,j}+U_{j,i})\text{in}S,

(28){\Sigma }_{ij}{n}_j={\overline{T}}_i\text{on}{S}_{\Sigma },

(29)U_i={\overline{U}}_i\text{on}{S}_U.

The constitutive equation of the homogeneous equivalent material, relating the macroscopic stress \Sigma and strain-rate E is given as

(30){\Sigma }_{ij}={\overline{G}}_{ijkl}{E}_{kl},

For an equivalent homogeneous matter, both stress and strain tensors are symmetric. So the constitutive law can be rewritten as

(31){\Sigma }_i={\overline{G}}_{ij}{E}_j,

where i, j = 1,…,6 for three-dimensional problems, and i, j = 1,…,3 for two-dimensional cases.

The Euler equation of the functional for the plane stress problem is written as

(32)\Phi ={\displaystyle {\int }_V{\overline{G}}_{ij}{E}_i{E}_{j}dV}-{\displaystyle {\int }_{S_{\Sigma }}{\overline{T}}_i{U}_{i}dS}-{\displaystyle {\int }_V\rho X_i{U}_{i}dV},

when minimized with respect to all compressible flow fields satisfying the prescribed boundary values {\overline{U}}_i on S_U.

(33)\delta \Phi ={\displaystyle {\int }_V\delta E_i\mathrm{2}{\overline{G}}_{ij}{E}_{j}dV}-{\displaystyle {\int }_{S_{\Sigma }}{\overline{T}}_i\delta U_{i}dS}-{\displaystyle {\int }_V\overline{\rho }{X}_i\delta U_{i}dV}=\mathrm{0.}

Numerical method

The substitution of the finite element shape functions and the differentiation of Eqs. (32) and (28) yields the virtual work expression in the matrix form as follows.

(34)\delta {\Phi }_e={\{\delta U\}}_e^T{[\overline{K}]}_e{\{U\}}_e-{\{\delta U\}}_e^T{\{F\}}_e=\mathrm{0.}

The assembly of all the element matrices yields the finite element matrix equation, written as

(35)[\overline{K}]\{U\}=\{F\}.

Numerical experiments

Before achieving the macroscopic compaction simulation for granular assemblies, the compaction to a microstructure unit (RVE) is implemented. Then by utilizing the relative results from microstructure numerical experiments, it is capable of continuing the next compaction experiment at macroscopic level. Therefore, from microscopic compaction experiments to macroscopic compaction experiments, the whole two-scale finite element compaction simulation for granular assemblies is accomplished finally.

Microscopic compaction experiments

First of all, the vertical compaction experiment to RVE is carried out for obtaining the microstructure compaction curves. Next, in order to get the macroscopic material modulus matrices of RVE ([\overline{G}]) at each step of the process of compaction, we take the vertical condensation, lateral condensation and shearing tests on RVE at each time step.

Vertical compaction experiment of RVE

As the microstructure develops throughout the process of compaction, the compaction curves for RVE are obtained. Meanwhile, we can get the number of time steps when the compaction is accomplished as well as the geometrical characters of RVE at each time step.

The compaction of time-dependent viscous granular materials is simulated step by step using the automatic adaptive mesh generation schemes (Fig. 6). Vertical loading conditions and boundary conditions of RVE are shown in Fig. 7, which is the same as those of the macroscopic level.

In the experiment, it is necessary to achieve three items. First the total number of steps when finishing compaction; Second the density ratio (V_S/V) as well as the corresponding value of vertical pressure at each step; Third the variation of distribution of particles and geometrical characters about RVE at each time step.

Lateral condensation, Vertical condensation and shearing tests on RVE

The lateral condensation, vertical condensation and shearing tests are implemented at each step of the procedure of compaction. All of them adopt specific displacement boundary conditions for the sake of obtaining the macroscopic material modulus matrix of RVE ([\overline{G}]). It is helpful to the following finite element numerical simulation at macroscopic level. So, eventually we can get the certain number of [\overline{G}] which reflects the material feature of RVE at every time step.

About the two-dimensional problem, [\overline{G}] has 3 rows and 3 columns. In lateral condensation tests, we assume RVE’s global strain-rate is one in x direction and zero in other directions by which we can get the first column elements ({\overline{G}}_{i\mathrm{1}}) of [\overline{G}]. In the same way, we get the second column elements ({\overline{G}}_{i\mathrm{2}}) by controlling RVE’s global strain-rate is one in y direction and zero in other directions in vertical condensation tests and the third column elements ({\overline{G}}_{i\mathrm{3}}) by keeping RVE’s global strain-rate is one in xy direction and zero in other directions in shearing tests.

Macroscopic compaction experiments

At the macroscopic level, the granular matter is regarded as linear viscous compressible continuum. But its viscous modulus matrix keeps changing throughout the process of compaction. Because of the assumption about that the viscous modulus matrix ([D]) at macroscopic level is equivalent to the global material modulus matrix of RVE ([\overline{G}]), we can learn the variation of material features at macroscopic level at every step during compaction by lateral condensation, vertical condensation and shearing tests proposed above. Here, we believe the compaction to RVE and the compaction of macroscopic level are synchronous, which is confirmed reasonably in the following numerical results. Thus, the compaction simulation for macroscopic level becomes easy. We just need to update the viscous modulus matrix step by step during compaction. The loading condition and boundary condition of the compaction to macroscopic level is shown in Fig. 7.

Numerical results

According to the contents of numerical experiments, the numerical results also can be classified into two parts. One is related to the microscopic level, i.e., microstructure cell (RVE). The other is involved with the macroscopic level.

Microscopic numerical results

Vertical compaction curves of RVE

We adopt the finite element method to simulate the process of compaction of RVE and the automatic adaptive mesh generation schemes to trace the variation of RVE in the respect of density ratio, pressure and geometry form.

Numerical results indicate the representative element has been close to the state of complete solid when the loading step reaches 109 as shown in Fig. 8. Then the density ratio and corresponding vertical pressures are presented and the compaction curves are obtained as shown in Fig. 9.

The tendency of compaction curves we get are in keeping with the analysis and experimental results (Fig. 1) of granular polyethylene glycol compaction tests made by Cuitino et al. [15]. In addition, the three stages in the compaction curve proposed by Cuitino et al. also can be found in our compaction curve for RVE as shown in Fig. 9. It verifies reliability and effectiveness of the finite element code RVEM as well as feasibility of the compaction numerical model for the microscopic cell.

Macroscopic material modulus matrix of RVE ([\overline{G}])

By lateral condensation, vertical condensation and shearing tests on RVE at each time step, we get [\overline{G}] of every time step. So there are totally 109 values for each element of [\overline{G}] throughout compaction. We can output the datum by the specific code and research their variation tendency and changing characteristics during compaction. At the same time we are able to gain the relationship between RVE’s global stress {\Sigma }_{ij} and corresponding macroscopic strain-rate E_ij during compaction. Figure 10 shows the variation of elements of matrix [\overline{G}] during compaction.

During the compaction of RVE, the variation tendency of {\overline{G}}_{\mathrm{11}}, {\overline{G}}_{\mathrm{22}} and {\overline{G}}_{\mathrm{21}} is similar. At the initial stage, they increase slowly, and then go up more and more rapidly. It is indicated that the stiffness of the specimen increasingly enlarges along with the reduce of porosity. In contrast, the value of {\overline{G}}_{\mathrm{33}} is much smaller and its tendency is also different. It demonstrates the change of stiffness of the granular system during the compaction. In addition, {\overline{G}}_{\mathrm{33}} shows the relationship between macroscopic shear strain-rate and macroscopic shear stress of RVE. According to its changing trend, it is found that the shearing strength of RVE increase unstably, and it is getting stable at the later stage of compaction. In contrast to other parameters, values and trends of {\overline{G}}_{\mathrm{13}} and {\overline{G}}_{\mathrm{23}} are tiny.

Macroscopic numerical results

When simulating the compaction at macroscopic level, we need to renew the viscous modulus matrix [\overline{G}] step by step. Finally, we can get the displacement and stress of macroscopic continuum at each time step. Figure 11 shows the effective stress of macroscopic level during compaction.

In the process of compaction simulation at macroscopic level, it is indicated that the gap of internal stresses of granular assemblies is getting less and less, and almost disappears at the109th time step. It is indicated that the compaction has been almost accomplished. Based on the effective stress obtained here, global behavior of the granular assembly can be predicted and investigated.

Conclusions

The present study is a significant research on the development of the two-scale FEM-FEM model to simulate the compaction of the granular assembly. After the numerical method is extended to the numerical example, it is indicated that the solutions obtained from the numerical model are in a good agreement with those of the experimental studies in the literature. Our present findings contribute to the field’s understanding of two-scale finite element modeling for granular matters. This method could be used in a variety of problems that can be represented using granular media, such as asphalt, polymers, aluminum, snow, and food products, but it does have some limitations. First of all, the spatial repetition of the representative volume elements is an ideal situation. Not all kinds of granular specimens are suitable for this method. The real granular medium is usually compounded by the geometry of different aggregates, making the analysis difficult. Secondly, the lack of knowledge about nonlinear regression seriously limits the usefulness of the numerical method on heterogeneous granular material. Particularly challenging is the generation of the macro constitutive relation from the micro level. Also, the 2D configuration from the experimental studies would limit the scope of the problem in two-dimension. However, the present study should have provided a solid foundation for additional researches in the future.

References

[1] Caglioti E, Loreto V, Herrmann H J, Nicodemi M. A tetris-like model for the compaction of dry granular media. Physical Review Letters, 1997, 79(8): 1575-1578

[2] Ribière PPhilippeP, Richard P, Delannay R, Bideau D. Slow compaction of granular system. Journal of Physics Condensed Matter, 2005, 17(24): S2743-S2754

[3] Forgerini F L, Figueiredo W. Random deposition of particles of different sizes. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2009, 79(4): 041602

[4] Cates M E, Wittmer J P, Bouchaud J P, Claudin P. Jamming, force chains, and fragile matter. Physical Review Letters, 1998, 81(9): 1841-1844

[5] Dintwa E, van Zeebroeck M, Tijskens E, Ramon H. Torsional stiffness of viscoelastic spheres in contact. European Physical Journal B, 2004, 39(1): 77-85

[6] Bartels G, Unger T, Kadau D, Wolf D E, Kertesz J. The effect of contact torques on porosity of cohesive powders. Granular Matter, 2005, 7(2-3): 139-143

[7] Hu W, Jin F, Zhang C, Wang J. 3D mode discrete element method with the elastoplastic model. Frontiers of Structural and Civil Engineering, 2012, 6(1): 57-68

[8] Lei Yang L, Jiang Y, Li B, Li S, Yang Gao Y. Application of the expanded distinct element method for the study of crack growth in rock-like materials under uniaxial compression. Frontiers of Structural and Civil Engineering, 2012, 6(2): 121-131

[9] Ludewig F, Vandewalle N, Dorbolo S. Compaction of granular mixtures. Granular Matter, 2006, 8(2): 87-91

[10] Uematsu K, Kim J Y, Miyashita M, Uchida N, Saito K. Direct observation of internal structure in spray-dried alumina granules. Journal of the American Ceramic Society, 1990, 73(8): 2555-2557

[11] Uematsu K, Miyashita M, Kim J Y, Kato Z, Uchida N. Effect of forming pressure on the internal structure of alumina green bodies examined with immersion liquid technique. Journal of the American Ceramic Society, 1991, 74(9): 2170-2174

[12] DuránP, Moure CJurado J R. Sintering and micro-structural development of ceria-gadolinia dispersed powders. Journal of Materials Science, 1994, 29(7): 1940-1948

[13] Zhang Y, Uchida N, Uematsu K. Direct observation of non-uniform distribution of PVA binder in alumina green body. Journal of Materials Science, 1995, 30(5): 1357-1360

[14] Naito M, Nakahira K, Hotta T, Ito A, Yokoyama T, Kamiya H. Microscopic analysis on the consolidation process of granule beds. Powder Technology, 1998, 95(3): 214-219

[15] ,Cuitiño A M, Alvarez M CRoddy M J, Lordi N G. Experimental characterization of the behavior of granular visco-plastic and visco-elastic solids during compaction. Journal of Materials Science, 2001, 36(22): 5487-5495

[16] Wu Y C, Thompson E G, Heyliger P R. The compaction of aggregates of non-spherical linear viscous particles. Computer Methods in Applied Mechanics and Engineering, 2003, 192(44-46): 4929-4946

[17] Wu Y C, Thompson E G, Heyliger P R, Yao Z H. The compaction of blended aggregates of non-spherical linear viscous particles. Computer Methods in Applied Mechanics and Engineering, 2004, 193(36-38): 3871-3890

[18] Wu Y C. Automatic adaptive mesh upgrade schemes of the step-by-step incremental simulation for quasi linear visco-plastic granular materials. Computer Methods in Applied Mechanics and Engineering, 2008, 197(17-18): 1479-1494

[19] Wu Y, Xiao J, Zhu C. The compaction of time-dependent viscous granular materials considering inertial forces. Acta Mechanica Solida Sinica, 2011, 24(6): 495-505

[20] Hughes T J R, Feijoo G R, Mazzei L, Quincy J B. The variational multi-scale method—a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 1998, 166(1-2): 3-24

[21] Wagner G J, Liu W K. Coupling of atomistic and continuum simulations using a bridging scale decomposition. Journal of Computational Physics, 2003, 190(1): 249-274

[22] Park H S, Liu W K. An introduction and tutorial on multiple-scale analysis in solids. Computer Methods in Applied Mechanics and Engineering, 2004, 193(17-20): 1733-1772

[23] Kadowaki H, Liu W K. Bridging multi-scale method for localization problems. Computer Methods in Applied Mechanics and Engineering, 2004, 193(30-32): 3267-3302

[24] Guennouni T, Francois D. Constitutive equations for rigid plastic or viscoplastic materials containing voids. Fatigue and Fracture of Engineering Materials and Structures, 1987, 10(5): 399-418

[25] Gărăjeu M, Suquet P. Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles. Journal of the Mechanics and Physics of Solids, 1997, 45(6): 873-902

[26] Danas K, Idiart M I, Castaneda P P. A homogenization-based constitutive model for isotropic viscoplastic porous media. International Journal of Solids and Structures, 2008, 45(11-12): 3392-3409

[27] Miehe C, Dettmar J, Zah D. Homogenization and two-scale simulations of granular materials for different micro-structural constraints. International Journal for Numerical Methods in Engineering, 2010, 83(8-9): 1206-1236

[28] Nitka M, Combe G, Dascalu C, Desrues J. Two-scale modeling of granular materials: A DEM-FEM approach. Granular Matter, 2011, 13(3): 277-281

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg