A CAS Approach to Handle the Anisotropic Hooke’s Law for Cancellous Bone and Wood
The present research entirely relies on the Computer Algebric Systems (CAS) to develop techniques for the data analysis of the sets of elastic constant data measurements. In particular, this study deals with the development of some appropriate programming codes that favor the data analysis of known values of elastic constants for cancellous bone, hardwoods, and softwood species. More precisely, a “Mathematica” code, which has an ability to unfold a fourth-order elasticity tensor is discussed. Also, an effort towards the fabrication of an appropriate “MAPLE” code has been exposed, that can calculate not only the eigenvalues and eigenvectors for cancellous bone, hardwoods, and softwood species, but also computes the nominal average of eigenvectors, average eigenvectors, average eigenvalues, and the average elasticity matrices for these materials. Further, using such a MAPLE code, the histograms corresponding to average elasticity matrices of 15 hardwood species have been plotted and the graphs for I, II, III, IV, V, and VI eigenvalues of each hardwood species against their apparent densities are also drawn.
The study of material symmetry of 3-dimensional space is of great interest due to having crucial theoretical as well as practical significance. This is because a symmetrical space includes crystals and all homogeneous fields without exceptions: electric, magnetic, gravitational, and so forth.
The variation of material properties with respect to direction at a stagnant point in a material is called material symmetry; for instance, if the material properties are same in all directions at some fixed point, they are called isotropic, whereas if the material properties show variation at the same point, they are called anisotropic . Of course, the familiarity with material symmetries is the best way to categorize the materials. However, according to , many materials are anisotropic and inhomogeneous due to the varying composition of their constituents. In such materials, it becomes ticklish to identify the symmetries or more particularly, the elastic symmetries. The variable composition method to identify material’s elastic symmetry becomes complicated and hence to overcome this difficulty, an approach is developed by , called “averaging anisotropic elastic constant data.” In this approach, the identification of elastic symmetries and method of variable composition are analyzed separately. With the aid of this fabulous approach, a sheer volume of research towards material symmetry has been put forward by various researchers, for example, estimation of the effective transversely isotropic elastic constants of a material from known values of the material’s orthotropic elastic constants , anisotropic Hooke’s law for cancellous bone and wood , the validation of generalized Hooke’s law for coronary arteries , constitutive relationships of fabric, density, and elastic properties in cancellous bone architecture , analysis of elastic anisotropy of wood material for engineering applications , a multidimensional anisotropic strength criterion based on Kelvin’s modes , spectral decomposition of a 4th-order covariance tensor: application to diffusion tensor MRI , a normal distribution of tensor valued random variables: applications to diffusion tensor MRI , and so forth.
Concerning computer assisted analysis of bone material, a few of the articles, like direct mechanical assessment of elastic symmetries and properties of trabecular bone architecture , relationships between morphology and bone elastic properties can be accurately quantified using high-resolution computer reconstructions , and intrinsic mechanical properties of trabecular calcaneus determined by finite element models using 3D synchrotron microtomography , have explored the subject material up to some great extent.
Now, as far as the present proposed research concerned, we turn our concentration towards cancellous bone and wood (soft and hard). According to , the type of material symmetry exposed by bone tissue is often ambiguous. So, in order to attain accuracy, we have to measure the elastic constants of bone specimen, which may be represented as orthotropic. However, the orthotropy of the bone specimen may or may not be of higher degree in all the directions with respect to a fixed point in the specimen. Thereby, for all practical purposes the bone specimen can be represented as transversely isotropic or even isotropic.
Now, the question arises is what is the need of studying the mechanical properties of cancellous bone and wood? The rigorous answer is that the knowledge about mechanical properties of cancellous bone is essential for the determination of bone fracture risk in osteoporosis and other pathological conditions involving impaired bone strength . Also, the mechanical properties of cancellous bone are determined by the properties of its bone tissue and its architecture [4, 14–21]. Further, the wood is a cellulosic, semicrystalline, cellular material and the mechanical properties (e.g., elasticity, strength and rheology, etc.) are generally higher along the bole of a tree than across the bole .
Mechanically, clear wood obeys the law of elastic orthotropic material like the bone does. In general the woody substance also exposes the high toughness and stiffness properties and these properties vary according to the type of wood and the direction in which the woody substance is examined, as the woody substance shows a higher degree of anisotropy.
In addition, to deal with anisotropic Hooke’s law, one should be familiar with the algebra of the 4th-order tensor. We would like to emphasize that the 4th-order tensor algebra is not only involved in the study of material symmetry, but its crucial appearance has been set up in the field of diffusion tensor MRI [9, 10], visualization and processing of tensor fields , biomechanics , tissue mechanics , geophysics , and so forth.
Though an adequate amount of research work has been dedicated towards the algebra of 4th-order tensors, then also, in the following section, we shall present a brief digest regarding this issue, as this issue hits the present study up to some great level.
2. The Algebra of Fourth-Order Tensor
The splendid word “Tensor” has been derived from the Latin word “Tensus”, which is the past participle of “tenděre” and stands for “Stretch.” This word was used in anatomy in the early 1704 to represent muscle and stretches. However, in Mathematics, it existed in 1846 when William Rowan Hamilton has explored his quaternion Algebra. Even though the Hamilton sense regarding tensor did not survive. The recent meaning of tensor is due to Voigt, who used the term “Tensortriple” in crystal elasticity around 1899.
Likewise the second-order tensor, a fourth-order tensor (say), is a linear map from to . That is, it assigns to each second-order tensor the second-order tensor ; that is, where is a vector space and is a linear space of dimension and is equivalent to Also, any 4th-order tensor has the representation where and the set of all four vectors, that is, the set of all linear maps from to is a linear space of dimension and is denoted by .
The outer (or open) product of the two 4th-order tensors is given by the composition In abstract index notations, we have Moreover, the transpose of is a 4th-order tensor and is defined by On , we can define the inner product as The induced norm is and the metric is Hence, with respect to the inner product given by (8), the set forms an orthonormal basis of .
2.1. Symmetries of a 4th-Order Tensor
In case of a second-order tensor, the only known symmetry represents an invariance under the mutual rotation of two indices, while for the 4th-order tensor, there are several notions of symmetry that represent invariance under exchanging pair of indices. Thus, for the 4th-order tensor generally we have three types of symmetries, namely, major, minor, and total symmetries, and these notions of symmetries widely play an important role in the theory of elasticity .
2.1.1. Major Symmetry of a Fourth-Rank Tensor
For , we say that is symmetric if , or in component form and we say that is skew-symmetric if or .
Moreover, we have the decomposition Then the set of all symmetric fourth-order tensors is is a subspace of of dimension .
In the theory of elasticity, this symmetry is often evoked as “major symmetry.” Since the major symmetry represents invariance under exchanging the pair and , then if , then and is an orthonormal basis of , such that The real number is called the eigenvalues of associated with the eigentensor .
Such a representation is called the spectral decomposition of . The trace and the determinant of the fourth-order tensor are defined, respectively, as
2.1.2. Minor Symmetry
The second type of symmetry is called “minor symmetry” and is defined by In component form, , .
Now, the invariance under the exchange of first pair of the indices is called the “first minor symmetry” and the invariance under the exchange of second pair of indices is called “second minor symmetry.”
The set of all 4th-order tensors that satisfy the minor symmetry is denoted by
2.1.3. Total Symmetry
A fourth-order tensor is said to have total symmetry if in addition to satisfying the major and minor symmetries, it also satisfies or in abstract index notations, for some permutation of .
The set of the fourth-rank tensors bearing total symmetry is denoted by Thus any fourth-order tensor can be decomposed into its totally symmetric part and its totally antisymmetric part , that is, .
The components of totally symmetric and antisymmetric parts of a fourth-order tensor are described as  Obviously, a fourth-order tensor is totally symmetric if or equivalently , where stands for a fourth-order null tensor. Likewise is antisymmetric if or equivalently .
Now, it is straightforward to show that if is totally symmetric and is a symmetric tensor, then .
Reference  has also shown that there is an isomorphism (i.e., a 1-1 linear map) between totally symmetric fourth-rank tensor and a homogeneous polynomial of degree four. Also, if a totally symmetric fourth-order tensor is traceless, that is, , then it is isomorphic to a harmonic polynomial of degree four. For this reason, a totally symmetric and traceless tensor is called harmonic tensor.
To meet the objectives of proposed research, let us turn our attention to the flattening (sometimes called unfolding) of the fourth-order 3-dimensional tensor with minor symmetries, as it is customary in the 3-dimensional elasticity theory.
In the following section, we discuss the notion of 9-dimensional representation of a 3-dimensional fourth-order tensor and then particularize this representation to a fourth-order elasticity tensor having minor symmetries due to the well-known anisotropic Hooke’s law.
3. A 9-Dimensional Exposition of a 3-Dimensional 4th-Order Tensor and the Anisotropic Hooke’s Law
In accordance with the evolution of tensor theory, the algebra of the 3-dimensional fourth- order tensor is still not fully developed. For instance, the techniques for calculating the latent roots an latent tensors of a tensor like , or the computation of its inverse (which is called a compliance ), are still not widely available. There are also some certain psychological constraints; that is, we are trained to tackle matrices but not trained to deal with the multiway arrays (sometimes called matrix of matrices). Though, recently Tamra [26–29] has developed a tensor toolbox that runs under MATLAB , in my view this tool still needs some enhancements to handle symbolic tensor data. Moreover, Constantinescu and Korsunsky  have developed some crucial packages, namely, “Tensor2Analysis.m,” “IntegrateStrain.m,” “ParametricMesh.m,” and “VectorAnalysis.m,” which are compatible with Mathematica and Mathematica programming; using these packages enables one to compute symbolic computations regarding elasticity tensors. Anyways, we shall come back to these CAS assisted issues in the subsequent sections, where we shall demonstrate the proposed objectives using CAS.
In order to study material symmetries and anisotropic Hooke’s law, it is customary to transform the 3-dimensional 4th-order tensor as a 9 × 9 matrix, that is, a 9-dimensional representation of a 3-dimensional fourth-order tensor. For this purpose, there is a basic notion of representing a fourth-order tensor as a second-order tensor . According to , if , , then we set Thus, is an orthonormal basis of .
Therefore, any fourth-rank tensor can be assumed as an -dimensional fourth rank tensor , or as an -dimensional second-order tensor , where , .
Again, we have from (8) that hence the two-way map is an isometry.
Let us apply this concept to anisotropic Hooke’s law by assuming that the anisotropic Hooke’s law given below is generalized one and also valid in the case where stress and strain tensors are not necessarily symmetric.
The anisotropic Hooke’s law in abstract index notations is often depicted as where are the components of an elastic tensor. Writing the stress, strain, and elastic tensors in usual tensor bases, we have Now, when working with orthonormal basis, one needs to introduce a new basis which should be composed of the three diagonal elements the three symmetric elements and the three asymmetric elements In this new system of bases, the components of stress, strain, and elastic tensors, respectively, are defined as where all the implicit sums concerning the indices range from 1 to 9 and is the compliance tensor.
Thus, for Hooke’s law and for eigenstiffness-eigenstrain equations, one can have the following equivalences: Using elementary algebra, we can have the components of the stress and strain tensors in two bases: where we have used the following notations:
Now in the new basis , the new components of stiffness tensor arewhere In case if we impose symmetries of stress and strain tensors, let us see what happens with generalized Hooke’s law (22).
In generalized Hooke’s law, when the stress-strain symmetries do not affect the stiffness tensor , the number of components of stiffness tensor in 3-dimensional space is equal to . Now, if we impose the symmetries of stress-strain tensors, that is, and , the will be like .
Moreover, imposing symmetrical connection, that is, , we would have only 21 significant components out of 81. Thus, if we flatten the stiffness tensor under the notion of Hooke’s law, we will definitely have a 6 × 6 matrix having only 21 independent elastic coefficients instead of 9 × 9 matrix having 81 elastic coefficients.
Here, we depict a figure (see Figure 1) which delineates the component reduction process for the stiffness tensor.
Now, with 21 significant independent components, the stiffness tensor can be mapped on a symmetric 6 × 6 matrix.
As the elasticity of a material is described by a fourth-order tensor with 21 independent components as shown in (Figure 1) and the mathematical description of elasticity tensor appears in Hooke’s law, now how to map this 3-dimensional fourth order tensor on a 6 × 6 matrix?
We are fortunate to have a long series of research papers concerning this issue. For instance, [2, 3, 32–39] and many more.
The very first approach to map 21 significant components of the elasticity tensor on a symmetric 6 × 6 matrix was introduced by Voigt  and after that various successors of Voigt have used his fabulous notions for flattening of fourth-order tensors. But Lord Kelvin  found the Voigt notations are inadequate from the perspective of tensorial nature and then introduced his own advanced notations, now known as Kelvin’s mapping. More recently,  have introduced some sophisticated methodology to unfold a fourth-rank tensor called “MCN” (Mehrabadi and Cowin’s notations).
Let us briefly go through these three notions of tensor flattening one by one.
3.1. The Voigt Six-Dimensional Notations for Unfolding an Elasticity Tensor
It is well known that the Voigt mapping preserves the elastic energy density of the material and elastic stiffness and is given by The Voigt mapping receives this relation by the mapping rules Using this rule, we have , , and .
This Voigt mapping can be visualized as shown in Table 1.
Thus, in accordance with Voigt’s mappings, Hooke’s law (22) can be represented in matrix form as follows: where the simple index conversion rule of Voigt (see Table 2) is applied.
But in the Voigt notations, many disadvantages were noticed. For instance,(1)the and are treated differently;(2)the norms of , , and are not preserved;(3)the entries in all the three Voigt arrays (see (35)) are not the tensor or the vector components and thus one may not be able to have advantages of tensor algebra like tensor law of transformation, rotation of coordinate system, and so forth.
3.2. Kelvin’s Mapping Rules
Likewise Voigt’s mapping rules, Kelvin’s mapping rules also preserve the elastic energy density of material under the following methodology: In this way, Naturally, Kelvin’s mapping preserves the norms of the three tensors and the stress and strain are treated identically. Also, the mappings of stress, strain, and elasticity tensors have all the properties of tensor of first- and second-rank tensors, respectively, in 6-dimensional space.
But the only disadvantage of this mapping is that the values of stiffness components are changed. However, there is a simple tool for conversion between Voigt’s and Kelvin’s notations. For this purpose, one just needed a single array (see Table 3). Thus,
3.3. Mehrabadi-Cowin Notations (MCN)
To preserve the tensor properties of Hooke’s law during the unfolding of fourth-order elasticity tensor,  have introduced a new notion which is nothing but a conversion of Voigt’s notations into Kelvin ones. This conversion is done using the following relation: Exploiting this new notion, (35) becomes Further, to involve the features of tensor in the matrix representation (40),  have evoked that (40) can be rewritten as where the new six-dimensional stress and strain vectors denoted by and , respectively, are multiplied by the factors . Also, is a new six-by-six matrix . The matrix form of (40) is given as The representation (42) is further produced in MCN pattern like below: To deal with all the above representations of a fourth-order elasticity tensor as a six-by-six matrix, a simple conversion table has been proposed by , which is depicted as in Table 4.
Even though the foregoing detailed description is interesting and important from the viewpoint of those who are beginners in the field of mathematical elasticity, the modern scenario of the literature of mathematical elasticity now involves the assistance of various computer algebraic systems (CAS), like MATLAB , Maple , Mathematica , wxMaxima , and some peculiar packages like Tensor Toolbox [26–29], GRTensor , Ricci , and so forth, to handle particular problems of the scientific literature.
However, in the following section, we shall delineate a CAS approach for Section 3.
4. Unfolding the Anisotropic Hooke’s Law with the Aid of “Mathematica”
“Mathematica” is a general computing environment intimated with organizing algorithmic, visualization, and user-friendly interface capabilities. Moreover, many mathematical algorithms encapsulated by “Mathematica” make the computation easy and fast [42, 46].
A fabulous book entitled “Elasticity with Mathematica” has been produced by . This book is equipped with some great Mathematica packages, namely, VectorAnalysis.m, Displacement.m, IntegrateStrain.m, ParametricMesh.m, and Tensor2Analysis.m.
Overall the effort of  is greatly appreciated for producing a masterpiece with Mathematica packages, notebooks, and worked examples. Moreover, all the aforementioned packages, notebooks, and worked examples are freely available to readers and can be downloaded from the publisher’s website http://www.cambridge.org/9780521842013.
We consider the following code that easily meets the purpose discussed in Section 3.
Here, kindly note that only the input Mathematica codes are being exposed. For considering the entire processing, an Electronic (see Appendix A in Supplementary Material available online at http://dx.doi.org/10.1155/2014/487314) is being proposed to readers.
First of all install the “Tensor2Analysis” package in Mathematica using the following command
In 1:= << Tensor2Analysis`
We have loaded the above package like Loading the package “Tensor2Analysis.m” by setting the following path:
Next is concerning the creation of tensor. The code in Algorithm 1 generates the stress and strain tensors using the index rules specified by .
Use the MakeTensor command to create tensorial expressions having components with respect to symmetric conditions. This command combines all the previous commands to do so (see Algorithm 2).
Now the next code deals with the index rules that is able to handle Hooke’s tensor conversion from the fourth-order to the second-order tensor or second-order Voigt form using indexrule (see Algorithm 3).
Afterwards, we have a crucial stage that concerns with the transformation of Voigt’s notations into a fourth-rank tensor and vice versa. This stage includes the codes like HookeVto4, Hooke4toV. Also the codes HookeVto2 and Hook2toV transform the Voigt notations into the second-order tensor notations and vice versa. Finally the code in Algorithm 4 provides us a splendid form of fourth-order elasticity tensor as a six-by-six matrix in MCN notations.
We have presented here a minimal code of mathematica developed by . However, a numerical demonstration of this code has also been presented by  and the reader(s), interested in understanding this code in full, are requested to refer to the website already mentioned above.
Let us now proceed to Section 5 which in our view is the soul of the present work. In this section, we shall deal with the eigenvalues, eigenvectors, the nominal average of eigenvectors and the average eigenvectors, and so forth, for different species of softwoods, hardwoods and some specimen of cancellous bone.
5. Analysis of Known Values of Elastic Constants of Woods and Cancellous Bone Using MAPLE
Of course, “MAPLE” is a sophisticated CAS and produced under the results of over 30 years of cutting-edge research and development, which assists us in analyzing, exploring, visualizing, and manipulating almost every mathematical problem. Having more than 500 functions, this CAS offers broadness, depth, and performance to handle every phenomenon of mathematics. In a nutshell, it is a CAS that offers high performance mathematics capabilities with integrated numeric and symbolic computations.
However, in the present study, we attempt to explore how this CAS handles complicated analysis regarding elastic constant data.
The elastic constant data of our interest, for hardwoods and softwoods, as calculated by Hearmon [47, 48] are tabulated in Tables 5 and 6.
For cancellous bone, we have the elastic constants data available for three specimens [2, 11]. These data are tabulated in Table 7.
Precisely speaking, to meet the requirements of proposed research, a MAPLE code consisting of almost 81 steps has been developed. This MAPLE code (in its minimal form) is appended to Appendix A, while its entire working mechanism is included in an electronic appendix (Appendix B) and can be accessed from http://dx.doi.org/10.1155/2014/487314. Particularly, as far as our intention concerns, analysis of elastic constants data regarding hardwoods, softwoods, and cancellous bone using MAPLE is consisting of the following steps.(i)Computating the eigenvalues and eigenvectors for all the 15 hardwood species, 8 softwood species, and 3 specimens of cancellous bone using MAPLE.(ii)Computing the nominal averages of eigenvectors and the average eigenvectors for all the 15 hardwood species.(iii)Computing the average eigenvalues for hardwood species.(iv)Computing the average elasticity matrices for all the 15 hardwood species.(v)Plotting the histograms for elasticity matrices of all the 15 hardwood species.(vi)Plotting the graphs for I, II, III, IV, V, and VI eigenvalues of 15 hardwood species against their apparent densities (see Figures 14–16; also see Figure 17 for a combined view).
5.1. Computing the Eigenvalues and Eigenvectors for Hardwoods Species, Softwood Species, and Cancellous Bone Using MAPLE
Here we present the outcomes of MAPLE code (which is mentioned in Appendix A) regarding the computation of eigenvalues and eigenvectors of the stiffness matrices concerning 15 hardwood species (see Table 5).
It is well known that the eigenvalues and eigenvectors of a stiffness matrix or its compliance matrix are determined from the equations  where is the 6 × 6 identity matrix and represents the normalized eigenvectors of or . Naturally, or being positive definite, therefore there would be six positive eigenvalues and these values are known as Kelvin’s moduli. These Kelvin’s moduli are denoted by and if possible are ordered in such a way that . Also, the eigenvalues of can simply de computed by taking the inverse of the eigenvalues of .
A MAPLE code mentioned in Step 16 of Appendix A enables us to simultaneously compute the eigenvalues and eigenvectors for all the 15 hardwood species. Also, the results of this MAPLE code are summarized in Tables 8 and 9. Further, by simply substituting the elasticity constants from Table 6 in Step 2 to Step 11 of Appendix A, and performing Step 16 again, one can have the eigenvalues and eigenvectors for the 8 softwood species. The computed results for softwood species are summarized in Table 10. Similarly, substitution of elastic constants data for cancellous bone from Table 7 into Step 2 to Step 11 and repeating Step 16 of the MAPLE code yields the desired results tabulated in Table 11.
5.2. Computing the Average Elasticity Tensors for Woods and Cancellous Bone Using MAPLE
In order to compute average elasticity tensors for the hardwoods, softwoods, and cancellous bone, let us go through a brief mathematical description regarding this issue , as this notion helps us developing an appropriate MAPLE code to have desired results.
Suppose, we have measurements of the elastic constants data for some particular species or material, and we want to construct a tensor representing an average elasticity tensor for that particular material. Then to do so, we need the following key algebra .(i)The eigenvalues , , and the eigenvectors , , for each of the th measurement of the particular species or material.(ii)The nominal average (NA), of the eigenvectors associated with the particular species. The nominal average is given by (iii)The average , which is obtained by the following formula: where stands for the inverse square root of the positive definite tensor ; that is, (iv)Now, we proceed to average the eigenvalue. For this, we need to transform each set of eigenvalues corresponding to each eigenvector to the average eigenvector ; that is, It is obvious that for all and , where is the nominal average of eigenvalues of material and is defined by (v)Eventually, we can now calculate the average elasticity matrix in such a way that Taking care of the above outlined steps, we have developed some MAPLE codes (See Step 17 to Step 81 of Appendix A) that enable us to calculate the nominal averages, average eigenvectors, average eigenvalues, and the average elasticity matrices for each of the 15 hardwood species tabulated in Table 5. In addition to this, Appendix A also retains few steps, for example, Steps 21, 26, 31, 36, 41, 46, 51, 56, 61, 66, 73, and 80, which are particularly developed to plot histograms of average elasticity matrices of each of the hardwood species as well as graphs between I, II, III, IV, V, and VI eigenvalues of all hardwood species and the apparent densities (see Table 5).
We summarize the above claimed consequences of each of the hardwood species one by one ss follows.
5.2.1. MAPLE Evaluations for Quipo
The nominal averages of eigenvectors for the two measurements (S. nos. 1 and 2 of Table 5) of Quipo are The average eigenvectors for Quipo areThe averages eigenvalue for the two measurements (S. nos. 1 and 2 of Table 5) of Quipo is 3.339500000.The average elasticity matrix for Quipo is
The histogram of the average elasticity matrix for Quipo is shown in Figure 2.
5.2.2. MAPLE Evaluations for White
The nominal averages of eigenvectors for the single measurements (S. no. 3 of Table 5) of White areThe average eigenvectors for White areThe average eigenvalue for the single measurement (S. no. 3 of Table 5) of White is 14.58799998. The average elasticity matrix for White is given as
The histogram of the average elasticity matrix for White is shown in Figure 3.
5.2.3. MAPLE Evaluations for Khaya
The nominal averages of eigenvectors for the single measurements (S. no. 4 of Table 5) of Khaya areThe average eigenvectors for Khaya are The average eigenvalue for the single measurement (S. no. 4 of Table 5) of Khaya is 16.15299998.The average elasticity matrix for Khaya is given as
The histogram of the average elasticity matrix for Khaya is shown in Figure 4.
5.2.4. MAPLE Evaluations for Mahogany
The nominal averages of eigenvectors for the two measurements (S. nos. 5 and 6 of Table 5) of Mahogany are
The average eigenvectors for Mahogany are The average eigenvalue for the two measurements (S. nos. 5 and 6 of Table 5) of Mahogany is 18.19400002.The average elasticity matrix for Mahogany is given as
The histogram of the average elasticity matrix for Mahogany is shown in Figure 5.
5.2.5. MAPLE Evaluations for S. Germ
The nominal averages of eigenvectors for the single measurement (S. no. 7 of Table 5) of S. Gem are