Journal of Computational Engineering
Volume 2016 (2016), Article ID 2753187, 9 pages
http://dx.doi.org/10.1155/2016/2753187

Research Article

A Monte-Carlo Algorithm for 3D Fibre Detection from Microcomputer Tomography

Robert Gloeckner,¹ Stefan Kolling,² and Christian Heiliger¹

¹Justus Liebig University, Institute for Theoretical Physics, Heinrich-Buff-Ring 16, 35392 Giessen, Germany
²Technische Hochschule Mittelhessen, Institute of Mechanics and Materials, Wiesenstr. 14, 35390 Giessen, Germany

Received 14 March 2016; Revised 15 July 2016; Accepted 11 August 2016

Academic Editor: Amine Ammar

Copyright © 2016 Robert Gloeckner et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A model-based approach to analyze fibre distributions in polymer composites applicable for high fibre content is suggested. The algorithm is a four-step iterative method using Monte-Carlo techniques in order to increase speed and robustness for fibre detection. Samples with up to 20% volume fraction of glass fibres and different matrix polymers (PP, PBT) have been analyzed regarding distributions of orientation and length and thickness of the fibres.

1. Introduction

Short fibre reinforced polymer composites are appealing to mid-sized industry, as their use promises low cost continuous production of parts with enhanced properties. Local fibre length distribution (FLD) and fibre orientation distribution (FOD) [1] determine mechanical properties [2–9] but are strongly affected by production processes [10]. The most commonly used manufacturing process for components made of short fibre reinforced plastics is injection moulding. Hereby, molten plastic material and glass fibres are filled into a mould cavity under high pressure. Particularly in thin walled structural parts, the fibres are oriented in tree layers: two boundary layers and a core layer whose thickness depends on the viscosity of the polymer matrix. The fibres in the boundary layers are mostly oriented in filling direction while the fibres in the core layer are oriented perpendicular to the filling direction. This influences strongly the anisotropy of the composite. Within the development process of moulded plastic parts, flow simulations are performed, for example, in order to avoid weld lines in highly stressed areas by optimization of the positions of the gates. Moreover, the local FOD can also be approximated within such flow simulations. This local FOD information can then be used in anisotropic finite element analysis for the structural part under mechanical loading. However, flow simulations are based on idealized fibre models and lead to uncertainty about FLD as well as FOD [11–13]. A comparison of the FOD data computed by flow simulation with real FOD test data is indispensable. Once the test data is available the material parameters, for example, the fibre interaction coefficient in [12], can be calibrated and flow simulation becomes much more predictive. Thus, the usage of short fibre reinforced polymers in safety related components demands effective and accurate quality testing.

While other methods either lack a complete characterization [14, 15] or require too much preparation and labour time [16–18], microcomputer tomography as a nondestructive testing method is a suitable alternative for those needs [19–22]. Common fibre contents of commercially available composites vary between 10% and 70% weight. Hence, due to clusters and accumulations, reliable methods will be needed that are able to differentiate and locate single fibres therein.

The system used in this work is a X-ray desktop microcomputer tomography (CT). Measurements were performed using a SkyScan 1172. Its microfocus sealed X-ray tube operated at 20–100 kV/0–250 A below a maximum power of 10 W. Typical volume picture dimensions used in this work are 1024 1024 900 cubic voxels of 1.7 m length. This CT together with the algorithm presented in this paper has recently been used in [23] to compute the elastic properties of a short-fibre plastic based on real measured microstructure data by homogenization. In the present work we set our focus on the numerical methodology of fibre detection and discuss the Monte-Carlo algorithm in detail.

2. Model-Based Algorithm

A common approach of model-based image analysis is to extract image primitives and to establish a correspondence between these and the model primitives. Probabilistic object pattern matching speeds up the scanning process, provides good scaling capabilities if done in parallel, and increases tolerance against noise and distortion. For an overview of image analysis of materials structures, the reader is referred to [24] and the references therein.

The fibre model consists of a chord of cylinders of constant radii and various lengths. The radius of scannable fibres is limited by the minimum cylinder length.

All Monte-Carlo pattern-matching processes run independent of each other. Their results are merged in a control-process which detects doublets and optimizes the vector data. Subsequently, this process modifies the voxel data as well as the parameters of the pattern-matching processes and initiates their restart.

2.1. Pattern Matching

Core algorithm steps are Monte-Carlo scanning for separate fibre parts and cylinder integral extrusion for fibre length detection. Scanning for separate fibre parts (index ) is done by calculating radial density functions in randomly chosen spherical volumes (centres ):

where

and is a chosen tolerance.

To approximate radial density, the spherical volume is divided into spherical shell lists of random point-symmetric coordinate pairs within nonequidistant radius-intervals:

where denotes the fibre radius and is the pattern radius of the final sphere. Using correlating point-symmetric coordinates according to (2) removes voxels from neighbour fibres as well as noise. Figure 1 shows the radial density as a function of the integral sphere radius and the cylinder radius , which represents the fibre. The density, that is, the probability that a random point within the sphere is also located within the fibre, is for . For an increasing sphere radius we have a decreasing probability and thus a decreasing density distribution which finally tends to zero.

Figure 1: Geometry detection of fibres using radial density distribution for different cylinder radii . Note the strongly decreasing hyperbolic function for .

The volume of the spherical shell defined by radii and is given by . And the corresponding volume of the spherical cylinder can be approximated by for . Thus, the probability that a randomly chosen point can be found within the cylinder shell is and thus

which is used as a convergence criterion. Only if this condition is satisfied, a separate fibre has been found and the centre of the scanned sphere volume is used for further fibre detection. Otherwise, a spherical volume around another randomly chosen within the sphere is checked.

After locating separate fibre parts at , the fibre orientation vector is calculated either via diagonalizing the matrix of second-order moments or using Jacobian matrices. The fibre length is determined by approximating cylinder integrals using circular coordinate pairs as defined in (6) and using ordinary parameters as angle and radial () and transverse () distances from the origin ():

where the normalized vectors , , and are perpendicular to each other and form a right hand system; that is, .

For a given starting point and direction, the fibre is made up of small cylinders , until the last partial sum falls below a certain fraction of the first partial sum , which determines the length of the fibre . The partial sum is given by

After solving the cylinder integrals, all necessary values (position, length, and orientation) for further fibre analysis are available. At this point it is instructive to conclude the complete algorithm for the fibre identification:

(1): chose randomly within the CT-image.(2)Generate concentric shells with radii at the centre .(3)Generate point-symmetric random points within the shells.(4)Compute mass and centre of gravity of all points within the fibre.(5)Use as a new starting point for generating spherical shells.(6); go to () until (7)Compute the radial density distribution of the random points.(8)If Equation (4) is fulfilled, a separate fibre is detected. Else and chose a new randomly within the shells and go to ().(9)Compute fibre length and orientation by cylinder extrusion.

2.2. Duplicate Detection and Match Optimization

Parallel independent running pattern-matching processes intentionally scan each fibre multiple times starting at random positions. In compounds with high filler fraction, matching quality strongly depends on the starting point of the scan. The control-process has to detect multiple scans of each fibre as duplicates and to select the best match of them.

The pattern-matching processes deliver lists of matching information: bottom coordinates , directions , length , radius , and included coloured voxels alias matching hits of (7).

The control-process uses the geometric matching information to generate vector representations. Each fibre-match is represented as a pair of nested prisms. The outer prism dimensions are equal to those of the fibre . Scaling down radius (index ) and length (index ) of the inner prism is done independently, keeping inner and outer prism-centre upon each other:

This prism-pairs are embedded into a 3D-space. To detect duplicate scan-results of a fibre, inner prisms are tested for overlapping. In order to prevent overlapping tests, all fibres are managed by a binning system, as seen in Figure 2. The binning system consists of blocks subdividing the test volume. Each binning-cell contains information about all objects touching it. This reduces the number of overlapping tests to the objects touching the same subvolumes.

Figure 2: Pairs of outer and inner prisms representing fibres embedded into a binning system.

In case of duplicates, is used as heuristic to sort out the best match.

2.3. Iterative Strategy

The goal of iterative scanning of the volume is to remove well-found fibres, in order to get residual fibres more separated and therefore easier to detect in subsequent scans speeding up the first step of the Monte-Carlo scanning for separate fibres. The procedure is somehow similar to separating bundles in the game of mikado.

As mentioned earlier, the model-based algorithm iterates through a loop of parallel distributed pattern-matching subprocesses, controlled by a centralized main process merging the match-data and modifying the volume image.

The central process also decides which fibres to delete from the volume image and how to adapt the parameters of the distributed pattern-matching before passing them to the subprocesses via network.

In the first iterations a descending fraction,

of possible voxels of a matched fibre is used.

In the case , radius boundaries and length are considered as well.

Depending on how many fibres do contribute to the recent deletion-criteria, parameters for pattern-matching subprocesses are adapted: decrease the pattern radius in order to find separate fibre parts accumulations, increase in (2) in order to tolerate bigger colour-differences, and increase in (6) to enhance tolerance against blurring or cloudy fibres.

The number of iterations needed to detect a specific fraction of fibres obviously depends on the image-quality, that is, the contrast between matrix and filler as well on relation resolution and fibre radius. Figure 3 shows for a voxel size of 1.7  and a fibre radius of 6  that, in case of poly-propylene (PP) containing up to 8%-volume glass fibres, 4 iterations are needed to detect 90% of the fibre volume, whereas poly-butylene terephthalate (PBT) compounds with 10%-volume of glass fibres need 7 iterations in order to detect at least 90% of the fibres-volume.

Figure 3: Iterative detection, fraction of matched voxels for different matrix-materials, and different filler content.

Figure 4 shows for PP with 4%-volume glass fibres the time fractions for each step of the pattern-matching of all subprocesses. Calculating the cylinder extrusion integrals is the most time consuming step. The time consumed is proportional to the amount of undetected fibres. The time needed for searching the volume image to find separate fibres is about constant in each iteration.

Figure 4: Iterative detection timing for different matching subprocesses: high numerical effort for solving cylinder integrals.

2.4. Algorithm Robustness

Reconstructed image volume quality depends on density contrast of filler to matrix, filler content, and image resolution. Low density contrast results in blurred object boundaries. High filler content yields to seemingly connected objects. Low image resolution leads to extra or missing voxels on object boundaries.

In order to separate seemingly connected objects, each object is scanned multiple times (see Figure 5). This automatically happens due to the fact that all pattern-matching processes are independent and randomly choose their starting points .

Figure 5: Separation of virtually overlapping fibres, Step 1: scanning for fibres (spheres) and direction analysis (blue arrows).

Accuracy of orientation- and length-detection may vary for each starting point (Figure 6) and so the controlling process decides upon heuristics regarding which combination of to keep and which one to drop (Figure 7).

Figure 6: Separation of virtually overlapping fibres, Step 2: calculating cylinder integrals and refining orientation.

Figure 7: Separation of virtually overlapping fibres, Step 3: optimizing results by removing duplicate detection.

In order to reduce the effect of blurred object boundaries, cylinder integral directions are varied and the maximum of is searched. Using point-symmetric coordinates allows us to separate fibres very efficiently.

In case of low image contrast between filler and matrix, can be reduced to tolerate missing or extra voxels within the cylinder-volume by increasing the acceptance of partial integral sum .

3. Results

First, the algorithm was tested using computed volume images in order to estimate error magnitudes and effectiveness under various idealized circumstances. After that, short glass fibre reinforced PP composites of 0.5% vol, 1.0% vol, 2.0% vol, 4.0% vol, and 8.0% vol to 14% vol. Finally, short glass fibre reinforced PBT composites including 20% vol were analyzed.

3.1. Fibre Detection Quality

According to the procedure proposed in [25], the algorithm was tested in several steps. It was applied to virtual generated fibre volume images and low filler content compounds and its results were compared to those of manual fibre tracking using slices of tomographic volume images as well as microscopic pictures of residual fibres after removing the matrix material by combustion.

Plotting virtual fibres into a volume (see Figure 8) was done to produce a perfectly known system. Additionally this allows us to analyze the effects of different picture defects like blurring and noise against detection efficiency separately.

Figure 8: Qualitative testing using virtual fibres plotted into a volume image. Detected fibres are displayed as prisms, surrounding their associated voxels.

Compounds with reduced filler content (Figure 9) were produced to monitor each step of the detection algorithm applied to real fibres. Due to the low filler content, the fibres can be separated by the human eye and this allows manual comparison of detection data to real fibre voxels. In contrast to virtual fibres, real fibre radius may not be constant and their curvatures may not vanish.

Figure 9: Qualitatively testing using low fibre content compounds (poly-propylene with 0.5% (vol.) glass fibres): scanning for fibres and orientation analysis.

As shown in Figure 10, aside from pieces of broken fibres, most fibres are detected correctly.

Figure 10: Qualitatively testing using low fibre content compounds (poly-propylene with 0.5% (vol.) glass fibres): cylinder integral calculation.

In the next step, 2D-slices ( layers) of tomographic images of commercial available fibre reinforced composites were used to manually track the fibres in order to compare relative magnitudes of the main components and of the orientation tensor in (11) between the tracking data and the algorithm results.

In Figure 11 results for are shown for several slices along the -axis. In order to check whether fibre length may be an issue for human tracking, the tracked fibres were sorted by length and groups of upper 10%, 20%, 30%, and 50% and all of them were taken into account to calculate the orientation tensor. Within Figure 11, different sizes of the marks “” have been used to illustrate the different fibre lengths. As can be seen, there is no significant influence on the fibre orientation. In order to assign automatically detected fibre data well defined to -slices, the centre point of each fibre was tested to lie within. To examine scattering of the orientation data, the thickness of the -slices was refined to be 100%, 33%, 20%, and 11% of the gap between the tomographic layer pictures.

Figure 11: Comparison of the main orientation tensor element between manual tracking (tomographic slices) and automatic detection of glass fibres (8%-vol) in a poly-amide compound.

Except some slices of the PBT data (Figure 11), orientation information of tracked fibres and automatically detected fibres are in good agreement.

Manually tracking fibres using microscope images of fibre content after probe combustion provides fibre length statistics to match.

3.2. Spatial Averaging

The fibre vector representations which are generated using geometric matching information can be analyzed to get spatial distribution of fibre information . One way to achieve this is to divide the test volume into regions. A fibre is taken into account for analysis if its centre lies within a specific region. The region can be any combination of polyhedron and/or spheres. In most cases, the volume is divided into equidistant boxes (-slices, see Figure 12) in order to examine intermediate layer orientation of fibres within a part.

Figure 12: Equidistant subdividing a volume into cuboids in order to spatialize fibre representations.

3.3. Orientation Analysis

In order to describe the fibre orientation state a second-order tensor,

is used [26]. Therefore, fibres are characterized as cylinders of equal and constant radius , equal length , starting points , and directions to determine the components of in (10) as

Taking into account variations in radius and length , this yields to a volume weighted form:

Figure 13 shows an example of the main orientation tensor components, which can be used directly in this form to validate numerical results from flow simulation.

Figure 13: Main orientation tensor components.

3.4. Length Analysis

For basic validation of length analysis results, several mixtures of virtual fibres with fixed lengths and correlated random orientation were plotted into a 3D-volume of voxels. Figure 14 shows an example detection of fibre-lengths with 200, 300, and 400 voxels. In general the peak of the detected lengths lies within an error of 10%, depending on . Overcompensated lowering of due to bad picture quality might yield to erroneous fibre elongation.

Figure 14: Example virtual fibre length analysis of a mixture of fibres with original lengths of 200, 300, and 400 voxels.

In order to validate the resultant lengths on real world probes obtained by the algorithm, we used a crucible to remove matrix material by combustion and viewed the residual fibres under an optical microscope (Figure 15) where the fibres were tracked manually. Thereby fibres touching the borders were not taken into account. Figure 16 shows the comparison of length analysis between manual fibre tracking and the computation. As can be seen the length distribution obtained by the proposed algorithm is in a good agreement to the experimental results.

Figure 15: Stitched microscope pictures of residual fibres after combustion of the polymer material.

Figure 16: Comparison of length analysis between manual fibre tracking of microscope pictures and automatic CT analysis.

3.5. Radius Analysis

Validating radius results was analogous to the procedure which was used to validate length detections. First virtual fibres were analyzed. Randomly chosen fibres from a small set of configurations (direction, length, and radius) was plotted at different positions inside a voxel-volume. Using the geometric representatives (nested prisms), it is possible to avoid crossings between the fibres. Inverted nested prisms (soft shell inside hard shell) are used in order to adjust minimal radial and tangential distances. Figure 17 shows detected diameter distributions for virtual fibres of 8, 12, and 16 voxels. In general, diameters of perfectly plotted virtual fibres are overestimated by an error of about one voxel size, whereas real fibre diameters are slightly underestimated about one voxel size, too.

Figure 17: Virtual fibre diameter analysis distribution of a mixture of virtual fibre configurations with original diameters of 8.0, 12.0, and 16.0 voxels. Detected fibre diameters are within an error of one voxel size.

Subsequently, different compounds which were reinforced with the same glass fibre material (original diameter 14 ) were analyzed. Figure 18 shows that the mean diameter generally underestimates the original diameter by half of a voxel size (0.9 ). The mean relative error of the radius detection of fibres approximately corresponds to the picture resolution. For a standard resolution of 1.8  and commercial available standard compounds (diameter 14 ) this results in a relative error of approximate <7%.

Figure 18: Several diameter analyses of fibre reinforced compound using the same glass fibre material. Original value of the fibre material is marked with a dashed line.

4. Reproducibility

To test the reproducibility of the results for a given tomographic volume image, a poly-amide compound with 8% volume fraction of glass fibres was analyzed 32 times and the probe volume was divided into 15 slices. For each of the 15 slices all elements of the orientation tensor were calculated and an average deviation of all elements for all slices of <0.015 was achieved. The deviation of all deviations was <0.0033.

Analyzing different tomographic volume images of the same specimen again using 15 slices lead to an average standard deviation for all orientation tensor elements of less than 0.025 with a standard deviation of all deviations of less than 0.011.

5. Conclusion and Outlook

Iterative model-based algorithms to analyze short fibre glass reinforced polymers yield a reproducible and robust method to characterize fibre morphology. The basic algorithm can be adapted to bad picture quality and high filler content, so that partially crossed fibres can be separated directly.

The model-based approach was successfully applied on computed volume graphics as well as on short glass fibre reinforced composites of poly-propylene, poly-amide, and poly-butyl-terephthalate up to commercial filler contents of 30 weight%.

In most cases 3 or 4 iterations are sufficient to detect at least 90% of the total fibre volume.

The results of the algorithm were validated by comparing them to manual tracked fibres within tomographic slices and microscopic images. Reproducibility was verified for repetitive analysis of tomographic volume images as well as for different tomographic volume images of the same specimen.

Accuracy and speed of the analysis depend on the aspect ratio of the fibres and on image contrast. Because its statistical basis, the algorithm is highly scalable. Due to its file based protocol, it can be used in heterogeneous clusters already.

In future, the geometric analysis will be extended to support semiautomatic damage examination. Therefore, spatial inhomogeneities of the fibre distribution as well as additional information about transversal and longitudinal distances of the fibres in combination with their length distribution will be used to mark regions of possible fibre fraction. These implementations will additionally be ported to GPUs using CUDA-techniques.

Competing Interests

The authors declare that they have no competing interests.

References

G. Fischer, “Quantitative Ermittlung der Orientierung von Kurzglasfasern mit der Bildanalyse,” Kunststoffe, vol. 77, no. 5, pp. 509–512, 1987.

H. L. Cox, “The elasticity and strength of paper and other fibrous materials,” British Journal of Applied Physics, vol. 3, no. 3, pp. 72–79, 1952.

J. C. Halpin and N. J. Pagano, “The laminate approximation for randomly oriented fibrous composites,” Journal of Composite Materials, vol. 3, no. 4, p. 720, 1969.

H. Fukuda and T.-W. Chou, “A probabilistic theory of the strength of short-fibre composites with variable fibre length and orientation,” Journal of Materials Science, vol. 17, no. 4, pp. 1003–1011, 1982.

I. M. Robinson and J. M. Robinson, “The influence of fibre aspect ratio on the deformation of discontinuous fibre-reinforced composites,” Journal of Materials Science, vol. 29, no. 18, pp. 4663–4677, 1994.

S.-Y. Fu and B. Lauke, “The elastic modulus of misaligned short-fiber-reinforced polymers,” Composites Science and Technology, vol. 58, no. 3-4, pp. 389–400, 1998.

B. Mlekusch, “Thermoelastic properties of short-fibre-reinforced thermoplastics,” Composites Science and Technology, vol. 59, no. 6, pp. 911–923, 1999.

S. W. Jung, S. Y. Kim, H. W. Nam, and K. S. Han, “Measurements of fiber orientation and elastic-modulus analysis in short-fiber-reinforced composites,” Composites Science and Technology, vol. 61, no. 1, pp. 107–116, 2001.

H. R. Lusti, P. J. Hine, and A. A. Gusev, “Direct numerical predictions for the elastic and thermoelastic properties of short fibre composites,” Composites Science and Technology, vol. 62, no. 15, pp. 1927–1934, 2002.

R. P. Hegler, G. Mennig, and C. Schmauch, “Phase separation effects in processing of glass-bead-and glass-fiber-filled thermo-plastics by injection molding,” Advances in Polymer Technology, vol. 7, no. 1, pp. 3–20, 1987.

G. B. Jeffery, “The motion of ellipsoidal particles immersed in a viscous fluid,” Proceedings of the Royal Society of London, Series A, vol. 102, no. 715, pp. 161–179, 1922.

F. Folgar and C. L. Tucker III, “Orientation behavior of fibers in concentrated suspensions,” Journal of Reinforced Plastics and Composites, vol. 3, no. 2, pp. 98–119, 1984.

U. Mohr-Matuschek, Auslegung von Kunststoff- und Elastomerformteilen mittels Finite-Elemente-Simulationen [Ph.D. thesis], Technische Hochschule Aachen, Aachen, Germany, 1992.

R. S. Bay and C. L. Tucker, “Fiber orientation in simple injection moldings. Part II: experimental results,” Polymer Composites, vol. 13, no. 4, pp. 332–341, 1992.

K. K. Kratmann, M. P. F. Sutcliffe, L. T. Lilleheden, R. Pyrz, and O. T. Thomsen, “A novel image analysis procedure for measuring fibre misalignment in unidirectional fibre composites,” Composites Science and Technology, vol. 69, no. 2, pp. 228–238, 2009.

L. Zedler, “Bildanalytische Rekonstruktion dreidmensionaler Füllstoffstrukturen,” AIF 10833B/VI Report, German Institute for Polymers, Darmstadt, Germany, 1998.

C. Eberhardt and A. Clarke, “Fibre-orientation measurements in short-glass-fibre composites. Part I: automated, high-angular-resolution measurement by confocal microscopy,” Composites Science and Technology, vol. 61, no. 10, pp. 1389–1400, 2001.

G. Zak, M. Haberer, C. B. Park, and B. Benhabib, “Estimation of average fibre length in short-fibre composites by a two-section method,” Composites Science and Technology, vol. 60, no. 9, pp. 1763–1772, 2000.

J. Goebbels, H. Heidt, A. Kettschau, and P. Reimers, “Computerized tomography of glass-fiber reinforced plastic components,” in Non-Destructive Testing—Proceedings of the 4th European Conference, vol. 3, pp. 2111–2113, NDT International, 1987.

C. N. Eberhardt and A. R. Clarke, “Automated reconstruction of curvilinear fibres from 3D datasets acquired by X-ray microtomography,” Journal of Microscopy, vol. 206, no. 1, pp. 41–53, 2002.

P. J. Schilling, B. P. R. Karedla, A. K. Tatiparthi, M. A. Verges, and P. D. Herrington, “X-ray computed microtomography of internal damage in fiber reinforced polymer matrix composites,” Composites Science and Technology, vol. 65, no. 14, pp. 2071–2078, 2005.

J. S. U. Schell, M. Renggli, G. H. van Lenthe, R. Müller, and P. Ermanni, “Micro-computed tomography determination of glass fibre reinforced polymer meso-structure,” Composites Science and Technology, vol. 66, no. 13, pp. 2016–2022, 2006.

V. Müller, B. Brylka, F. Dillenberger, R. Glöckner, S. Kolling, and T. Böhlke, “Homogenization of elastic properties of short-fiber reinforced composites based on measured microstructure data,” Journal of Composite Materials, vol. 50, no. 3, pp. 297–312, 2016.

J. Ohser and K. Schladitz, 3D Images of Materials Structures: Processing and Analysis, Wiley-VCH, Weinheim, Germany, 2009.

H. Shen, S. Nutt, and D. Hull, “Direct observation and measurement of fiber architecture in short fiber-polymer composite foam through micro-CT imaging,” Composites Science and Technology, vol. 64, no. 13-14, pp. 2113–2120, 2004.

S. G. Advani and C. L. Tucker III, “The use of tensors to describe and predict fiber orientation in short fiber composites,” Journal of Rheology, vol. 31, article 751, 1987.