Abstract
Computational model using continuous source functions along the fibre axis is presented for simulation of temperature/heat flux in composites reinforced by short fibres with large aspect ratio. The aspect ratio of short fibres reinforcing composite material is often as large as 103 : 1–106 : 1, or even larger. 1D continuous source functions enable simulating the interaction of each fibre with the matrix and also with other fibres. The developed method of continuous source functions is a boundary meshless method reducing the problem considerably comparing to other methods like FEM, BEM, meshless methods, or fast multipole BEM formulation.
1. Introduction
Short fibres are widely used as reinforcing materials in most advanced composites [1]. In last decades, fibre-reinforced composites have been widely used in engineering applications due to the superiority of their electrothermomechanical (ETM) properties over the single matrix. Because of these properties very large gradients are localized in all ETM fields along the fibres and in the matrix. Particularly, composite materials reinforced by short fibres/tubes (CRSF) are often defined as materials of future. Understanding the physical behaviour of these fibre-reinforced composites is essential for structural design.
Accurate numerical simulation of the fields is important for correct analysis and design of the material behaviour. Among the existing numerical methods, finite element method (FEM) [2, 3], boundary element method (BEM) [4, 5], and meshless methods [6] can be used for simulating the CRSF behaviour. However, these classical numerical methods are suitable for the lowest scale simulation only. All the methods mentioned above require millions or even billions of equations after numerical discretization to obtain a sufficiently accurate solution of this kind of computer simulations.
Recently, the fast multipole method (FMM) [7] was developed to increase the efficiency of numerical models. The fast multipole boundary element method (FMBEM) [8], which was used for simulation of these problems, reduces the computational cost for the far field interaction simulations; however, the classical BEM has to be used for simulating the near field interactions and the supercomputers are necessary to solve the problem with good accuracy. The FMBEM is the only numerical method which enables solving representative volume element (RVE) of the matrix containing two to five thousand short fibres (three to eleven million degrees of freedom (d.o.f.)). The problems have been solved on a supercomputers, or on clusters of PCs [9–11]. If a matrix reinforced by carbon nanotubes (CNTs) having 1% volume of CNTs 1 mm long and 20 nm diameter will contain 3.2 × 107 tubes and the CNTs can have general form (curvature and orientation) and can be randomly distributed, then the control volume element (CVE) for homogenisation would contain several millions of CNTs for correct assessment of homogenized properties of the composite.
More recently, method of continuous source functions (MCSF) developed in [12–14] using 1D continuous source functions distributed along the fibre axis enable simulating interactions of a fibre with matrix and with other fibres very effectively. 1D continuous force or heat sources and force or heat dipoles distributed along fibre axes simulate the interaction among each pair of fibres with the matrix in some CVE for linear elasticity or thermal problems. The inter-domain continuity between fibre and matrix is satisfied at collocation points on the fibre boundaries and so the MCSF is a meshless method, which has some common features with the method of fundamental solutions (MFS) [15, 16].
In the present formulation the problem is defined as follows. The infinite homogeneous material of the matrix is subjected to constant heat flow. A patch of regularly distributed fibres parallel to each other and to the temperature gradient and overlapping in their length in this domain is subjected to simulate the influence the fibres to the heat flow. Heat conductivity of fibres is much larger than that of the matrix. The aim of the simulation is to evaluate increase of conductivity of the matrix by the fibres.
The total solution is split into two parts: the first is homogeneous problem without fibres and the second is solving the interaction of fibres with the matrix. The problem is solved iteratively. In the first step, it is supposed that the fibres are superconductors; that is, the temperature of each fibre is constant. Temperatures in centres of fibres are found from the energy balance and distribution of temperature along fibres is solved in each iteration step for finite conductivity of fibres.
In previous research, the MCSF models were compared with computational models of Trefftz-FEM (T-FEM) and adaptive cross approximation BEM (ACA BEM) with the aim to show possibilities of reduction the problem (see more in [17]). All three methods imply considerable problem reduction and the largest one is achieved by the MCSF.
A brief outline of the paper is as follows. Detail description of the model is described in Section 2. Computational results are presented in Section 3 and conclusions about modelling and behaviour of composite material are given in Section 4.
2. Method of Continuous Source Functions for Simulation of Fibres
The basic assumptions are homogeneous and isotropic matrix materials and fibres and the dimensions of the matrix are infinite. To simulate heat conduction behaviour of CRSF by the MCSF with finite conductivity fibres (several orders larger than matrix), the temperature and heat flow in all fibres should be obtained. Heat flow can be computed by the source functions of unknown intensity, which can be solved from the interdomain continuity (on the fibre-matrix interface).
For the simplicity of modelling all fields are split into two parts, the homogeneous stage corresponding to the homogeneous problem of the matrix without fibres and the local part containing the influence of interactions of fibres with the matrix.
In the present model, the source functions (heat sources and heat dipoles) are 1D continuously distributed along the fibre axis (Figure 1) to simulate the interactions of fibres with the matrix and with other fibres in some CVE. Heat source is a scalar quantity. Heat dipole is a heat source and heat sink acting at the same point by approaching each other in some coordinate direction. Mathematically, it is a derivative of the heat source in corresponding direction (i.e., the heat dipole is a vector).
Temperature field induced by a unit heat source acting in arbitrary point of infinite domain is the fundamental solution for heat problems and it is given by where is the distance of the field point and source point , where the heat source is acting at.
Temperature field induced by a unit heat dipole in direction is
Heat flow in direction by the unit heat source (1) is given by where is heat conductivity of the matrix. Usually the heat conductivity of fibre is several orders higher than that of the matrix.
Similarly the heat flow by the dipole (2) is obtained from the second derivative of (1): where is Kronecker delta.
Heat flow through a surface with the normal is defined as
Collocation points (Figure 1) are located on the interface between fibre and matrix. Then the intensities of the source functions can be computed by satisfying the interdomain boundary (continuity) conditions in collocation points on the fibre-matrix interface boundaries. In the present method, four collocation points are used on each cross section of fibre for satisfying the boundary (continuity) conditions in perpendicular direction of fibre-matrix interface boundary. It should be stressed that the ends of a fibre should be in the form of half spheres or cylinders. It is important to satisfy the boundary conditions (b.c.) in these parts as well. Without considering these b.c., the source functions located along the fibre axis may lead to incorrect results in evaluation of the heat conduction behaviour.
In order to find the unknown intensities of the source functions, we have to solve 1D quasi-singular integral equation in the form: where is kernel function, which is corresponding source function in our case, is the unknown intensity of the source function, is the function prescribing the b.c., and Γ is the 1D integration curve along fibres’ axes.
The well-known numerical solution of integral equations is by the BEM and by its more effective form, FMBEM [4, 7–11]. However, no one formulation was used for composites reinforced by fibres with large aspect ratio and authors of this paper do not know any publications to problems like those shown below. Also the MFS [15, 16] was used to solve problems expressed by the boundary integral equations (BIE). In all these formulations, however, very large number of equations is necessary to obtain the solution with good accuracy.
The source functions are presented by 1D quadratic nonuniform rational B-splines (NURBS) shape functions defined by discrete points along the fibre axis. Numerical integration along 1D elements and corresponding nonuniform rational B-splines (NURBS)shape functions [18, 19] are the integrands in (6). They are quasi-singular and complicated for analytical evaluation. After numerical experiments on behaviour of source functions in composite material and on their influence to physical fields, the following conclusions are drawn for generation of collocation points and numerical integration: Very fine points defining NURBS for source functions are to be chosen in the end parts of the fibres because of large gradients in all fields in fibres and matrix. As there are very large differences in denominator (distance between source and collocation points) the whole integration path is split into integration elements. In our models the smallest element (closest to the collocation point) has to be equal to the fibre diameter and the others are about as large as the distance of its closest point from the collocation point. In this way same number of Gauss integration points can be used for all elements with about equal numerical error from all integration elements in the model. The five-point Gauss quadrature is enough to give good accuracy in this case.
All convenient distribution of nodal points along integration curve (fibre axes), order of integration (number of Gauss integration points), and collocation points (along fibre-matrix interface) are decisive for numerical stability and accuracy of the model. All these problems were in detail studied and discussed in our previous publications [12–14, 20–22] and presented in several conferences on computational mechanics.
The boundary conditions along fibres are unknown as we know the b.c. in infinity only. So, they are specified as follows. Temperature in the middle of each fibre is equal to zero; fibres are specified as superconductors; that is, the temperature is constant along fibres and opposite gradient to that in homogeneous material is prescribed along fibres. Heat sources along fibres’ axes serve to satisfy the b.c., but they will give some temperature gradient in cross section direction (CSD) of fibres. Heat dipoles in CSD will serve to give constant temperature, thereby prescribing temperature differences in opposite points of corresponding CSD equal to zero. The temperatures in middle of fibres are obtained by satisfying heat balance along each fibre. It is realized by additional right-hand sides (r.h.s.) in the system of resulting equations that temperature in corresponding fibre is prescribed to one and zero to the other fibres. In this way we have r.h.s. for fibres in a patch. Resulting system of equations is solved in the least squares (LS) sense as we have more equations than unknown source function intensities (we have one equation in each collocation point for temperature and half of it for temperature differences in opposite points). The intensities of heat functions and dipoles are defined by the intensities of shape functions in NURBS for all fibres and for all r.h.s.
The energy balance is received by integrating the heat source function along each fibre, which gets system of equations bringing us the temperatures in the fibres’ midpoints.
3. Results
Three examples showing distribution of temperature and heat flow in fibres are presented. The first example is a patch of 14 fibres and the other two are patches with 38 fibres arranged symmetrically and alternately in 3 layers (3 × 3 × 3, resp., 5 × 5 × 3) as shown in Figure 2. The heat conductivity of fibres is 103 times larger than that of the matrix in the first two examples and 105 times larger in the third example (38 fibres, so as in example 2). The aspect ratio of fibres is 50 : 1 (radius of fibres is equal to 1 and length 100). All units are dimensionless as the problem is linear. The gap between fibres in the fibre direction, , is 40, and the distance between fibre layers in and direction is 10.
The next three examples (number 4 to 6, Figure 3) give results for matrix reinforced by fibres with aspect ratio equal to 500 : 1 (radius of fibres is again equal to 1) and there are 3 × 3 × 3, 3 × 3 × 5, and 5 × 5 × 3 layers of fibres aligned as given in Figure 2 so that we have (13, 23, and 37 fibres in corresponding patch of fibres). Distance between layers is 100 in and directions and the gap between fibres in direction is 40. The heat conductivity of fibres is 5 × 104 times larger than that of the matrix. The configuration is different from that in the first three examples.
Heat flow is in the fibre direction and the gradient of temperature is equal to 1 so that the temperature in the centre of fibres is −70, 0, and 70, respectively, in the lower, middle, and upper -layers of fibres in the first three examples and −1040, −520, 0, 520, and 1040 in the last three examples, which correspond to the coordinate in Tables 1–4.
Figure 4 gives the heat flow in fibres and Figure 5 temperature change to the temperature in the centre of fibre along the fibres.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
The larger examples 2 and 3 gave by the numerical solution a system of 13224 × 3762 linear equations which were solved in least square sense once and only r.h.s. was changed during the iteration process. Recall that there were r.h.s. for the problem with fibres. The largest problem, example 2, was solved using homemade program in MATLAB in notebook within 2600 seconds.
There were six iteration steps used in the first two examples and only two in the last example (very large difference in conductivities of the matrix and fibres). Full lines correspond to the last iteration step. Figures 4(b) and 5(b) (second example) show the graphs for the fist (largest values), the second (lowest values), and the last iteration steps.
Table 1 contains temperature increase of the centres of fibres obtained in the last iteration step in the lowest layer of fibres from the temperature of homogeneous material (−70). The problems are symmetric in and direction and antisymmetric in direction.
Maximal values of heat flow in fibres are given in Table 2.
Tables 3 and 4 contain temperature increase and maximal heat flow in the centres of fibres obtained in the last iteration step in the middle layer of fibres from the temperature of homogeneous material (0).
4. Discussion and Conclusions
Recall that the infinite matrix of homogeneous material with the patch of its reinforced part with fibres is simulated in the present models. The b.c. are those defined in infinity by temperature gradient in fibre direction. Such b.c. are simple and do not require additional collocation points on the control volume boundaries of the patch. So, the problem is similar to problem with a homogeneous material with another material with different conductivity in the patch. We can observe different heat flow and temperature distribution in different fibres and influence of the change of conductivity of fibres on the conductivity of the matrix from the figures and tables, especially the following.(1)The convergence of iteration process to find the temperature in the centre of fibres and temperature distribution along fibres is fast and two iteration steps are enough to obtain the results. It is because the temperature of each fibre is near constant in such case. Of course, the conductivity in fibres’ direction of the composite is largest in such case.(2)There is strong interaction not only of the fibres with matrix, but also of the fibres with closest fibres. One can observe this effect from the distribution of temperature along the fibres.(3)Numerical experiments have shown that, if the aspect ratio is very large or the difference between conductivity of matrix and fibres is smaller, there is a “saturation” of the heat flow and the heat flow in the middle part of fibres is constant which correspond to the case of infinite fibres reinforcing the matrix.
One can find by observing the heat flow in figures and maximal values of the heat flow in fibres that the heat flow in fibres is much higher than that in the homogeneous material (recall that it is equal to 1) and conclude how the fibres contribute to the increase of the conductivity of the composite material in the fibre axis direction.
Note that the heat in fibre divided by the heat conductivity of the fibre gives the heat flow in the matrix along the fibre boundary. The heat flow in the first two and also in the last three examples is much large in the first iteration step (corresponding to superconductive material of the fibre) and decreases in the iteration process when the finite conductivity is taken into account by changing the interdomain b.c. The b.c. in the next iteration step are changed very little in the third example; that is, the fibre behaviour is similar to the case of superconductive fibre. Although the heat conductivity in the last three examples is only a half of that case, but the aspect ratio of fibres is ten times larger, the finite conductivity influences the b.c. much more that in the third example.
Note also, especially by comparing the last three examples, that the boundary conditions, caused by the finite number of fibres in infinite matrix, influence both the heat flow and temperatures in fibres quite considerably. This effect is to be expected also in other finite regions by composite materials.
The present model can be used also for homogenization but only for regularly distributed fibres, because the heat and temperature are not homogeneous in the patch. For general case, it is necessary to define the b.c. along the control volume and the collocation points have to be located so that they will satisfy the b.c. in good accuracy because of large gradients in corresponding fields there.
Recall that the fibres are parallel and straight in the present model and quite general problems with irregularly distributed, curved fibres require using CVE containing large number of fibres to define correctly the homogenized material properties. We have also to realize that the dimension of fibres with large aspect ratio results in different material behaviour in the regions with local shape changes to the behaviour of homogeneous material with similar material properties.
Mechanical properties of composite material reinforced by short fibres are even more complicated to simulate and generally require three times more equations to solve the problem [12] as the basic source functions are forces vectors (in the heat flow the heat source is scalar) and its derivatives are tensor functions. However, the behaviour of computational model is very similar to the presented problems.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The support of DSSI, corp., Bratislava and APVV bilateral project of Slovak and Czech Republic of applied research of 7AMB12SK126 (SK-CZ-0008-11) “Computational Models and Applications of Composite Materials Reinforced by Short Fibres” is gratefully acknowledged.