PARTICLES INTERACTIONS IN COMPOSITES REINFORCED BY FIBRE AND SPHERICAL INCLUSIONS PARTICLES INTERACTIONS IN COMPOSITES REINFORCED BY FIBRE AND SPHERICAL INCLUSIONS

we a new Method of Continuous Source Functions (MCSF) to modelling of such like composites by ﬁnite length ﬁbres with a large aspect ratio and composites reinforced by spherical inclusions. The source (forces dipoles) are continuously distributed along the ﬁbre axis (i. e. outside of the domain, which is the domain of the matrix) and their intensities are modelled by 1D quadratic elements along the axis in order to satisfy continuity conditions between the matrix and ﬁbre. The spherical inclusions are modelled by a triple dipole located in the centre of the particle and the intensities of the dipole can be computed using a small number of collocation points on the particle


Introduction
Composite materials reinforced by stiff particles or fibres are important materials possessing excellent mechanical and also thermal and electro-magnetic properties. Such composites contain huge number of reinforcing elements with large gradients in all fields in small parts of the matrix (in micro scale) around the reinforcing elements and accurate computational models are important for homogenization of material properties in macro scale (adjustment of local stiffness of such material) and for evaluation of material strength.
Mechanical behaviour of composites under loading is extremely complex and can only be understood if the observed behaviour is interpreted in terms of micro-mechanical or macro-mechanical analyses. The fibres can be in form of cylinders or with hemispheres in the tips. The tips cause large gradient in all displacement, strain and stress fields not only in a close distance to the fibre, but also in a relatively far distance perpendicular to the axis of the fibre. It is very important to accurately satisfy the continuity conditions in the tips of the fibre. Correct simulation of these fields is important for simulation of interaction of fibres and for evaluation of stiffening effect. Mechanical properties and possible failure modes of these composites can be predicted early during the design stage using modelling techniques [1].
A suitable method for solution of the above problems is a multiregion approach that often leads to inaccurate results, particularly, when there is a large difference between the material properties of the matrix and that of the fibre resulting in coefficients in the system matrix differing by orders of magnitude. However, with a relatively large number of fibres in a given problem, this type of approach is not very feasible because a very large amount of computing resources as well as substantial modelling efforts are necessary. It is well known that using a volume element approximation such as FEM, hundreds or thousands of elements are necessary to achieve a required accuracy even for a simple problem.
The classical Eshelby solution [2] was obtained for an elastic isotropic inclusion in an infinite elastic matrix. The treatment of the RVE as an infinite space implies that the inclusion concentration is diluted and, therefore, a direct application of these results to the case of finite inclusion concentration is only approximate. An improved model was suggested by Mori and Tanaka [3]. Their method also assumes the absence of all inhomogenities, but it includes a certain effect of inhomogeneity by taking an average strain in the matrix phase when all the inhomogenities are present. Recently, Sauer [4] solved the elastic field of an idealized, spherical, finite RVE embedded in an infinite, homogeneous, isotropic medium using Boundary Integral Equations (BIE). A solution is found which satisfies the continuity of displacements and traction fields across the RVE/composite interface. However, the model is simplified and does not take into account the interaction of discretely distributed particles in the matrix and calculates the Eshelby tensor from simplified Dirichlet and Neumann boundary conditions.
In our presentation we will use the method of continuous source functions (MCSF) and will show how the Trefftz Radial Basis Functions (TRBF), i. e. RBF satisfying the governing equations (which can be the fundamental solutions, or more general functions, dipoles, dislocations, etc.) can be used to increase the efficiency of simulations. The most efficient methods will be those which will best approximate both domain variables and boundary conditions. The TRBF are source functions having their source points outside the domain. Special attention will be given to the application of the TRBF in the form of dipoles to the simulation of composites reinforced by particles and/or short fibres.
Compared to the Method of Fundamental Solutions (MFS) [5,6] which does not require any integration, the MCSF requires an integration along 1D and 2D element. The integrals are quasisingular and quasi-hyper-singular and the numerical integration is computationally cumbersome and inefficient. The analytic integration using a symbolic manipulation is a very efficient tool used in the models. A relatively small number (usually fewer than 10) of elements in a fibre is necessary to obtain a good accuracy also by a large aspect ratio (e.g. 1:100).
The proposed method is not fully meshless for these particular models as it requires 1D elements outside the 3D matrix domain, however, the model presents a significant reduction (even by several orders) of the resulting system of equations comparing to FEM, BEM, or other known meshless methods and can be qualified as a Mesh Reducing Method (MRM).

Solution method
The boundary integral equation for analysis of elastic domain containing rigid inclusions is used (1) where K is the kernel function for force intensity or dipole intensity along the fibre, g is a boundary condition on the fibre surface. A lower index s denotes the source point where force is acting and f is the field point where the displacement is introduced. Kernel function is represented by kernel functions for displacement and traction components in the fundamental solution (Kelvin´s solution) respectively, which can be found in Appendix A. For fibre reinforced composites the source functions (forces and dipoles) are continuously distributed along the fibre axis (i. e. outside of the domain which is the domain of the matrix) and their intensities are modelled by 1D quadratic elements along the axis in order to satisfy continuity conditions between the matrix and fibre. A 2D distribution of source functions is selected in the parts of the fibres where large gradients appear.
The RBF can be used also for simulation of the interaction between the matrix and particles with a very large aspect ratio such as composites reinforced with short fibres where the aspect ratio can be 1000:1 or even larger. The inter-domain boundary conditions can be simulated by 1D distribution of the TRBF's (source functions) along the fibre axis. When the TRBF's are approximated by polynomials then the problems lead to evaluation of the following integrals x is the coordinate along the fibre axis, the subscripts s and f denote the source and field point and exponents are integer numbers and y is the distance of the field point from the source point. For computational purpose the integral (2) is transformed to The numerical integration of such integrals would be computationally very laborious because of the quasi-singularities and quasi-hyper-singularities in the integrals; however, analytic evaluation of the integrals containing the kernel function and polynomial approximation of the unknown function is a very elegant way of numerical evaluation of the integrals, if the axis of the fibre is straight, i.e. the value of y is constant in the integrals above.
If the ideas of the TRBF and MFS are used then a simple and efficient formulation can be developed. In the case of composite material reinforced by spheres or particles with an aspect ratio not very different from 1, a particle can be modelled by a triple dipole located in the centre of the particle and the intensities of the dipole can be computed by using a small number of collocation points on the particle boundary (Fig.1). The method of discrete dipoles is very simple and details can be found in [7].
Because of the large aspect ratio, continuity of strains between a matrix and a fibre can be simulated by continuously distributed  [8,9] as they are known from the potential theory) along the fibre axis (Fig. 2). The continuous source functions enable to simulate the continuity conditions with much reduced collocation points along the fibre boundary.

Some Applications and Results
Example 1: First example simulates the interaction of fibres with the matrix and also the interaction of fibres: 1) a patch of non-overlaying rows of fibres as shown in Fig. 3 on the left and 2) a patch of overlaying rows of fibres according to Fig. 3 on the right. In the examples the modulus of elasticity of the matrix was E ϭ 1000 and Poisson ratio ν ϭ 0.3. The matrix was reinforced by a patch of straight rigid cylindrical fibres. The length of fibres was L ϭ 100 and L ϭ 1000 and the radius R ϭ 1. The distance between fibres was ∆1 ϭ ∆2 ϭ ∆3 ϭ 16 and for longer fibres also ∆3 ϭ 200 in the fibre direction. The fibres in the patch contain approximately 1 % of the volume of the composite material. The domain is supposed to be loaded by far field stress σ 33∞ ϭ 10 in the direction (x 3 ), which is also parallel to fibres' axes. The model of the fibre used in these examples contained fewer than 100 unknown parameters (intensities of the source functions) and about 200 collocation points. The problem is solved by the least square (LS) method. The variation of forces for longer fibres is given in Fig. 4.
As the fibres are long and thin, they are much stiffer in the axial direction than in bending and the satisfaction of continuity of displacements, strains and tractions on the surface between the matrix and fibres and corresponding displacements and strains along the fibre would require a very large number of TRBF (source points) to simulate the interaction.
Moreover, in the end parts of a fibre the fields have very large gradients, which increases the difficulties with accuracy and numerical stability of the solution. In our models a continuous distribution of the source points is used for simulation of the interaction. It is possible to use both distributed forces and distributed dipoles along the fibre axis (1D distribution) and oriented in the axis direction in the model. Their role is mainly to satisfy continuity in the fibre axis direction. Continuity in directions perpendicular to the fibre axis is served mainly by the continuous dipoles along the fibre axis, but directed perpendicularly to the fibre axis. Recall that continuously distributed dipoles are derivatives of continuously distributed forces. The distribution is approximated by piecewise quadratic functions with C 0 continuity between the elements. More about the model can be found in [10].
The forces in the fibres are much greater if the fibres overlap then in the fibres without the overlap (Fig. 4) and thus the stiffening effect is considerably influenced by the overlapping. The forces can lead to axial stresses which can exceed the stresses in the matrix by several orders and can cause fracture of the fibres in tension or loss of stability in compression.
Extreme shear forces between the fibre and the matrix can lead to de-bonding of the fibre or to de-cohesion and re-cohesion at the ends and also in the middle of a fibre close to another fibre in materials reinforced with nanotubes, which are typical and very efficient novel reinforcing materials. Note that the large gradients in shear stresses arise not only at the ends of fibres but also in the parts perpendicular to the ends of neighbour fibres and the models are very sensitive to 1D distribution functions along the fibre axis and can cause numerical errors [11].

Example 2:
A Rigid Sphere in Elastic Medium In this example composite material is reinforced by spheres or particles with an aspect ratio not very different from 1. The particle is modelled by a triple dipole located in the centre of the particle and the intensities of the dipole can be computed by using a small number of collocation points on the particle boundary. The interaction of 2 particles modelled by two triple dipoles in their centres and using only 6 collocation points on the particle boundaries in an eigenstrain field is shown in the next figures. Deformation, radial and tangential components of tractions are given as computed from the simplified models for radius of smaller particle equal to one fifth of the radius of a larger one and with the distance between particles equal to one fifth of the radius of the smaller particle (Figs. 5a-d). The undeformed form of particles is red and the green and blue are corresponding fields for a corresponding particle. Intensity of tractions is given by radial distance of the corresponding circle from the undeformed form. Recall

Conclusion
TRBF are shown to be very efficient interpolation functions, which satisfy governing equations inside the domain but they can also satisfy boundary conditions in some part of the domain boundary. The TRBF can be introduced by the fundamental solution (unit force acting on infinite continuum), its derivatives (dipoles, couples, dislocations) in mechanics, thermal, or other source functions in other field problems. The TRBF correctly simulate the decay of field variables and so they can efficiently model any concentrators in field variables. They can be also source functions acting in other domains (Boussinesq-Cerutti solution for half space which can be used for effective modeling of effect of local loading [13,14], analytic solution for layered structures, etc.).
It is demonstrated that the TRBF can be used in connection with boundary collocation methods to simulate a microstructure reinforced by particles using ideas similar to MFS with only single triple dipoles located into centres of the particles to simulate the interaction of particles with matrix and with other particles, as well. No meshing and no integration are necessary. The ideas of Fast Multipole Method (FMM) are also possible to formulate using mechanical principles instead of Taylor series expansions by the formulations. The far field interaction is then introduced by resulting dipole taking into account the force and moment equilibrium.
For simulation of a microstructure reinforced with short fibres 1D continuous distribution (it is the TRBF, too) of source functions was developed by the authors. It can reduce the model comparing to other numerical models by many orders. The forces or dipoles can be used for simulation of interdomain continuity in a fibre axis direction and continuous dipoles in perpendicular directions. The models can be further augmented to simulate composites reinforced by imperfect or curved fibres by using a continuous distribution of couples along the fibre axis in order to keep the moment equilibrium of the fibre reinforcing effect. In this way the fibres with a large aspect ratio like carbon nanotubes (CNT) which are very stiff in the fibre axis direction, but much more flexible in bending can be correctly simulated for interaction with the matrix and with other fibres, too. The examples show how important the correct simulation of all interactions is for a global behaviour assessment.
Numerical models can take into account different topologies (size and distribution of particles) of composite, different materials of each particle and can be a part of multiscale computational models. The reinforcing particles can be on the surface only and they can form surface layers. From the computational point of view the models can define the Functionally Graded Material (FGM) from microstructural changes of material properties in the surface layer.
The paper presents the MCSF for composite materials reinforced by stiff fibres. For modelling the interactions by MCFS we used functions dipoles as a source. The boundary conditions on the stiff parts can be defined in the form of rigid body displacements and by strains. Different source functions define different relations between the components of deformation and stress and thus the satisfaction of all boundary conditions is decisive. In some situations different source functions can contribute to a better numerical stability. The source functions are quasi-singular along the fibre boundary. The analytic evaluation of integrals is simple also for a higher order polynomial approximation of the intensities of the source functions.

Appendix A
This appendix provides the details of kernel functions used in the presented MCSF formulation. The field of displacements in an elastic continuum by a unit force acting in the direction of the axis x p is given by Kelvin solution  where δ ij is the Kronecker's delta.
The displacement field of a dipole can be obtained from the displacement field of a force by differentiating it in the direction of the acting force, i.e.   The displacements (A.1) by a force are weak singular, the displacement gradients, strains and stresses are strong singular. The fields defined by a dipole have one order higher singularity (strong singularity in the displacement field and hyper-singularities in the strain and stress fields). The derivatives in the perpendicular direction to the force define the force couple. The displacement field for the couple is