Research Article  Open Access
A CAS Approach to Handle the Anisotropic Hooke’s Law for Cancellous Bone and Wood
Abstract
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 fourthorder 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.
1. Introduction
The study of material symmetry of 3dimensional 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 [1]. Of course, the familiarity with material symmetries is the best way to categorize the materials. However, according to [2], 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 [3], 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 [4], anisotropic Hooke’s law for cancellous bone and wood [2], the validation of generalized Hooke’s law for coronary arteries [5], constitutive relationships of fabric, density, and elastic properties in cancellous bone architecture [6], analysis of elastic anisotropy of wood material for engineering applications [7], a multidimensional anisotropic strength criterion based on Kelvin’s modes [8], spectral decomposition of a 4thorder covariance tensor: application to diffusion tensor MRI [9], a normal distribution of tensor valued random variables: applications to diffusion tensor MRI [10], 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 [11], relationships between morphology and bone elastic properties can be accurately quantified using highresolution computer reconstructions [12], and intrinsic mechanical properties of trabecular calcaneus determined by finite element models using 3D synchrotron microtomography [13], 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 [4], 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 [6]. 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 [7].
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 4thorder tensor. We would like to emphasize that the 4thorder 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 [22], biomechanics [13], tissue mechanics [1], geophysics [23], and so forth.
Though an adequate amount of research work has been dedicated towards the algebra of 4thorder 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 FourthOrder 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 secondorder tensor, a fourthorder tensor (say), is a linear map from to . That is, it assigns to each secondorder tensor the secondorder tensor ; that is, where is a vector space and is a linear space of dimension and is equivalent to Also, any 4thorder 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 4thorder tensors is given by the composition In abstract index notations, we have Moreover, the transpose of is a 4thorder 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 4thOrder Tensor
In case of a secondorder tensor, the only known symmetry represents an invariance under the mutual rotation of two indices, while for the 4thorder tensor, there are several notions of symmetry that represent invariance under exchanging pair of indices. Thus, for the 4thorder 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 [24].
2.1.1. Major Symmetry of a FourthRank Tensor
For , we say that is symmetric if , or in component form and we say that is skewsymmetric if or .
Moreover, we have the decomposition Then the set of all symmetric fourthorder 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 fourthorder 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 4thorder tensors that satisfy the minor symmetry is denoted by
2.1.3. Total Symmetry
A fourthorder 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 fourthrank tensors bearing total symmetry is denoted by Thus any fourthorder 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 fourthorder tensor are described as [25] Obviously, a fourthorder tensor is totally symmetric if or equivalently , where stands for a fourthorder 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 [25] has also shown that there is an isomorphism (i.e., a 11 linear map) between totally symmetric fourthrank tensor and a homogeneous polynomial of degree four. Also, if a totally symmetric fourthorder 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 fourthorder 3dimensional tensor with minor symmetries, as it is customary in the 3dimensional elasticity theory.
In the following section, we discuss the notion of 9dimensional representation of a 3dimensional fourthorder tensor and then particularize this representation to a fourthorder elasticity tensor having minor symmetries due to the wellknown anisotropic Hooke’s law.
3. A 9Dimensional Exposition of a 3Dimensional 4thOrder Tensor and the Anisotropic Hooke’s Law
In accordance with the evolution of tensor theory, the algebra of the 3dimensional 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 [30], in my view this tool still needs some enhancements to handle symbolic tensor data. Moreover, Constantinescu and Korsunsky [31] 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 3dimensional 4thorder tensor as a 9 × 9 matrix, that is, a 9dimensional representation of a 3dimensional fourthorder tensor. For this purpose, there is a basic notion of representing a fourthorder tensor as a secondorder tensor [22]. According to [22], if , , then we set Thus, is an orthonormal basis of .
Therefore, any fourthrank tensor can be assumed as an dimensional fourth rank tensor , or as an dimensional secondorder tensor , where , .
Again, we have from (8) that hence the twoway 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 eigenstiffnesseigenstrain 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 stressstrain symmetries do not affect the stiffness tensor , the number of components of stiffness tensor in 3dimensional space is equal to . Now, if we impose the symmetries of stressstrain 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 fourthorder 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 3dimensional 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 [32] and after that various successors of Voigt have used his fabulous notions for flattening of fourthorder tensors. But Lord Kelvin [40] 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, [39] have introduced some sophisticated methodology to unfold a fourthrank 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 SixDimensional 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 secondrank tensors, respectively, in 6dimensional 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. MehrabadiCowin Notations (MCN)
To preserve the tensor properties of Hooke’s law during the unfolding of fourthorder elasticity tensor, [39] 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), [39] have evoked that (40) can be rewritten as where the new sixdimensional stress and strain vectors denoted by and , respectively, are multiplied by the factors . Also, is a new sixbysix matrix [39]. 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 fourthorder elasticity tensor as a sixbysix matrix, a simple conversion table has been proposed by [39], 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 [30], Maple [41], Mathematica [42], wxMaxima [43], and some peculiar packages like Tensor Toolbox [26–29], GRTensor [44], Ricci [45], 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 userfriendly 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 [31]. 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 [31] 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:
SetDirectory[C:/Program Files/Wolfram/Research/Mathematica/8.0/AddOns/Packages].
Next is concerning the creation of tensor. The code in Algorithm 1 generates the stress and strain tensors using the index rules specified by [31].

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 fourthorder to the secondorder tensor or secondorder Voigt form using indexrule (see Algorithm 3).

Afterwards, we have a crucial stage that concerns with the transformation of Voigt’s notations into a fourthrank tensor and vice versa. This stage includes the codes like HookeVto4, Hooke4toV. Also the codes HookeVto2 and Hook2toV transform the Voigt notations into the secondorder tensor notations and vice versa. Finally the code in Algorithm 4 provides us a splendid form of fourthorder elasticity tensor as a sixbysix matrix in MCN notations.

We have presented here a minimal code of mathematica developed by [31]. However, a numerical demonstration of this code has also been presented by [31] 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 cuttingedge 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.
 
Source: this data is consulted from [47, 48]. The number of repetitions of the particular species in this table shows the number of measurements done for that particular wood species. The units of (bulk or apparent density) are g/cm^{3} and the units of elastic constants are GPa = dynes/cm^{2}. 
 
Source: this data is consulted from [47, 48]. The number of repetitions of the particular species in this table shows the number of measurements done for that particular wood species. The units of (bulk or apparent density) are g/cm^{3} and the units of elastic constants are GPa = dynes/cm^{2}. 
For cancellous bone, we have the elastic constants data available for three specimens [2, 11]. These data are tabulated in Table 7.
 
Source: the data for cancellous bone has been taken from [2, 11]. 
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 [3] 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.
