INTRODUCTION

The development, growth, and remodeling of skeletal structures is a highly regulated process beginning with mesenchymal stem cells condensations in the early embryo and finishing with the homeostatic skeleton of the adult. It is widely accepted that both genetic and epigenetic

factors determine the final shape and strength of the skeleton, and many authors have specifically proposed an epigenetic role for mechanical forces (1), (2), (3). Equations have been proposed to describe how mechanical forces modulate growth where a mechanobiological remodeling rate is superimposed on a baseline biological growth rate (4).

In 1976 a rigorous mathematical bone remodeling model was proposed by Cowin and Hegedus which is so-called adaptive elasticity theory. In the adaptive elasticity theory, it is assumed that the rate of change in bone mass is correlated with the history of mechanical strain (5). Skalak et al. later formulated a continuum model of growth (6). Carter (7) proposed a bone remodeling model which was targeted to produce a homeostatic level of an effective stimulus and at the same time, Huiskes et al. proposed that the process of bone remodeling is aimed at producing a homeostatic value of the strain energy density in the tissue (8).

Bone growth models incorporating both biological and mechanobiological influences have been proposed by Van der Meulen and colleagues (9) for modeling the cross-sectional growth of long bones and by Stevens et al. for modeling endochondral growth using a finite element model (4). An alternative approach is that remodeling of bone is in response to microdamage, either to regulate microdamage to a homeostatic level or as the stimulus for the activation of the coupled responses of osteoclast and osteoblast cells (10), (11), (12).

Recently, Vahdati and co-workers incorporated a cellular accommodation effect into the Huiskes et al.'s semi-mechanistic model (13). They showed that the model is sensitive to temporal sequence in which loading is applied. In a very recent effort, Vahdati and Rouhi proposed another modification on the Huiskes et al's model (13). They included the effects of both microcracks and disuse on activation of resorption in one unifying formulation based on latest experimental findings besides considering cellular accommodation effect. These findings have not been published, yet. In 2004, Rouhi and co-workers proposed a modification on the adaptive elasticity theory by replacing volume fraction with free surface density in the constitutive equations (14). In another attempt, Rouhi et al. incorporated a microcrack factor in their first model and showed that not only mechanical stimuli, but also their rate and history are effective and at play in the bone remodeling process (15). Considering the great importance of bone resorption process in the osteoporotic cases, Rouhi et al. proposed a separate model of bone resoprtion by using mixture theory with chemical reactions. In their bi-phasic mixture model, they found that not only mechanical factors are at play in the resorption process, but also chemical and biological factors have crucial effects (16). Considering three different constituents of bone, i.e. bone matrix, bone fluid, and bone resorbing cells; Rouhi proposed a tri-phasic model of bone resorption using mixture theory with chemical reactions (17). It is concluded that rate of bone resorption is a function of apparent density of bone matrix and bone fluid, fluid velocity, momentum supply to the fluid phase, and internal energy densities of different constituents, in the former model.

Frost demonstrated that not only mechanoregulation diagram is not plateau in the lazy zone for growing bone, but also it has a considerable time dependent positive slope. This positive slope represents modeling process (Fig. 1) (18).

[FIGURE 1 OMITTED]

MATERIALS AND METHODS

It is widely accepted that immature bone structure can be altered via two major mechanisms: (1) Metabolic base growth, which is dependent on genetic patterns, nutrition, etc., (2) Mechanoregulatory mechanisms. The former is called Genetic Factor while the latter is called Epigenetic Factor (19). Generally speaking, it can be applicable superimposing Genetic to Epigenetic effects (20).

Computational simulations are used to investigate the mechanoregulation of bone growth. The bone growth model employed here was proposed by Van der Meulen et al. (9). It simulates the growth of the cross-section of a long bone from an embryonic bone collar to maturity, where the rate of bone apposition, [[rho].sub.t], is equal to the sum of the baseline biological rate, [[rho].sub.m], and the rate due to mechanobiological effects, [[rho].sub.m], as defined in Eq. (1):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The baseline biological growth rate is a decaying exponential function of time that reduces to approximately zero by the age of 18 years, defined as follows (9), (20):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [tau] = 3.6 years, [[rho].sub.b0] is the initial biological growth rate in [mu]gr.c[c.sup.-1] .[day.sup.-1], arisen from Eq. (3) in which initial r is taken to be 10 [mu]g.[day.sup.-1] in this model, where r is the rate of change in trabecular thickness, [S.sub.v] is the specific surface area and pt is the density of bone matrix which is assumed to be constant and equal to 2 g.c[c.sup.-1]. Figure 2 shows how r changes as an individual ages (9).

[FIGURE 2 OMITTED]

Figure 3 shown the specific surface area ([S.sub.v]), as a function of porosity (p) with coefficients in m[m.sup.2].m[m.sup.-3] as follows (3):

[S.sub.v] = 0.02876[p.sup.5] - 0.10104[p.sup.4] + 0.13396[p.sup.3] - 0.09304[p.sup.2] + 0.03226p (4)

[FIGURE 3 OMITTED]

The finite element stress analysis is used to calculate stress and strain within bone. The strain energy density, SED, can be expressed as follows:

[SED.sub.e] = [1/2] [[epsilon].sub.ij.sup.e][[sigma].sub.ij.sup.e] (5)

where [[sigma].sub.ij] and [[epsilon].sub.ij] are the stress and strain tensors, respectively and e represents the number of each element in the model. The resulted SED is assumed to act as the mechanical stimulus of the remodeling process. The bone remodeling set of equations was introduced by Ruimerman et al. in the following form (22):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

As stated earlier, since there are two mechanisms contributing in the bone adaptation process before maturation, i.e. Genetic and Epigenetic; we hypothesized here to superimpose their effects. Thus, adaptation equation will consist both growth rate, [[rho].sub.b]and the remodeling rate, [[rho].sub.m]. It is performed by inserting a linear alternative, instead of plateau, into the lazy zone of common remodeling models, appeared in Eq. (6). The superposition assumption takes place at mid point of lazy zone at which the magnitude of d.[[rho].sub.e]/dt is equal to the calculated growth rate from Eq. (3). This model also interpolates the growth rate over the lazy zone to determine the magnitude of d.[[rho].sub.e]/dt for any other point within the domain. Thus, updated governing equations will be as follows (Fig. 4):

[FIGURE 4 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

where [u.sub.R] is the lower or resorption SED threshold, uD is the upper or deposition SED threshold, [[rho].sub.e] represents elemental density and c is a growth rate constant. The values of [u.sub.R] and [u.sub.D] are age-dependent and will approach the mature values with linear approximate variation within 18 years. The value of [u.sub.R] starts from zero at birth and increases until reaches its mature value, i.e. 0.002 J m[m.sup.3] Based on Ruberg's study (23), [u.sub.R] is considered as 80 percent of midpoint SED in lazy zone, while [u.sub.D] is assumed to be 20 percent greater than midpoint SED in the lazy zone. By this assumption, [u.sub.D] is also a linear function of time during first 18 years of life.

The superposition process also alters the resorption and deposition zones. This is carried out via replacing common models' constant, i.e. c, with its alternative for growing bone (24), (25). As it is quite evident in Fig. 1, the proportionality constant c in immature bone is not the same as in mature bone. Thus, we have considered a linear approximate temporal variation, like the ones for lazy-zone thresholds, which starts from its initial value, i.e. c = 0.01 g.day.m[m.sup.-3]., and approaches its mature value, i.e. 0.02 g.day.m[m.sup.-3], at the age of 18 years (23). In other words, the constant c is, in fact, a time dependent coefficient which affects the rate of bone adaptation under the effect of external load placed on the bone. A higher value of c will lead to more expansion and a thicker bone; it also acts as a handle for controlling the convergence (20).

Concentrating on Eq. 7, it is obvious that the second term imposes a positive, age-decaying slope to mechanoregulation diagram over the lazy zone. Assuming that the maturity happens at the age of 18 years, this model converges to common remodeling models proposed by others (e.g. Jacobs (28), Carter et al.(19), Huiskes et al.(12)), for mature bone.

As a physical run-time constraint, the elemental densities were assumed not to outrage the following range (0.001, 1.74 g/cc). The lower bound is taken 0.001 instead of zero due to some computational considerations. Zero density will lead to zero stiffness, internal stress shielding, checkerboard configuration, and even singularity in global stiffness matrix of structure in finite element solution (26). Thus, the following constraint was set:

0.001[g/cc]<[[rho].sub.e]<1.74[g/cc] (8)

The Van der Muelen et al.'s model uses age-dependent loading cases to simulate the loading history for an individual bone (Fig. 5). This force was obtained from body weight during growth, and it is assumed that the force and age are linearly proportional. The load starts from approximately zero values and gradually increases with a mean slope of 1.11 Newton per year, until reaches 20 Newtons at the age of 18 years (27).

[FIGURE 5 OMITTED]

The relation between the elastic modulus and the bone density is assumed to be in the following form as proposed by Jacobs in 1994 (28):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

where [[rho].sub.e] refers to the density of each element. In addition to elastic modulus regulation, the model proposed here uses the following equation, also proposed by Jacobs (28), to update elemental Poisson's ratios:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

In this study, the simulation was performed utilizing a MATLAB code. Bone specimen was assumed a 2-dimensional 0.1 x 0.1 m[m.sup.2] square containing a uniformly distribution of bone sensor cells, i.e. osteocytes (1600 Osteocytes per cubic millimeters). The aforementioned square geometry was meshed with linear constant strain triangular elements. The minimum appropriate size of mesh was obtained 25 x 25 nodes from independency test, while independency test variable was defined by Eq. (11) as follows:

F = [2304.[sigma](e=1)][[epsilon].sub.ij.sup.e].[[sigma].sub.ij.sup.e] (11)

where, F is target function and e is the summation index number over elements. First and second degrees of freedom of first node at left-bottom corner and first degree of freedom of 25th node at right-bottom corner were restrained as minimum required boundary conditions. We also prevented overconstrained condition in applying boundary conditions. Applied boundary conditions are shown in Fig. 7.

[FIGURE 7 OMITTED]

These two factors also converge to zero at 25 nodes per side and higher.

Post-processing was performed to calculate stresses and strains at gauss integration points. Then, the state of remodeling, i.e. d.[[rho].sub.T]/dt, is calculated using Eq. 7. Finally, at the end of each iteration new elemental densities are updated with the following equation:

[[rho].sub.e](t + [delta]t, n) = [[rho].sub.e] (t, n) + [delta]t [[[sigma][[rho].sub.e](t, n)]/[delta]t] (12)

where [[rho].sub.e(t+[delta]t)] is the updated density of nth element which will lead to updated young's moduli and Poisson's ratios.

The iterative solution procedure includes two loops which one fulfils temporal simulation and the other stands for convergence testing.

The convergence loop is the inner and the temporal loop is the outer one. Number of temporal iterations was considered 75 which approximately reveals 20th year of age, assuming 365 days each year, and 100 days each iteration.

[FIGURE 6 OMITTED]

Ruberg defined volume-averaged change in bone volume fraction, e, as follows:

e = [[[integral] (v)][delta][v.sub.b]dV]/[[integral](v)]dV]] (13)

where, V represents geometry volume and [v.sub.b] shows bone volume fraction. Stoppage threshold of e was set between [10.sup.-4] and [10.sup.-5]. Considering [v.sub.b] = [V.sub.B]/[V.sub.T], in which [V.sub.B] is volume of bone tissue and [V.sub.T] represents total bone specimen volume, this formulation leads to Eq. (14). Algebraic procedure can be found in (23).

If d[[rho].sub e]/dt > [10.sup.-4] gr/cc.day then continue! (14)

where e represents an arbitrary element within the geometry. Fig.8 shows a general flowchart of our solution.

[FIGURE 8 OMITTED]

RESULTS AND DISCUSSIONS

Four different simulations were performed to determine whether our model shows reasonable behavior under various loading conditions. The conceptual initial architecture was considered as a uniform density distribution with a density of 0.7 g/cc (29). In the figures 9-12, final configurations of density distributions within our model can be seen. The density increases on contours from blue, the lowest, to red, the highest.

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

In the first simulation, we have applied a 1 Newton force on the upper side of square which was oriented 45[degrees] clockwise with respect to horizon (Fig. 9a).

Then, we changed the orientation of force to 45[degrees] counterclockwise, with same magnitude. Figures (9b-e) show the adaptation process due to change in external force orientation. It is obvious that hard tissue mass is gradually transferred from initial major strut to second major strut. It is strongly in agreement with Huiskes et al. findings in (12). Both final architectures show similar configurations and mass distributions (Fig. 9).

In the second test, we applied a uniformly distributed shear force on the right side of the square. As can be shown in Fig. 10, model has regulated itself by producing two struts in response to loading; first just beneath loaded nodes and second is developed with 45[degrees] angle oblique to the first. The first strut is developed in the direction of maximum shearing stress and the second is produced in the direction of principal direction of normal stresses. Based on Wolff's law, some researchers have interpreted that bone reinforces itself in orientation which is maximally stressed (7), (8), (28). Our 2nd simulation is in agreement with the former interpretation of Wolff's law.

In the third simulation, the bending condition was simulated by applying a bi-ramp on the upper side of the model. As can be shown in Fig. 11, the final configuration of this case has a strong qualitative correlation with distal femur macroscopic architecture. Asymmetricity found in final configuration of this loading case is interpreted in relation with asymmetric boundary conditions.

Trabecular thickness, which is one of the most important geometric parameters relevant to osteoporosis and also osteopenia (30), is the next issue we are going to discuss. In the forth and last qualifying round of tests, initially we have loaded the model with an oblique 1 Newton force. Then, we have over- and under-loaded the model by 25% of the original force (1N). Graphical presentation of such test can be seen in Fig. 12. As it is obvious in the Fig. 12, over-loading condition induces an elevation and under-loading or disuse, induces a reduction in trabeculae thicknesses. It is widely expressed by authors that the former (overloaded one) causes a hypertrophy, and the latter (the under-loaded and/or disused case) will lead to an atrophy and osteoporotic condition in a healthy bone(31), (32).

In this model, nothing was said about growth yet. We have introduced our growth included model (GIM) by Eq. 7. We also mentioned growth excluded model (GEM) in Eq. 6. After applying these two models on loading cases, specimens' masses were drawn as a function of age. All the diagrams show qualitative resemblance, but not numeric. We have included the diagrams of first test results for both GIM and GEM in Fig. 13, for instance. As the figure shows, both models start from initial mass in the order of 0.61 grams. As individual ages, the GIM accelerates more than GEM due to the growth term. This situation continues until the growth is stopped at the age of 18. The difference between GIM and GEM at the age of 18 is about 80mg. This difference represents approximately 9.3% of weight for homeostatic mature GIM in first loading case, which is 0.856 grams. The aforementioned difference can be considered as the error included in the model due to growth effect cancellation assumption.

[FIGURE 13 OMITTED]

Figure 14 contains a diagram which is calculated by subtracting mass of the growth excluded model (GEM) from the mass of growth included model (GIM) which is called as "impact of growth assumption" by the authors. Based on this figure, convergence takes place at the age of 18. This shows that bone mass calculated using the GIM will approach the one of derived employing the GEM after maturation. In addition to that, this diagram resembles a step response of a 2nd--order closed-loop control system which can be studied considering the magnitude, nature and physiologic interpretation of its particular System parameters, i.e. damping ratio and natural frequency. Some authors have interpreted system parameters of such system as a gauge of bone health and youth and its ability to response to remodeling stimulus (33).

[FIGURE 14 OMITTED]

The rate of change in bone mass diagram is also studied in this research (Fig. 15). Figure 15 shows the rate of change of bone mass as individual ages. The values of this parameter are 0.0436g/day and 0.0674g/day at birth for GEM and GIM, respectively. This difference is relevant to the initial rate of growth.

[FIGURE 15 OMITTED]

Figure 16, also shows temporal variation in the difference between GEM and GIM. The parameter decays exponentially with time and the numeric analysis shows that the time constant is equal to that of the growth rate equation, introduced in Eq. 2.

[FIGURE 16 OMITTED]

CONCLUSIONS

In this research, the authors are trying to put an emphasis on the point that growth and remodeling can be considered simultaneously to predict homeostatic structures of growing bone during the first two decades of life. It is well accepted that in a growing bone, or immature bone, both genetic and epigenetic factors are at play. In other words, both modeling and remodeling are active. Thus, in this study, the authors superimposed bone modeling equation on bone remodeling equation, in the hope of modeling immature bone adaptation.

From the simulations performed in this study, it is concluded that this model can produce strut- and platelike structures found in substructures of spongy bone; as was resulted in other studies (1), (23). All the homeostatic configurations show a fair agreement with Wolff's law (3); topologic optimization methods (34); bone micro-structure; and also other researchers' works (22), (35). This study also shows that neglecting growth effect in immature bone models can cause a noticeable error into the final results. In this study, neglecting growth effect in bone adaptation process could cause an error of about 9.3% in final bone mass.

It is evident that the biomechanics society is suffering from the lack of a mechanistic model of bone adaptation. One reason of this shortage comes from the fact that bone adaptation is a multidisciplinary problem which is related to different fields (e.g. cell biology; mechanics; and chemistry). One of the most important challenges in the field of bone adaptation is mechanotransduction of bone cells. In other words, understanding how mechanosensors, i.e. osteocytes, react to different mode of loadings. It is hoped that by closer collaboration among researchers from medicine and engineering more insight will be gained in the near future, and more lights will be shed on this very essential and complex problem of biomechanics.

ACKNOWLEDGEMENT

Special thanks to Amirkabir University of Technology and Azad University of Medical Sciences, Mashhad, Iran, and University of Ottawa, Canada.

REFERENCES

(1.) Pauwels F., 1980, Biomechanics of the locomotor apparatus. Berlin, Springer.

(2.) Thompson D.W., 1992, On growth and form. New York, Dover Publications.

(3.) Wolff J.L., 1892, The law of bone remodeling, Springer, Berlin.

(4.) Stevens S.S., Beaupre G.S. and Carter D.R., 1999, Computer model of endochondral growth and ossification in long bones: Biological and mechanobiological influences. J. Orthopedia 17:646-653

(5.) Cowin S.C. and Hegedus D.H., 1976, Bone remodeling I: a theory of adaptive elasticity. J. of Elasticity, 6: 313-326.

(6.) Skalak, R., Dasgupta G., Moss M.,Otten E., Dullemeijer P. and Vilmann H., 1982, Analytical description of growth. J. Theor. Biol. 94 (3): 555-577.

(7.) Carter D.R., 1987, Mechanical loading history and skeletal biology. Biomechanics, 20:1095-1099.

(8.) Huiskes R., Weinans H. and Grootenboer H.J., 1987, Adaptive bone remodeling theory applied to prosthetic design analysis. J. of Biomech., 20:1135-1150.

(9.) Van der Meulen, M.C.H., Beaupre G.S. and Carter D.R., 1993, Mechanobiologic influences in long bone cross-sectional growth. Bone 14 (4): 635-642.

(10.) Martin R.B., 2000, Towards a unifying theory of bone remodeling. Bone, 26 (1): 1-6.

(11.) Prendergast P.J. and Taylor D., 1994, Prediction of bone adaptation using damage accumulation. J. Biomech., 27 (8): 1067-1076.

(12.) Huiskes R., Ruimerman R., Van Lenthe G.H. and Janssen J.D., 2000, Effects of mechanical forces on maintenance and adaptation of form in trabecular bone, Nature, 405: 704-706.

(13.) Vahdati A., Rouhi G., Ghalichi F., and Tahani M., 2008, Mechanically induced trabecular bone remodeling including cellular accommodation effect: a computer simulation, The CSME Transactions, 32(2).

(14.) Rouhi G., Firoozbakhsh K., Epstein M., Herzog W. and Sudak L., 2004, Free surface density instead of volume fraction in the bone remodeling equation: theoretical considerations. Forma, 19 (3): 165-182.

(15.) Rouhi, G., Epstein M., Herzog W. and Sudak L. 2006, Free surface density and microdamage in the bone remodeling equation: theoretical considerations. Int. J. Eng. Sci. 44 (7): 456-469.

(16.) Rouhi G., Epstein M., Sudak L. and Herzog W., 2007, Modeling bone resorption using mixture theory with chemical reactions. J. of Mech. of Mat. and Struct., 2(6): 1141-1156.

(17.) Rouhi G., 2007, Theoretical investigations on bone resorption process: a tri-phasic mixture model, J. Mech. Beh. Biomed. Mat. (Accepted subject to revision).

(18.) Frost H.M., 1990, Skeletal structural adaptations to mechanical usage (SATMU): 2. Redefining Wolff's law: the remodeling problem. Anat. Rec. 226: 414-422.

(19.) Carter D.R., Mikic B. and Padian K., 1998, Epigenetic mechanical factors in the evolution of long bone epihyses. Zool. J. Linn. Soc. 123 (2): 163-178.

(20.) Niamh C., Patrick and Prendergast J., 2005, Evolution of mechanoregulation of bone growth will lead to non-optimal bone phenotypes. J. of Theor. Biology, 235: 408-418.

(21.) Martin R.B., 1994, Porosity and specific surface of bone. Chapter 3 in: Critical reviews in biomedical engineering, CRC Press 10: 179-222.

(22.) Waide V., Cristofolini L., Stolk J., Verdonschot N. and Toni A., 2003, Experimental investigation of bone remodeling using composite femurs, Clinical Biomechanics, 18: 523-536.

(23.) Ruberg T., 2004, PhD dissertation, Technische Universitat Braunschweig, Germany.

(24.) Uhthoff H.K. and Jaworski Z.F.G., 1978, Bone loss in response to long term immobilization. J. Bone Jt Surg., 60-B: 420-429.

(25.) Jaworski Z.G.F., Liskova-Kiar M. and Huthoff H.K., 1980, Effect of long term immobilization on the pattern of bone loss in older dogs. J Bone Jt Surg., 62-B: 104-110.

(26.) Zidi M. and Ramtani S., 2000, Stability analysis and finite element simulation of bone remodeling model, J. of Biomech. Eng., 122: 677-680.

(27.) McCammon R.W., 1970, Human growth and development, Thomas, Springfield, IL, USA.

(28.) Jacobs C.R., 1994, Numerical simulation of bone adaptation to mechanical loading, PhD Dissertation, Stanford University.

(29.) Janko D. Jovanovic, Miomir L.J.. Jovanovic, 2004, Biomechanical model of vertebra based on bone remodeling, Med. and Biol. 11(1): 35-39.

(30.) Guo X.E. and Kim C.H., 2002, Mechanical consequence of trabecular bone loss and its treatment: a three-dimensional model simulation, Bone, 30(2): 404-411.

(31.) Buckwalter J.A., Glimcher M.J., Cooper R.R. and Recker R., 1995, Bone Biology, Part I, Structure, blood supply, cells, matrix, and mineralization, J. of Bone and Joint Surg., 77: 1276-1275.

(32.) Buckwalter J.A., Glimcher MJ., Cooper R.R. and Recker R., 1995, Bone Biology, Part II, Formation, form, modeling, remodeling, and regulation of cell function, J. of Bone and Joint Surg., 77: 1276-1289.

(33.) Mullender M.G., Huiskes R. and Weinans H., 1994, A Physiological Approach to the Simulation of Bone Remodeling as Self Organizational Control Process, J. Biomech., 27(11): 1389-1394.

(34.) Sigmund O., 2001, A 99 line topology optimization code written in matlab. Struct. Multidisc.Optim. 21: 120-127.

(35.) Treharne R.W., 1981, Review of Wolff's law and its proposed means of operation. Orthop Rev, 10: 35-47.

(1) Seyyed Amir Hooshiar Ahmedi, (2) Gholamreza Rouhi, (1) Hamidreza Katouzian and (3) Seyyed Ali Hooshiar Ahmedi

(1) Department of Biomedical Engineering, Amirkabir University of Technology, Tehran, Iran

(2) Department of Mechanical Engineering, University of Ottawa, 161 Louis Pasteur, Room A017 Ottawa, Ontario, Canada K1N 6N5

(3) Azad University of Medical Sciences, Mashhad, Iran

Corresponding Author: Gholamreza Rouhi, Department of Mechanical Engineering, University of Ottawa, 161 Louis Pasteur, Ottawa, Ontario, Canada K1N 6N5, Tel: 001-(613) 562-5800/6268, Fax: 001-(613) 562-5177