Numerical Implementation

Analytical solutions to a set of boundary integral equations are rare, even with simple geometries and boundary conditions. To make any reasonable progress, a numerical technique must be used. There are basically four issues that must be discussed in any numerical scheme dealing with integral equations. The first and most basic one is how numerical integration can be effected, together with an effective way of dealing with singular kernels of the type encountered in elastostatics. Numerical integration is usually termed numerical quadrature, meaning mathematical formulae for numerical integration. The second issue is the boundary discretization: when integration over the whole boundary is replaced by a sum of the integrations over the individual patches on the boundary. Each patch would be a finite element, or in our case, a boundary element on the surface. Obviously a high-order integration scheme can be devised for the whole domain, thus eliminating the need for boundary discretization. Such a scheme would be problem dependent and therefore would not be very useful to us. The third issue has to do with the fact that we are constrained by the very nature of the numerical approximation process to search for solutions within a certain subspace of L2, say the space of piecewise constant functions in which the unknowns are considered to be constant over a boundary element. It is the order of this subspace, together with the order and the nature of the interpolation of the geometry, that gives rise to the names of various boundary element schemes. Finally, one is faced with the task of solving a set of linear algebraic equations, which is usually dense (the system matrix is fully populated) and potentially ill-conditioned. A direct solver such as Gauss elimination may be very efficient for small- to medium-sized problems but will become stuck in a large-scale simulation, where the only feasible solution strategy is an iterative method. In fact, iterative solution strategies lead naturally to a parallel algorithm under a suitable parallel computing environment. This chapter will review various issues involved in the practical implementation of the CDL-BIEM on a serial computer and on a distributed computing environment.

Completed Double Layer Boundary Element Method

Despite the linearity of the Navier equations, solutions to complex boundaryvalue problems require substantial computing resources, especially in the so-called exterior problems, where the deformation field in the space between the inclusions to infinity must be calculated. In the traditional spatial methods, such as finite difference, finite element, or finite volume, this space must be discretized, perhaps with the help of "infinite" elements or a truncation scheme at a finite but large distance from the inclusions (Beer and Watson, Zienkiewicz and Morgan). There are two important limitations of spatial methods. The first is the mesh generation problem. To be numerically efficient, we must use unstructured mesh and concentrate our effort on where it is needed. Efficient two-dimensional, unstructured, automatic mesh generation schemes exist-one good example is Jin and Wiberg-but unstructured three-dimensional mesh generation is still an active area of research. The second limitation is much more severe: even a moderately complicated problem requires the use of supercomputers (e.g., Graham et al.). Since we are concerned with the large-scaled simulations of particulate composites, with the aim of furnishing constitutive information for modeling purposes, our system will possibly have tens of thousands of particles, and therefore the spatial methods are out of the question. We have seen how the deformation field can be represented by a boundary integral equations, either by a direct method, which deals directly with primitive variables (displacement and trciction) on the surface of the domain, or by the indirect method, where the unknowns are the fictitious densities on the surface of the domain. When the field point is allowed to reside on the surface of the domain, then a set of boundary integral equations results that relates only to the variables on the boundary (displacement and traction, or fictitious densities), and this is the basis of the boundary element method. The boundary is then discretized, and the integrals are evaluated by suitable quadratures; this then leads to a set of algebraic equations to be solved for the unknown surface variables.

Faxén Relations and Ellipsoidal Inclusions

A number of problems involving Stokes flows past rigid particles can be solved by using the well-known Faxén relations and multipole representation of the velocity field. Despite the analogy between the Stokes and Navier equations, the method does not seem to be widespread in the theory of elasticity. In this chapter, we derive Faxén relations and describe some applications involving spherical and ellipsoidal inclusions in an elastic matrix.

Some Applications of CDL-BIEM

This chapter presents some selected three-dimensional applications of the CDL-BIEM in elasticity and Stokes flows, especially to particulate solids for which the method is devised. It is paramount that any numerical method should be validated against known analytical solutions. The method will therefore be benchmarked against known simple solutions of the type reported in chapters 2 and 5. Some selected nontrivial examples, where no analytical solutions are available, will also be presented. The translating sphere is a simple problem with known analytical solution and smooth bounding surface; it is a popular benchmark problem for boundary element codes. Here a rigid spherical inclusion of radius a, centered at x = 0, is displaced by either (1) a constant vector U or (2) acted on by a force F, and we seek the force in the case of problem (1), or the rigid displacement in problem (2).

Load Transfer Problem and Boundary Collocation

The acoustic and mechanical properties of cemented granular materials such as sedimentary rocks are directly related to the load transfer problem between two granules (Stoll). The theoretical description of granular materials has been based on the Hertzian contact problem between two elastic spherical inclusion in an elastic matrix, or its modifications; a review of the contact problem can be found in Johnson. In essence, the deformation problem resulting from a relative displacement between two nearby spherical elastic inclusions is studied, and the load transfer between the two is used to construct a constitutive theory for the particulate solid. In particular, Dvorkin et al. studied the deformation of an elastic layer between two spherical elastic grains, using a two-dimensional plane strain analysis similar to those of Tu and Gazis and Phan-Thien and Karihalo. They concluded that the elastic properties of the cemented system can depend strongly on the length of the cement layer and the stiffness of the cement. The main problem with the method is the assumption that the contribution to the load transfer between the granules comes from the region near contact. The assumption is well justified in the case where the Poisson’s ratio of the cement layer is 0.5 (incompressible), in which case the problem is equivalent to the corresponding Stokes flow problem where exact and asymptotic solutions are available (see, for example, Kim and Karrila). The Stokes asymptotic solution shows that the leading term in the load transfer is of O(є-1), where є is the dimensionless thickness of the cement layer. In the case where the Poisson’s ratio of the elastic layer is less than 0.5, it is not clear that the load is still strongly singular in є, and therefore a local stress analysis in the region of near contact may not necessarily yield an accurate answer, unless є is extremely small. The load transfer problem is pedagogic in that it allows us to demonstrate an effective technique often used in Stokes flow known as the reflection method, which has its basis in Faxén relations (discussed in the previous chapter).

Fundamental Equations

There is a need for theoretical and computational tools that provide macroscopic relations for a composite continuum, starting from a description of the composite microstructure. The outlook for this viewpoint is particularly bright, given current trends in high-performance parallel supercomputing. This book is a step along those directions, with a special emphasis on a collection of mathematical methods that together build a base for advanced computational models. Consider the important example of the effective bulk properties of fiberreinforced materials consisting of fibers of minute cross section imbedded in a soft elastic epoxy. The physical properties of such materials is determined by the microstructure parameters: volume fraction occupied by the fibers versus continuous matrix; fiber orientations; shape of the fiber cross sections; and the spatial distribution of fibers. Hashin notes that “While for conventional engineering materials, such as metals and plastics, physical properties are almost exclusively determined by experiment, such an approach is impractical for FRM (fiber-reinforced materials) because of their great structural and physical variety,” The analysis of warpage and shrinkage of reinforced thermoset plastic parts provides yet another example of the important role played by computational models. The inevitable deformation of the fabricated part is influenced by the interplay between constituent material properties, the composite microstructure and macroscopic shape of the component. Computational models play an important role in controlling these deformations to minimize undesired directions that lead to warpage and shrinkage. The strength, stiffness, and low weight of these materials all result from the combination of a dispersed inclusion of very high modulus imbedded in a relatively soft and workable elastic matrix. It thus appears reasonable, as a first approximation, to consider a theory for the distribution of rigid (infinite modulus) inclusions in an elastic matrix, reserving the bulk of our efforts for the study of the role of inclusion microstructure. A framework for computational modeling has been established for materials processing, using models of microstructure with simplified rules for the motion of the inclusions.

Multipole Expansion and Rigid Inclusions

In the previous chapter, we showed how various integral representations arise in the theory of elasticity. All of these representations can be thought of as a surface distribution of various types of singularity solutions. One can go a step further and consider either a finite pointwise or a continuous line, surface, or volume distribution of some suitable singularity solutions; the points on which the singularity solutions reside may not necessarily coincide with points on the surface of the domain. In fact, one may wish to include in the distribution any solution of the Navier equations, not necessarily singular. Such an approach is known as the method of fundamental solution, or simply the singularity method. It has been well developed in Stokes flow, mostly in relation to slender body motion, but to a lesser extent in the theory of elasticity. The method lacks a rigorous theoretical foundation, but is easy to implement numerically, since there is no singular integral to be considered. In the boundary element literature, it is sometimes known as Trefftz method, or the indirect discrete method (Patterson et al.). In this chapter, we will look at some of the elements of the method, namely, the singularity solutions, the far-field expansion (multipole expansion) and some related topics. In the search for suitable solutions of the Navier equations, we may seek guidance from general representations (i.e., general solutions of the Navier equations); the most well known is the Papkovich-Neuber solution.