Vannessa Duarte C.^1.2 y Miguel Cerrolaza^2,3

1 Aula CIMNE-INABIO, Universidad Central de Venezuela. Caracas, Venezuela. E-mail: vannessa.duarte@inabio.edu.ve

2 Instituto Nacional de Bioingeniería, Universidad Central de Venezuela. Caracas, Venezuela.

3 Centro Internacional de Métodos Numéricos en Ingeniería, Universidad Politécnica de Cataluña

RESUMEN

Este trabajo se orientó en cómo las propiedades eléctricas del hueso influyen en el modelado óseo y se explicó cómo el proceso de mecano-regulación y la estimulación eléctrica contribuyeron con el crecimiento, curación y reparación del tejido óseo. En este estudio se usó el Método de Elementos de Contorno (MEC) para materiales piezoeléctricos y se modificaron las ecuaciones para incluir la densidad del hueso, simulando dos ejemplos numéricos: una diáfisis con una fractura y un ostión en progresión. Se analizó cómo se originan cargas eléctricas debido a la influencia mecánica y cómo el proceso de absorción ósea se generó en zonas con cargas positivas, mientras que con cargas negativas fueron asociadas con la formación de hueso. Se pudo concluir que la ausencia de cargas mecánicas hizo que no existan potenciales eléctricos y, a su vez, si el proceso de modelado no se llevó a cabo, el hueso perdió su capacidad de recuperación y su densidad.

Palabras clave: Piezoelectricidad, Modelado óseo, Elementos de contorno, Potencial eléctrico, Cargas mecánicas.

Piezoelectric bone modeling by using the boundary element method

ABSTRACT

Keywords: Piezoelectricity, Bone modeling, Boundary element method, Electric potential, Mechanical loads.

Recibido: febrero 2013 Recibido en forma final revisado: octubre 2013

INTRODUCCIÓN

La adaptabilidad de los huesos al ser sometidos a cargas mecánicas ha sido ampliamente estudiado por autores como Wolff (1892) y Martínez et al. (2006), cuyos trabajos tratan de modelar y simular el comportamiento del hueso ante solicitudes mecánicas y su proceso de modelado, remodelado o curación de fracturas. El mecanismo de control por el cual el proceso de modelado se lleva a cabo, es poco conocido, pero puede atribuirse en gran medida al efecto piezoeléctrico en el hueso (Marino & Becker, 1970).

La piezoelectricidad en el hueso fue descubierta por Fukada & Yasuda (1957), quienes dieron a conocer el mecanismo por el cual los esfuerzos mecánicos resultan en polarización eléctrica (efecto indirecto); y la aplicación de un campo eléctrico causa deformación (efecto converso) (Figura 1). La curación y crecimiento óseo son controlados por el porcentaje de deposición de hidroxiapatita que es atribuida al trabajo de las células óseas como los osteoblastos, éstos son atraídos por la diferencia eléctrica producida por la piezoelectricidad o debido a la formación de hueso (especialmente de colágeno) (Ramtani, 2008).

Los modelos teóricos asociados con el modelado de huesos han demostrado que el potencial eléctrico se genera cuando él mismo es sometido a cargas o esfuerzos mecánicos, dicho potencial puede influenciar la respuesta ante solicitudes funcionales (Norton & Moore, 1972); es así como Guzelsu & Demiray (1979), Marino & Becker (1970) y Norton & Moore (1972) han postulado que el potencial negativo generado por la deformación mecánica en la superficie del hueso mejora la formación de nuevo tejido óseo, contribuyendo en el proceso de curación y reparación.

Este trabajo se basa principalmente en un modelo electromecánico presentado por Fernández et al. (2012) quienes consideraron un hueso isótropo elástico asumiendo que el coeficiente de Poisson ν es independiente de la densidad aparente ρ caracterizando el comportamiento óseo y garantizando que el campo eléctrico se incremente con la densidad. A su vez, nuestro modelo introduce la utilización de cualquier grado de anisotropía para el estudio del modelado óseo y utiliza para ello el método de elementos de contorno cuyas ecuaciones son resueltas usando la transformada de Radon y muestra como el método resulta ser una herramienta de gran utilidad en el estudio de procesos biológicos.

MATERIALES Y MÉTODOS

Las propiedades de los tejidos biológicos se consideran elásticas, eléctricas, dieléctricas y de comportamiento electromecánico que dependen de una densidad ósea variable que permite predecir la resistencia (Guzelsu & Demiray, 1979; Sierpowska et al. 2003)

La caracterización de las propiedades mecánicas del hueso vienen dadas por el coeficiente de Lame cuyo arreglo matricial es similar al de un material transversalmente isótropo. De acuerdo con Weinans et al. (1992), el módulo de Young y la densidad aparente son tomados como:

donde: M y β son constantes positivas que caracterizan el comportamiento del hueso y cuyo valor son tomados como M=3790, β = 3 y el coeficiente de Poisson ν=0,3.

La estructura ósea como órgano vivo es complicada ya que 2/3 de la masa del tejido óseo compacto lo constituye material inorgánico (hidroxiapatita) y el resto del hueso consta de material orgánico, principalmente colágeno (compuesto macromolecular, proteína fibrosa que posee alta elasticidad). Debido a esta composición y tomando en cuenta de que el hueso se comporta similar a un cristal con simetría hexagonal o trigonal (que cristalográficamente se describe mediante una celda hexagonal), donde el eje polar es el eje del hexágono, la matriz de coeficientes piezoeléctricos resultante tienen muchos coeficientes nulos y de los no nulos sólo se toman los valores independientes (Gjelsvik, 1973).

A pesar que la hidroxiapatita presenta simetría hexagonal no exhibe respuesta piezoeléctrica; sin embargo, el material compuesto formado por las estructuras fibrosas de colágeno y la inclusión de cristales de hidroxiapatita sí presenta dicha propiedad (Rodríguez et al. 1995).

En el material óseo el eje de simetría juega un papel importante en el subsecuente análisis y por conveniencia es referenciado como la dirección del material (Figura 2) (Gjelsvik, 1973) Si el eje coordenado es colocado en el eje z coincidente con la dirección del material se consideran solo 4 constantes y la matriz piezoeléctrica toma la siguiente forma:

Modelo piezoeléctrico

El efecto piezoeléctrico es una de las principales causas de la actividad bioeléctrica del tejido óseo, comúnmente descrito como el acople constitutivo de los campos mecánicos y eléctricos, cuyo vector general de desplazamiento U_k está constituido por los desplazamiento u_k y el potencial eléctrico φ. La relación constitutiva (Ecuación (4)) consiste en las ecuaciones de movimiento y carga, cuyos coeficientes C_ijkl, e_ikl y ε_il son las constantes elásticas, piezoeléctricas y dieléctricas respectivamente, σ es el tensor de esfuerzo, γ el tensor de deformación, E el campo eléctrico y D es el desplazamiento eléctrico (Denda & Wang, 2009).

Para garantizar que el campo eléctrico se incremente con la densidad del hueso se introduce dentro del modelo una extensión de la dependencia electro-mecánica clásica, regulando el acople entre los campos mecánicos y eléctricos por una función (ρ ⁄ ρ* )^ζ (Fernández et al. (2012).

La relación de desplazamiento-deformación y campo eléctrico-potencial viene dada por la Ecuación (5):

Las variables elásticas y eléctricas son combinadas en una simple ecuación constitutiva (Ecuación (6)) por la introducción de índices en letras mayúsculas en un rango de 1 a 4 e índices en letras minúsculas en un rango de 1 a 3.

Introduciendo una notación matricial compacta que consiste en el reemplazo de pares de índices ij o kl por índices simples p o q, donde i,j,k,l pueden tomar valores de 1 a 4 y p o q toman valores de 1 a 9 (Yang 2005), la ecuación constitutiva piezoeléctrica (Ecuación (6)) puede expresarse en forma matricial (Ecuación (9))

Formulación de Elementos de Contorno

Moviendo el punto de carga ξ al contorno Γ de un cuerpo piezoeléctrico tridimensional, se obtiene la ecuación integral de contorno piezoeléctrica (Ecuación (10))

donde: C_KJ(ξ) denota el coeficiente de término libre en ξÎΓ que contiene la contribución de la parte infinita y la singularidad fuerte de la solución fundamental de tracción en la integral de contorno. La solución fundamental es obtenida resolviendo la ecuación en diferencias parciales (Ecuación (11)),

La ecuación integral de contorno (Ecuación (10)) es el punto de inicio para la discretización numérica. Para ello, el contorno Γ es subdividido en un número finito de elementos de contorno Γ^(e), que son asignados dentro de elementos bidimensionales de un problema tridimensional como sigue:

donde: X_n^(e) denota los valores nodales, N_n son las funciones de forma elementales, L^(e) es el número de nodos locales del elemento (e), y ξ, η son las coordenadas locales del elemento (Beer (2001). En este trabajo se utilizaron elementos lineales isoparamétricos: interpolación lineal en la geometría e interpolación lineal en la variable (desplazamientos, potencial, tracción, superficie de carga).

Para solucionar el sistema de ecuaciones usando una aproximación y aplicando la Ecuación (10) sucesivamente por cada nodo del contorno discretizado, se obtiene un sistema de ecuaciones de esta manera:

donde: {u} y {t} contienen los valores nodales de las variables tanto elásticas como eléctricas.

Soluciones Fundamentales para un cuerpo piezoeléctrico

Una solución similar a una anisótropa elástica fue presentada por Gaul et al. (2003) y Thoeni (2005) la cual puede ser extendida a piezoelectricidad estática usando la transformada de Radon como una representación integral de la solución fundamental. Introduciendo la función tensorial:

donde: a_i, b_l, como vectores arbitrarios, la solución fundamental viene dada por:

El contorno de integración de la Ecuación (17) ha sido desarrollado numéricamente por la intersección de la esfera unitaria con el plano Z_ir_i^o = 0 (Figura 3), que es el argumento de la distribución de Dirac.

Con la finalidad de calcular la solución fundamental para las tracciones T_MJ*, la primera derivada de la solución fundamental de la Ecuación (17) es necesaria,

con los tensores:

Finalmente, la solución fundamental T_MJ* puede ser calculada por:

Evaluación de las soluciones fundamentales piezoeléctricas anisótropas

En coordenadas esféricas (r,θ₁,θ₂) son definidos por:

donde: el vector unitario r_i⁰ sólo depende de los ángulos θ₁ y θ₂ y el vector z_i está expresado a lo largo del círculo unitario sobre el cual la integración tiene lugar. Se determinan los dos vectores mutuamente ortogonales sobre el plano α_i y β_i, entonces α_i, β_i y r_i⁰ forman las bases del sistema ortonormal para lo cual se escogen:

y el vector β_i es dado por:

Con esto, se obtiene z_i:

donde: los ángulos θ₁ y θ₂ son constantes durante la integración sobre f.

Entonces se puede expresar la función del tensor dada (Ecuación (16)) en términos de f. La evaluación de la solución fundamental se obtiene de:

Escena de interpolación para las soluciones fundamentales piezoeléctricas

Ya que G_MK^u (θ₁,θ₂) son funciones de sólo dos variables; el ángulo polar 0 ≤ θ₁ ≤ π y el ángulo ortogonal 0 ≤ θ₂≤ 2π y todos los posibles valores de θ₁ y θ₂ son localizados dentro del área rectangular del plano θ₂ θ₁, como se muestra en la Figura 4, se subdivide el área dentro de una retícula rectangular y se almacenan los valores calculados de GMKu en cada punto de la retícula θ₁^m y θ₂^m en el correspondiente arreglo. De esta manera, cualquier valor de G_MK^u arbitrario θ₁,θ₂ puede ser obtenido por interpolación lagrangiana de los valores almacenados.

Múltiples Regiones

El contorno de cada región (∂Ω) se puede descomponer en Γ_j y Γ_I, que tales Γ_j υ Γ_I = ∂Ω_j, j=1,2, donde: Γ_I es la interface del contorno correspondiente a Ω₁ y Ω₂.

Luego el sistema de ecuaciones para cada región viene dado por:

Las condiciones de continuidad y de equilibrio sobre la interface Γ_I pueden ser expresadas como:

Al introducir estas condiciones en la Ecuación (27) se puede obtener el siguiente sistema de ecuaciones sobre el dominio Ω.

Por último al aplicar las condiciones de borde sobre la Ecuación (29) se obtiene un sistema de ecuaciones algebraico lineal (A_u=b) (Beer, 2001).

La formación y reabsorción ósea puede deberse a las cargas eléctricas en la superficie del hueso, originadas por las contribuciones producidas por la carga mecánica Dε. Para mostrar numéricamente este fenómeno y entender el comportamiento mecánico, se han empleado los ejemplos de Fernández et al. (2012), una diáfisis de un hueso largo y un ostión con cargas a compresión y torsión utilizando el método de elementos de contorno.

El código empleado está basado en el MEC desarrollado por Beer (2001) en lenguaje Fortran y modificado con la adición de algoritmos que permiten simular materiales con cualquier grado de anisotropía, piezoelectricidad utilizando el concepto estático y multirregión para resolver el sistema de ecuaciones presentado.

Las constantes del material están dadas de acuerdo con Fotiadis et al. (1999), considerando los coeficientes piezoeléctricos (C/mm²) y dieléctricos (C²/N mm²) como sigue:

e₃₁= 1,50x10^-9 ,       e₃₃=1,87x10^-9,

e₁₅= 3,57x10^-9,        e₁₄=17,88x10^-9,

ε₁₁= 88,54x10^-12      ε₃₃=106,24x10^-12

Diáfisis de un hueso largo con una fractura

Se considera un segmento óseo (Diáfisis) fracturado a lo largo en una posición angular de 10º para su curación, en dirección vertical. La diáfisis tiene un diámetro 27 mm y largo de 150 mm, tomando en cuenta sólo la mitad de diáfisis a lo largo debido a la simetría presente.

Se presenta una parte interna de hueso esponjoso o trabecular con una densidad de 0,8 g/cm³ (Fernández et al. 2012) y dos segmentos de hueso cortical con densidad mayor de 1,6 g/cm³. Aplicándose compresión de 2,5 MPa, y con restricción de movimiento y potencial con un valor de cero en la parte inferior del hueso (Figura 6).

En esta situación ideal, el hueso tiende a fortalecerse y la formación del mismo ocurre en las superficies cóncavas mientras que se reabsorbe en las zonas convexas.

Para la simulación se analizaron los valores de los nodos ubicados en la interfaz de las regiones a diferente tamaño de retícula (Tabla 1), en donde por razones de gasto computacional, almacenamiento y calidad de los resultados se escogió el uso de 100x100 ya que presenta un resultado acorde para un material isótropo, cambiando solo el tamaño de la retícula a mayor grado de anisotropía.

Los valores de desplazamiento, potencial, tracciones y superficie eléctrica son obtenidos de la interpolación langrangiana para cada nodo ubicado dentro del plano θ₁ θ₂ (Figura 7 y 8), tanto para hueso compacto como esponjoso en una malla de 637 elementos cuadriláteros con dos regiones heterogéneas y elásticamente isótropas.

Ostión en progresión

El ostión o sistema Harversiano es la unidad funcional fundamental del hueso compacto cuya densidad es de 1,6 g/cm3 y presenta una estructura cilíndrica fuerte de 700μm de diámetro externo con un túnel cilíndrico y un corte esférico cóncavo de 200μm, que representa el área de reabsorción de la cavidad del hueso cortical.

El ostión es sometido a compresión de 20MPa a lo largo de la dirección longitudinal y torsión de 0,0084 Nmm de un máximo de carga durante un ciclo de caminata (Figura 9).

Se fija el potencial a cero en todo el contorno externo y el desplazamiento se restringe en todo el área inferior a cero.

Para la simulación se realizó un análisis de convergencia de malla (Figura 10), comenzando con una malla burda de 58 elementos hasta una malla refinada de 2328 elementos. En la Figura 10 se puede apreciar como la diferencia en los valores de desplazamientos de nodos en la parte superior (N₁), en medio del cono interno (N₂), y en el borde interno inferior (N₃), entre la malla de 190 nodos y las más finas, es poco significativo convergiendo a medida de que el número de nodos aumenta, pudiendo tomar la malla de 190 nodos por la reducción en tiempo de computo y la obtención de un buen resultado. Sin embargo tomando en cuenta el potencial eléctrico, la diferencia de resultados hace que la malla de 770 nodos presente un valor bastante acorde al compararlo con la más fina. Por lo anterior, en esta simulación se usó una malla de contorno del ostión con 768 elementos cuadriláteros y 770 nodos y una retícula de 100x100 para la evaluación e interpolación de las soluciones fundamentales (Figura 7).

En cuanto al efecto de borde se hicieron cálculos variando las condiciones en ±0,0001; ±0,001;± 0,1 en el ostión y el resultado ante estas variaciones no superaba el umbral de error de 1%, por lo cual el problema es estable ante variaciones de las condiciones de borde.

RESULTADOS

Diáfisis de un hueso largo con una fractura

Para evidenciar la influencia que tiene las condiciones mecánicas sobre la carga eléctrica producida en el hueso, en la Figura 11 se muestra la sección media de la diáfisis en la cual se puede notar superficie de carga eléctrica negativa asociada con procesos de formación de hueso (las células óseas son atraídas por la electronegatividad) y una absorción ósea en áreas con carga positiva.

Por la inclinación del segmento óseo perpendicular al eje z (eje de formación de hueso), las zonas que sufren compresión (partes cóncavas de un hueso) se cargan negativamente, y en ellas se deposita tejido óseo, y en las zonas que sufren tensión (partes convexas) se generan cargas positivas, y en ellas se reabsorbe el tejido óseo (Figura 12).

Debido a que el vector normal a la superficie es perpendicular a la dirección del eje z, los desplazamientos en esa dirección son afectados por el factor eléctrico en la superficie. Las Figuras 12 y 13, muestran cómo el potencial eléctrico y el desplazamiento mecánico están relacionados, por lo cual la estimulación eléctrica adicional ayuda a la deformación del hueso contribuyendo a la formación del mismo y/o disminuir el tiempo de curación después de una fractura.

Ostión en progresión

Para este ejemplo lo primero que se hizo fue probar como el uso de MEC sobre otros métodos, en este caso Método de Elementos Finitos (MEF) puede ofrecer una herramienta de cálculo y simulación de procesos biológicos adecuada y que a su vez el gasto computacional sea menor al solo necesitar la discretización del contorno de la pieza a estudiar. La Figura 14 muestra (a) el resultado de MEF obtenido de un software comercial (Comsol Multiphysics 3.4) y (b) el obtenido del código computacional desarrollado en Fortran. En MEF se usó una malla de aproximadamente 2300 elementos cuadriláteros y en MEC una de tan solo 770, obteniendo como resultado una diferencia de 0,06; que en sólidos complejos como un hueso largo (fémur) podemos considerar despreciable.

El remodelado óseo ocurre en grupos locales de osteoblastos y osteoclastos lo que se llama unidades multicelulares de hueso (BMUs por sus siglas en inglés), donde cada unidad se organiza dentro de un cono hueco de osteoclastos que reabsorben hueso seguido de osteoblastos que lo rellenan. Esta conducta se puede observar en un comportamiento electromecánico en el cual se puede apreciar la manera como el hueso responde a potenciales negativos que hacen que las células migren a zonas de depósito óseo y consuman hueso en zonas positivas (Figura 15 (a)), contribuyendo al proceso de modelado del hueso; así mismo, se puede notar como una carga adicional aplicada de aproximadamente 3μV en la superficie del ostión hace que el potencial aumente lo que también lleva a un incremento en el desplazamiento mecánico.

Por otra parte, el mecanismo por el cual las fuerzas que actúan sobre el hueso determinan la formación o la reabsorción del tejido óseo en una zona determinada, se debe a la generación de potenciales eléctricos, los cuales se forman en el componente orgánico del hueso, como consecuencia de las tensiones a la que es sometido (Becker & Marino, 1982) El colágeno tiene propiedades piezoeléctricas, y al igual que determinados cristales, es capaz de transformar las fuerzas que lo deforman en impulsos eléctricos, de este modo, el arqueamiento por tracción o presión, que sufren las largas fibras de colágeno, genera pequeños potenciales eléctricos. En la Figura 15 se puede observar que los potenciales eléctricos a largo plazo pueden hacer que se llenen espacios de hueso faltante y adapten su forma a las fuerzas ocasionadas por un ciclo de caminata normal.

DISCUSIÓN

Para resolver el problema utilizando MEC se recurrió al uso de la transformada de Radon y un esquema de retícula que satisface el problema de la integral de dominio que queda en la integral de contorno piezoeléctrica, el tamaño de dicha retícula debe ser estudiada para cada caso de anisotropía de modo de obtener resultados acordes a los ofrecidos por otros métodos numéricos.

El estudio se basó en el modelado de huesos para un instante de tiempo y la diferencia de densidades de diferentes tipos de hueso, ya que con ello se puede obtener una función estimulo que permita estudiar el proceso de remodelado óseo y la influencia de superficie de carga eléctrica en la curación de fracturas logrando predecir la cantidad de estimulación eléctrica necesaria para generar mayor formación de hueso y reducir el tiempo de recuperación.

Ante la ausencia de estímulos mecánicos. no se obtiene una estimulación eléctrica que es la responsable de la migración celular y cuya función es el modelado de la estructura ósea, ocasionando el debilitamiento y la capacidad de recuperación que poseen los huesos. Sin embargo, con la aplicación de un estímulo eléctrico se puede lograr que se llegue a un modelado casi normal del hueso conservando en cierta medida la densidad ósea.

AGRADECIMIENTOS

Está investigación se realizó gracias al apoyo de la Academia de Ciencias Físicas, Matemáticas y Naturales de Venezuela y el proyecto ALFA-ELBENET de la comunidad Europea.

REFERENCIAS

1. Becker, R. O. & Marino, A. (1982). Electromagnetism & Life. State University of NY Press.        [ Links ]

2. Beer, G. (2001). Programming the Boundary Element Method. Wiley.        [ Links ]

3. Denda, M. & Wang, C. (2009). 3D BEM for the General Piezoelectric Solids. Comput. Methods Appl. Mech. Engrg, 198, 2950-2963.        [ Links ]

4. Fernández, J., García-Aznar, J., Martínez, R. (2012). Piezoelectricity could predict sites of formation/resorption in bone remodelling and modelling. J Theor Biol (292), 86-92.        [ Links ]

5. Fotiadis, D. I., Foutsitzi , G., Massalas, C. V. (1999). Wave propagation modeling in human long bones. Acta Mechanica (137), 65- 81.        [ Links ]

6. Fukada, E., & Yasuda, I. (1957). On the piezoelectric effect of bone. J Phy Soc Japan, 12 (10), 1158-1162.        [ Links ]

7. Gaul, L., Kögl, M., Wagner, M. (2003). Boundary Element Methods for Engineers and Scientists. Springer.        [ Links ]

8. Gjelsvik, A. (1973). Bone remodeling and piezoelectricity I. J. Biomechanics, 6, 69-77.        [ Links ]

9. Guzelsu, N. & Demiray, H. (1979). Recent Advances: Electromechanical Properties and related models of bone tissues. Int. J. Engng Sci, 17, 813-851.        [ Links ]

10. Marino, A. & Becker, R. (1970). Piezoelectric Effect and Growth Control. Nature, 228, 473-474.        [ Links ]

11. Martínez, G., García, J., Doblaré, M., Cerrolaza, M. (2006). External bone remodeling through boundary elements and damage mechanics. J. of Math. and Comp in Simul. (73), 183-199.        [ Links ]

12. Norton, L. A. & Moore, R. R. (1972). Bone Growth in Organ Culture Modified by. Journal of Dental Research, 51 (5), 1492-1499.        [ Links ]

13. Ramtani, S. (2008). Electro-mechanics of bone remodelling. Int J Eng Sci, 1173-1182.        [ Links ]

14. Rodríguez, R., Otero, J., Bravo, J., Rivalta, M. (1995). Vibraciones torsionales de un cilindro piezoeléctrico. Revista Mexicana de Física, 41 (6), 865-874.        [ Links ]

15. Sierpowska, J., Töyräs, J., Hakulinen, M. A., Saarakkala, S., Jurvelin, J. S., Lappalainen, R. (2003). Electrical and dielectric properties of bovine trabecular bone—relationships with mechanical properties and mineral density. Phys Med Biol, 48, 775-786.        [ Links ]

16. Telega, J., & Wojnar, R. (2002). Piezoelectric Effect in Biological Tissues. J Theor Appl Mech. 3 (40), 723-759.        [ Links ]

17. Thoeni, K. (2005). Effiziente Berechnung anisotroper Fundamentallösungen für die Methode der Randelemente. Master Tesis. Graz University of Technology, Graz, Austria.

18. Weinans, H., Huiskes, R., Grootenboer, H. (1992). The Behavior of Adaptive Bone-Remodeling Simulation Models. J. Biomechanics, 25 (12), 1425-1441.        [ Links ]

19. Wolff, J. (1892). Das Gesetz der Transformation der Knochen. A. Hirsehwold.

20. Yang, J. (2005). An introduction to the theory of piezoelectricity (9). New York: Springer.        [ Links ] ]]> 1 1982 2 2001 3 2009 198 2950-2963 4 2012 (292) 86-92 5 1999 6 1957 12 10 10 1158-1162 7 2003 8 1973 6 69-77 9 1979 17 813-851 10 1970 228 473-474 11 2006 of Math. and Comp in Simul. 73 73 183-199 12 1972 51 5 5 1492-1499 13 2008 14 1995 41 6 6 865-874 15 2003 48 775-786 16 2002 3 40 40 723-759 18 1992 25 12 12 1425-1441 20 2005