An efficient SPH methodology for modelling mechanical characteristics of particulate composites

Z.J.Zheng .S.Kulsegrm .P.Chen .Y.Q.Chen ,*

aDepartment of Mechanics.College of Mechanical Engineering.Chongqing University of Technology.Chongqing.400054.China

b Institute of Applied and Computational Mechanics.Cardiff School of Engineering.Cardiff University.Queen’s Buildings.The Parade.Cardiff.CF24 3AA,W ales.UK

cDepartment of Mechanics and Engineering Science.College of Engineering.Peking University.Beijing,100871.China

Keywords: Particulate composites SPH Essential boundary condition Multi-point constraints Master-slave method

ABSTRACT Particulate composites are one of the widely used materials in producing numerous state-of-the-art components in biomedical.automobile.aerospace including defence technology.Variety of modelling techniques have been adopted in the past to model mechanical behaviour of particulate composites.Due to their favourable properties,particle-based methods provide a convenient platform to model failure or fracture of these composites.Smooth particle hydrodynamics (SPH) is one of such methods which demonstrate excellent potential for modelling failure or fracture of particulate composites in a Lagrangian setting.One of the major challenges in using SPH method for modelling composite materials depends on accurate and efficient way to treat interface and boundary conditions.In this paper,a masterslave method based multi-freedom constraints is proposed to impose essential boundary conditions and interfacial displacement constraints in modelling mechanical behaviour of composite materials using SPH method.The proposed methodology enforces the above constraints more accurately and requires only smaller condition number for system stiffness matrix than the procedures based on typical penalty function approach.A minimum cut-off value-based error criteria is employed to improve the computational efficiency of the proposed methodology.In addition.the proposed method is further enhanced by adopting a modified numerical interpolation scheme along the boundary to increase the accuracy and computational efficiency.The numerical examples demonstrate that the proposed master-slave approach yields better accuracy in enforcing displacement constraints and requires approximately the same computational time as that of penalty method.

1.Introduction

The application of meshless methods has demonstrated numerous advantages in many engineering problems concerning crack propagation/growth [1-4].erosion in projectile impacting[5,6].debris generation and propagation [7,8].large/finite deformation of materials[9-12]and moving boundaries[13-15]where the traditional finite element or boundary element method may encounter difficulties.Compared to mesh-based methods.meshless methods have less mesh-dependency and may be more suitable for performing self-adaptation analysis.In addition.the meshless methods can facilitate the construction of higher order interpolation functions [16-18]to improve accuracy.

It is well known that in conventional meshless methods.for example.in element free Galerkin (EFG) method [1-4]and SPH method [19-22].the shape functions do not naturally satisfy the Kronecker-δ property.which prevents direct implementations of essential boundary conditions (EBCs).Several numerical methods including singular kernel function method [23-25].boundary transformation method [26,27].Lagrangian multiplier method[28,29].penalty function (PF) method [30,31].alternatives of Nitsche’s method [32,33].master-slave method [34]and modified collocation method [35]have been developed in the past to implement EBCs in meshless methods.The modified collocation method proposed by Zhu [35]is widely used in meshless local petrov-Galerkin methods (MLPG) [36]and extended by other researchers to enforce essential boundary condition in various elasto-static and dynamic problems [37,38].And among these numerical methods the PF method is the most widely used due to its simplicity and efficiency.Hu [39]developed a novel threedimensional meshless Galerkin method for horizontal well reservoir simulation using PF method.Based on EFG method and penalty function procedure.a numerical approach is proposed by He et al.[40]for topology optimization of geometrically nonlinear structures.Zhang analysed heat transfer in orthotropic structure by using EFG and PF method [41].Ren et al.proposed a dual-support SPH which makes the particle collocation more flexible.and employed PF to implement the hourglass force[42].To handle the variable discontinuity in meshless as well as mesh-based methods,Ren et al.proposed a nonlocal operator method and employed PF to implement Dirichlet boundary condition and the operator energy functional [43].However.in PF method.the accuracy of the constraints depends on the penalty coefficient or penalty factor and it is difficult to easily predict the optimum value of penalty factor in advance.Furthermore.the PF method may introduce truncation errors and can increase the condition number of the system stiffness matrix [44].

In this paper.to develop an accurate and efficient modelling platform for simulating mechanical properties of particulate composites,a master-slave approach is proposed for imposing EBC and other linear displacement constraints in the context of SPH method.A Lagrangian kernel with corrected SPH format (CSPH).which provides high accuracy [20,21]and eliminates tensile instability[20,21,31,45].is employed in the present approach.The kernel function and its gradient can be corrected to meet the completeness condition by various procedures[20,21,31,42].Although these corrections ensure that the computational model passes patch test and thus enable the solution to be more stable.the corrected SPH formulation may still suffer from spurious mode due to rank deficiency [20].Many techniques.such as stress points [22,45].least square stabilization [20,21].and introduction of hourglass energy[42]or differential operator energy [43].have been proposed to address these issues.These modifications may lead to some variations in the structure of tangential stiffness matrices.However,the difference in governing equation does not essentially influence the treatment of the proposed method provided that the formulation of the nodal displacement approximation is the same.Therefore.for simplicity.the fundamental form of Lagrangian CSPH [20,21,31]is employed in the present analysis.In the proposed master-slave approach the displacement constraints are treated as linear multi-point constraints and the master-slave (MS) method [46]is utilized to build linear displacement relationships between master and slave degrees of freedom (DOFs).Firstly.the total DOFs are appropriately categorized into master DOFs and slave DOFs.and subsequently the slave DOFs are represented as linear combinations of the master DOFs.The global equilibrium equations are then congruently transformed by employing the linear displacement relations.Finally,the master DOFs are solved from the transformed equilibrium equations and slave DOFs are obtained from the displacement relationships between master and slave DOFs.Numerical tests show that the proposed approach enforces the displacement constraints exactly,and does not significantly change the condition number of the stiffness matrix [34,46,47].The numerical simulation of eigen-strain problem and debonding of particles in particulate composites are investigated to demonstrate the advantages of using the proposed MS approach in the context of SPH method.

This paper is organized as follows.Section 2 presents a brief summary of MS method.In Section 3.the MS method based governing equations for the implementation of EBC and the interface displacement continuity conditions are formulated in terms of a Lagrangian corrected SPH (CSPH) form for modelling the matrix and inclusion of particulate composite.In Section 4.two optimization procedures are proposed to improve the computational efficiency of the developed algorithms.Section 5 summarises the concluding remarks.

2.A brief description of master-slave method

In general.the static or quasi-static equilibrium equations in solid mechanics can be expressed as,

where K.u and f are the stiffness matrix,the displacement vector,and the load vector,respectively.Usually,equilibrium condition(i.e.Eq.(1)) is subjected to linear constraints,

where L is a linear(differential)operator matrix,and p is a possible non-homogeneous term.If the constraints are well defined,Eq.(2)can be solved.for example.by Gaussian elimination procedure according to the following equation,

in which umand usare the master and slave DOF vectors.respectively; u=um∪us.um∩us= ?; q is a constant vector of the constraint relationship.Here the subscripts m and s denote the master and slave DOFs respectively.Apparently.the number of slave DOFs are equal to rank(T)as well as rank(L).Eq.(3)is known as master-slave(MS)relation,and the matrix T is called MS relation or transformation matrix.In general.the partitions of master and slave DOFs and the corresponding relation matrix are not unique.The convenience of implementation.numerical efficiency and stability are the main factors to be considered while determining the partitions of master and slave DOFs.

By partitioning ofK,u and f into corresponding master and slave DOFs,the displacement vector u in Eq.(1)can be written in terms of um.T and q as

Now.substituting Eq.(4) into Eq.(1).and pre-multiplying the resulting equation byleads to,

Expanding Eq.(5) yields,

The above expressions can also be obtained through minimizing the potential energy functional

subjected to Eq.(4).Thus,MS approach transforms the equilibrium Eq.(1) into Eq.(6) from which the master displacements can be directly solved.and the slave displacements can be obtained then from Eq.(3).

From Eq.(7) it can be realised that the computation oftransformations TTKsmand TTKssT is efficient if T is sparse and has only a few rows.In other words.the transformation is efficient if the number of slave DOFs is small and each slave DOF is constrained by a small number of master DOFs.

This condition is usually satisfied in traditional mesh-based methods.however this property may not be satisfied in meshless methods.as demonstrated later in this paper.Besides.in mesh based methods.the congruent transformation represented by Eq.(6) is expected to be performed during the assembly of stiffness matrix.However.in meshless methods the transformation is generally performed after the assembly of the stiffness matrix.This is due to the fact that the stiffness matrix constructed in particle or meshless method includes more particles (or points) and consequently the resulting stiffness matrix is not as sparse as that of stiffness matrix obtained in mesh-based methods.

3.Master-slave method in the context of Lagrangian CSPH

3.1.EBC of the Lagrangian CSPH model

The SPH method is chosen here to demonstrate the effectiveness of MS method in enforcing multi-freedom constraint in the context of a meshless method of similar nature.In SPH interpolation.the interpolation value of a nodal displacement of particle a is expressed as,

where Wb(Xa，h)is the corrected and compactly supported weight function of particle b which lies within the domain of influence of particle a[20,21]; Vbis the statistical volume of particle b.The half radius (or smoothing length) h of the domain of influence is conventionally assumed to be constant unless special treatment is conducted as in Ref.[42].The entry (or element) of the stiffness matrix involving particle a and b based on Lagrangian SPH formulation can then be expressed as [31].

The subscript“Γi”is used here to indicate that the particle i lies on a boundary.As stated earlier,DOFs can be partitioned in various ways to build up the corresponding MS relation matrix T.In this investigation,the DOFs are partitioned into two parts:DOFs on the boundaries are treated as slave DOFs and the rest as master DOFs.Hence Eq.(11) can now be re-written as

where a and b represent particles corresponding to master and slave DOFs,respectively;Eq.(12)can now be written in matrix form as,

in which.the subscript “?！?indicates the boundary nodes and “I”indicates the inner domain nodes.

If ΦΓΓis non-singular.uΓ can be obtained from Eq.(13) as:

Noting that the corrected weight functions have hat-like shape,the diagonal entries of the matrix ΦΓΓare generally the largest in each column or row for the above partition of DOFs.From the numerical investigations it was found that this partition is most likely to produce higher level of accuracy and efficiency in determining the MS relation matrix

3.2.Imposing interface conditions

MS method can also be used for imposing other displacement constraints.for example.the interface displacement continuity conditions.which are often encountered in numerical simulations of mechanical behaviour of composite materials.

The interface has significant effect on the macroscopic performance of particulate composites.The thickness of the interface is often ignored for convenience in the numerical simulations.The interface stress model and the linear spring model are often used to model the stress or displacement jump across the interface.In most of the cases.the displacement continuity conditions across interfaces are considered to be linear[31,48-53]as described by Eq.(2).Therefore,the displacement continuity conditions can be easily implemented with the proposed MS method.

One of the ways in which interfaces can be modelled in meshless method is by representing the interface using particle pairs as illustrated in Fig.1.

The two particles of an interface pair are located opposite to each other as indicated in Fig.1.For simplicity.the interface pair particles are assigned same coordinates.If the interface is intact(or undamaged),the displacements of the two particles of the interface pair are identical.i.e.

Here the superscript“+”denotes the particle on the matrix side of the interface.while “- ” denotes the particle on the inclusionside.The SPH discretization of Eq.(15) can then be written as,

Fig.1.Interface and particle pairs.

where N+and N-are the numbers of particles in the supporting domain of the particle in matrix and inclusion,respectively.Eq.(16)represents the discretized form of the linear multi-freedom constraints on the displacement of the interface.By grouping the particles involved in Eq.(16)into three different categories,i.e.,the matrix interfacial particles (denoted as “B+“).the matrix interior particles(denoted as “I+“),and the particles in inclusion(denoted as “N-“).Eq.(16) can be expanded as,

or in matrix form,

For simplicity.in the present analyses.it is assumed that the interfaces in the particulate composites under consideration are not in the vicinity of the boundary.In that case.the compact characteristics of support domain of the weighting function W(X，h) leads to,

Therefore MS relations for interface displacement continuity condition (Eq.(19)) and those for EBC (Eq.(14)) are not coupled.Hence.they can be combined in the following manner:

Let [uI+uN-]Tbe um.The final equation can be obtained by substituting Eq.(21) into Eq.(7) as,

where,

Fig.2.Inclusion in an infinite plate with constant Eigen-strain.

For simplicity,a failure criterion based on maximum normal(or threshold)stress is employed in the case of particulate composites to simulate the failure of an interface bond.This failure criterion assumes that whenever the normal stress in the interface particle pair is larger than a threshold value,the interface particle pair fails or in other words debond.When an interface pair debonds.the corresponding multi-freedom constraint will be removed subsequently.

For simplicity.Eq.(22)above is solved in a quasi-static manner by using Euler method.

3.3.Numerical examples

Example 1.Consider an infinite plate with a circular inclusion of radius R=20 (see Fig.2).A constant dilatation eigenstrain ε*=0.01is assigned to the inclusion.The elastic moduli and Poisson ratios of the inclusion and the matrix are Ei= 1000,νi= 0.28.Em= 900,νm= 0.33.respectively [29,31].A regular distribution of particles as illustrated in Fig.3 is used in the numerical computation for stable solution and better accuracy.

The analytic solutions of displacement and strain for this example are given by,

where the superscript “i” and “m” denote the inclusion and the matrix.respectively.

From the analytical expression Eq.(24).it can be noted that there is a sudden jump in radial strain across the interface.Generally,it is difficult for a meshless method to accurately capture this phenomenon.Though the distributions of the displacement u(r) obtained by both methods are close to the analytic solution given by Eq.(23)(Fig.4(a)),the strain obtained along the radius by MS method is closer to the analytical predictions (Fig.4(b)).

Fig.3.Local particles distribution for an infinite plate with one inclusion(one quarter).

The accuracy of the implementation of interface condition by both methods were further investigated by measuring the relative error in displacement(between each particle of the interfacial pair),defined as,

As shown in Fig.5,the relative error obtained with MS method is close to the value of double precision,while those by PF method are much higher.Another comparison is made to investigate whether the radial displacements of particles located along X-and Y-axis are identical as predicted by the analytic solutions.This radially symmetric nature of the displacement can be used to evaluate the accuracy of the solution.For this purpose.an expression for the relative error is defined as

The results of the computed error on the symmetricity of the displacement are shown in Fig.6.It can be seen from the figure that the displacement computed in MS method is highly symmetric and the values of the relative errors are accurate up to the level of double precision.The errors computed in PF method display larger dispersion and are 2-4 orders higher than that of MS method.

Example 2.Consider a 100 μm × 100 μm × 0.02 μm plate with several randomly distributed elliptical inclusions (see Fig.7).The Young’s moduli and Poisson’s ratios of the inclusion and the matrix material are 150 MPa,1 MPa and 0.3,0.4,respectively[29,31,50].In this example,4664 particles(with 9328 DOFs)are used for the SPH simulation.The top and bottom edges of the plate are subjected to a symmetric tensile displacement of 0.05 μm.The supporting domain of the kernel or interpolation function is set to be a circle with a radius h = 2a.where a is the statistical radius of a particle obtained from the average volume of particle.After several load steps,most of the interface pairs debonded and the rest were intact.Fig.7 illustrates the distribution of the particles after a certain number of load steps.

The displacement vectors obtained for the particles within the inclusions using both MS and PF methods are shown in Fig.8.

There is a significant difference between the displacement fields observed in the inclusions marked with “A” and “B” (compare Fig.8(a)and(b)).In Fig.8(b),the inclusion A and B are subjected to unexpected rotations during the deformation,which are false mode of deformation produced by PF method.In Fig.8(c) and (d).there are no such visible rotations in the deformation.but the displacement vectors observed are still significantly different from that of Fig.8(a).

It was observed that PF method does not provide adequate constraints for accurate deformation if the number of remaining undamaged interface pairs is very few.For example.at least two undamaged interface pairs are needed to constrain an inclusion to deform in a more natural manner; otherwise the inclusion may demonstrate spurious rigid rotations(Fig.9(a)).However,if the two interface pairs are very close to each other.the small distance is most likely to cause the local system to become unstable in a numerical method.Consequently.the constraints provided by PF method cannot fully prevent the inclusion from rotating (see Fig.9(b) and (c)).which is not the case in MS method.Even with three undamaged interface pair particles,the inclusion A in Fig.8(b)undergoes a significant rotation.It is possible to provide strongerconstraints in PF method by increasing the coefficient or penalty factor.however this will increase the condition number of the system stiffness matrix and lead to accumulation of the numerical error.This explains why in Fig.8 the largest penalty coefficient in tern provides the most exaggerated false rotation.In other words,it is difficult to predict the possibility of rigid rotation in PF method and hence the error in the solution of debonding simulation is not well controlled.

Fig.4.Comparisons of (a) radial displacement (b) radial strain between SPH solutions and analytic solutions.

Fig.5.Relative errors of interfacial pairs.

Fig.6.Relative errors between the radial solutions along two axes.

Fig.7.Particles distribution for Example 2.

From the numerical investigations it has been noted that the MS method provides better results with no rigid rotation.This is due to the fact that MS method can impose displacement constraints more accurately and yield a well-defined stiffness matrix.It is demonstrated from the results obtained that the proposed MS method exhibits better computational stability and accuracy than PF method.

The computational efficiencies of both MS and PF methods are further investigated and the results are listed in Table 1 and Table 2.

Though the MS method demonstrates better stability and accuracy.it lacks in computational efficiency compare to PF method.Therefore,further optimization of the solution procedure is needed to improve the practical application of the MS method.

Fig.8.Displacement vector fields.(a)MS method(b)PF method (penalty coefficient =106)(c)PF method (penalty coefficient =104)(d)PF method (penalty coefficient = 102).

Fig.9.The constraints of an inclusion.

Table 1 The computational time for Example 1 (2.6 GHz CPU,4 GB RAM,RELEASE mode).

Table 2 The computational time for Example 2 (2.6 GHz CPU.4 GB RAM.RELEASE mode).

4.Additional remarks on the use of MS method in SPH modelling

4.1.Investigation of properties of MS transformation matrix T in SPH models

To further investigate the implementation of displacement constraints using MS method in the context of SPH method,Eq.(14)and Eq.(19) are re-written as,

Since the weight functions have compact support,Φssand Φsmare sparse matrices,butis in general a dense matrix.Therefore T may contain many non-zero entries.The matrix T(=Φsm)in Example 1 has more than 400 non-zero entries(in each row)to be taken into account while implementing a single constraint.Fig.10 illustrates the distribution of the magnitude of the coefficient (or penalty factor) for a slave particle located at Xs= (0，20)T.The implementation of a single constraint.for example.us-=us+may involve many particles and most of them are within the supporting domain of kernel function located at Xs.However.some of these particles may be very far from Xs.Hence.the transformation as indicated by Eq.(6)becomes computationally very demanding and as a result.T will not be a sparse matrix.This will introduce many non-zero entries into the stiffness matrix during the assembly and will slow down the numerical computation.

Fig.10.Distribution of master degrees of freedoms of particle Xs(0， 20) and logarithms of their coefficients.

In this example.there are less than 30 particles in the supporting area of the particle located at Xs=(0，20)T(blue solid circle in Fig.10).And only these 30 particles are used to implement the constraint us-=us+when the PF method is employed.

Fig.11.Computational time and relative error versus threshold.

By investigating the magnitude distribution of the nonzero entries in Fig.10,one may note that the most of those entries are very small.For example,roughly more than 3/4 of them are smaller than 10-15and contribute very little to the accuracy of the constraint due to numerical truncations.Hence.it is not necessary to take all the non-zero entries of the relation matrix T into account.

4.2.Determination of optimum T matrix by threshold or cut-off value

A natural way to speed up MS method is by employing a cut-off threshold value in computing the entries of the MS relation matrix T.To investigate the performance of such cut-off threshold method,different cut-off threshold values were tested in Example 1.The relative error introduced from cut-off value is measured by,

As the cut-off threshold value increases,the computational time decreases and the relative error increases.Using a cut-off threshold procedure yields faster MS transformation and decreases the total time taken for the whole solution procedure.This is due to the fact that less nonzero entries are introduced into the assembled stiffness matrix by the cut-off threshold.Moreover,the error caused by the cut-off treatment increases with the cut-off threshold value and can be controlled by choosing an optimum cut-off threshold value without compromising the computational efficiency of the solution procedure.

Further.considering the radial decay of the hat-like shape or weighting function used in SPH method.the farther the master degree of freedom is located from the slave.the smaller is its corresponding contribution to coefficient.Thus,the distance between master and slave DOFs can serve as a parameter to determine whether their contribution to coefficients(or elements)of relation matrix T can be neglected or not.As commonly done in numerical linear algebra.a cut-off symbolic matrix of T*can be formed by employing a criterion based on distance to label the retained coefficients before performing the numerical treatment to introduce cut-off threshold to relation matrix T.

4.3.Singular interpolation on boundary

An alternative method to speed up the computational procedure is by modifying the interpolation technique used in SPH method and applying MS method as it is implemented in mesh-based methods.In Ref.[44].the methods to implement boundary conditions in meshless method are divided into two groups:(1)methods that modify the weak form and (2) methods that modify the interpolation shape functions.The proposed approach can be interpreted as a combination of the two approaches.

Further examination of Eq.(27)and Fig.10 reveals that there are number of overlapping support domains of different slave degrees of freedom, and thus several connected domain distributions are formed.Ineach connected domainis full andthe MSrelation matrix Tisnot assparseasexpected.Thismakes the transformation in Eq.(6)less efficient.If no overlap occurs between the supporting domains of salve DOFs, both Φssand its inverse are block diagonal matrices and T will be sparse and can be determinedefficiently.This can be simply achieved by setting the slave degrees of freedom “invisible” to each other.

Fig.12.Supporting particles of a slave freedom.(a) Original supporting particles set (b) modified supporting particles set.

For example.consider a particle indicated by a grey solid circle as shown in Fig.12(a).In its supporting domain the boundary particles (slave degrees of freedom) are indicated by black solid circles and the interior particles (master degrees of freedom) are represented by unshaded circles.There are two slave particles which are located within the supporting domain of the grey particle.Therefore.the supporting domains of slave particles will overlap in this case.Our proposed method here is to set the two slave particles(i.e.solid black particles)“invisible”to the grey slave particle,thus removing the two particles from the summation over particles in the supporting domain of the grey particle (see Fig.12(b)).This can be achieved by modifying the meshless interpolation[1-3].This procedure would enable the MS relation matrix T to be directly obtained without the inverse operation of Φss.

However.the above treatment might cause the supporting domain singular and can yield additional errors as fewer particles are used to reproduce the deformation field.A possible remedy to avoid singularity is to introduce additional particles on the boundary (the unshaded circle with darker outline in Fig.12(b))near the grey particle to compensate for the loss of supporting particles,especially when the slave particles in supporting area are too few (for example.less than 4 or 5).Similar corrections are generally applied to guarantee the accuracies of interpolation[20,21].By using the above proposed approach.the corrected weighting function satisfies,

It is clear from Eq.(29)that the non-diagonal entries ofbe zero.which will significantly speed up the MS method.Compared to the singular kernel function method(SKF)[23,24],which ensures the Kronecker-δ property at boundary.i.e.,

it is much easier to construct a kernel function that satisfy Eq.(29).

4.4.Numerical examples using the modified MS method

To validate the proposed optimization procedures.the debonding process of inclusions from matrix of a composite material is simulated [29,31].A representative volume with a few randomly distributed elliptical inclusions(same as that in Example 2) is considered.Both the EBCs and the interface conditions are tested by using MS methods with and without optimizations.andstandard PF method for comparison.For convenience.the MS method with cut-off threshold applied is denoted as MS1 and,the method with cut-off threshold applied and with modified boundary interpolation treatment is denoted as MS2.

Fig.13.Particles distribution for a composite plate with 9 inclusions.

Fig.14.Macroscopic stress versus macroscopic strain (Example 3).

Fig.15.The stress fields (YY component) after the (a) 10th.(b) 20th.(c) 30th step (blue circle indicate the undamaged interface pairs).

Example 3.In this example.the plate geometry described in Example 2 is considered again.All interfaces are initially undamaged.Interface pairs debond when the normal stress exceeds 0.04 MPa.The top and bottom edges are subjected to incremental tensile displacements with step size equal to 0.05 μm.The particles distribution is identical to that of Example 2 (see Fig.13).

The relation between macroscopic stress and strain for the first 30 steps is shown in Fig.14.

It can be seen from Fig.14 that the MS method and its variants provide similar results.In contrast.the penalty method yields significantly different stress versus strain behaviour.More details ofthe strain and stress fields at the 10th,20th,30th step are displayed in Fig.15.It can be seen that the strain and stress fields,as well as the debonding characteristics obtained by PF method are much different from those obtained by MS methods.This again shows that the PF method may lead to unfeasible debonding simulation in complex scenarios.

Table 3 Time costs of Example 3 (40 steps,2.6 GHz CPU.4 GB Rom.RELEASE mode).

The results obtained by variants of MS methods are generally similar to each other.They yield almost same stress and strain values at locations where the composite has very large deformations.The sequential order in which the particles debond along the interfaces is also identical in the case all MS variants.However.small differences near the interfaces were observed between MS2 and other MS approaches due to the type of SPH interpolation employed at the interface.As a result,the debonding of interfacial pairs by different methods may not exactly coincide at the same displacement step.It can be noted from the last row of Fig.15 that in (a) and (b).the inclusion denoted by “A” debonds completely after the 30 steps,whereas the same inclusion denoted by “B”in(c)does not debond completely.In fact,MS method with modified boundary interpolation treatments predicts that the same inclusion will completely debond at the 31st step.

The individual computational time taken for various steps in the solution process while using all MS and PF methods is listed in Table 3.

In this example.MS method without optimizations is significantly time consuming.The optimizations make MS method more efficient.in all cases considered.MS2 performs almost as efficient as PF method.

5.Conclusion

In this paper.to develop an accurate and efficient modelling platform for simulating mechanical properties of particulate composites,a master-slave method is adopted within SPH formulation for imposing the essential boundary conditions and other linear displacement constraints.Two optimization methods are also proposed to improve the efficiency of the developed master-slave method.

In the past[46,47],the master-slave method was used in meshbased methods and it offered a number of advantages such as (1)high accuracy.(2) no empirical parameters.(3) providing symmetric stiffness matrix.(4)little influence on the matrix condition number [46,47].These advantages are also reflected in the above examples which were numerically simulated using SPH method.However.it may require more computational time to obtain the master-slave relation(or transformation)matrix and if the relation matrix has many nonzero entries [46].the transformation and the solution of the equilibrium equation may be computationally more demanding.This is due to the fact that the displacement is obtained by interpolating more nodal variables than in the case of meshbased method.and the SPH interpolation functions do not naturally satisfy the Kronecker-δ property.As a result.the MS relation matrix possesses too many non-zero terms with very small absolute values.

To improve the computational efficiency of the proposed MS method for SPH computations.two optimization procedures are proposed.One is to employ a cut-off threshold value in determining the MS relation matrix.The effects of this procedure on the accuracy.computational time and the truncation errors are investigated.Generally.the MS method works more efficiently with higher threshold value.however the truncation error becomes larger with the higher threshold value.It was also found that the magnitudes of the non-zero terms of the MS relation matrix partly depend on the spatial distance between the particles (or nodes).Thus.the symbolic counterpart of MS relation matrix can be preconstructed according to the distance of involved particles before any matrix is formulated.

It is shown that the cut-off optimization with proper threshold values can significantly speed up the MS method without reducing the accuracy.

By making use of the flexibility of meshless interpolation techniques,another optimization procedure proposed here is to remove some supporting particles of the slave particle so that the supporting domains of slave particles do not overlap each other.By this approach.it is convenient to obtain the relation matrix as the inversion matrix in Eq.(27)becomes a diagonal matrix.In this case,the MS relations become local and cover smaller supporting domain of the slave particle.which ensures much smaller number of non-zero terms in the MS relation matrix T.Thus.MS method with this optimization is more efficient and its computational time is similar to that of PF method.This treatment is also applicable to other problems in mesh-less methods.To improve the accuracy of the simulation.additional particles are assigned as compensation particles near the slave particle to ensure that there are sufficient number of particles within the supporting domain of each salve particle after removing the particles responsible for overlapping supporting domains.In essence,when the particles are sparse,the insertion of compensation particles is important to improve accuracy.

The examples presented above demonstrated that the optimized MS methods can provide better accuracy for the displacement constraints and the optimization with modified boundary interpolation technique(i.e.MS2)is computationally as efficient as PF method.Therefore,the proposed MS method with optimizations can be a better choice for implementing EBC and other displacement constraints to overcome the potential problems encountered by PF method.

Declaration of interests

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Acknowledgments

This work is supported by National Key R&D Program of China(No.2018YFC0809700 .No.2017YFC0803300 ) and National Natural Science Foundation of China (No.71673158 .No.11702046).