Application of the Stiffness-Force method in the analysis of the fiber-matrix adherence in composite laminates

I. Adarraga, N. Insausti, J. Urruzola, F. Mujika

GRUPO MATERIALES + TECNOLOGÍAS / MECÁNICA DE MATERIALES

Departamento de Ingeniería Mecánica

Escuela de Ingeniería de Gipuzkoa

UPV/EHU

Aplicación del método Rigidez-Fuerza en el análisis de adherencia fibra-matriz en laminados composites
Historia del artículo: Recibido 5 de Mayo 2017 En la versión revisada 5 de Mayo 2017 Aceptado 31 de Mayo 2017 Accesible online 21 de Junio 2017	En el presente trabajo se ha realizado el análisis de la adherencia de una fibra a tracción embebida en la matriz mediante el método Rigidez-Fuerza (SFM) aplicado a los elementos finitos. En la formulación del método SFM se tienen en cuenta las ecuaciones de equilibrio y las de compatibilidad de las deformaciones, obtenidas éstas a partir de los principios del trabajo virtual y del trabajo virtual complementario. Para el análisis por elementos finitos mediante el método SFM se ha realizado una subrutina de usuario UEL de un elemento cuadrilátero de cuatro nodos axisimétrico, y se ha implementado en el programa comercial Abaqus. Primero se ha realizado un análisis mediante el programa Abaqus que utiliza el método de la rigidez en su formulación. A continuación se ha analizado el mismo modelo mediante la subrutina de usuario UEL que utiliza el método SFM. Finalmente se han comparado los resultados de ambos análisis con diferentes grados de discretización.
Palabras clave: SFM Elementos finitos Tracción
Application of the Stiffness-Force method in the analysis of the fiber-matrix adherence in composite laminates
Keywords: SFM Finite elemens Traction	This work shows the analysis of the adherence of a fiber embedded in the matrix by the Stiffness-Force method (SFM) in Finite Elements. The formulation of the SFM method takes into account not only the equilibrium equations but also the compatibility conditions, obtained from the principles of the virtual work and complementary virtual work. For the finite element analysis using the SFM method a user subroutine (UEL) of a four node axisymmetric quadrilateral element has been formulated, and has been implemented in the commercial code Abaqus. In a first step the specimen has been analyzed by Abaqus, which uses the stiffness method in its formulation. Then the same model has been analized by the user subroutine which uses the SFM method. Finally, the results obtained by both analysis have been compared with varying degrees of discretization.

1 Introducción

El método de la rigidez (SM) en elementos finitos puede interpretarse como una forma aproximada de resolver las condiciones de equilibrio, derivadas éstas del principio de los trabajos virtuales. En el método de la rigidez la ecuación de gobierno es

${\textstyle \left[K\right]\left\{x\right\}=}$$\left\{P\right\}$ (1)

El método de las fuerzas integrado dual (IFMD) tiene la misma ecuación de gobierno de la ecuación (1), sin embargo, en este método además de las ecuaciones de equilibrio, se tienen en cuenta las ecuaciones de compatibilidad, derivadas éstas del principio del trabajo virtual complementario [1].

2 El método Rigidez-Fuerza (SFM)

Las ecuaciones del IFMD para un elemento continuo discretizado en elementos finitos se obtienen a partir del principio del trabajo virtual y del trabajo virtual complementario.

${\textstyle \int _{V}^{}{\left\{\delta \epsilon \right\}}^{T}\left\{\sigma \right\}dv=}$${\left\{\delta {a}^{i}\right\}}^{T}\left\{{P}^{i}\right\}$ (2)

${\textstyle \int _{V}^{}{\left\{\delta \sigma \right\}}^{T}\left\{\epsilon \right\}dv=}$${\left\{\delta {P}^{i}\right\}}^{T}\left\{{P}^{i}\right\}$ (3)

{aⁱ} y {Pⁱ} son los vectores de desplazamientos y de fuerzas nodales respectivamente. Con los desplazamientos nodales, se conocen los desplazamientos en cualquier punto del elemento [2]:

${\textstyle \left\{u\right\}=\left[N\right]\left\{{a}^{i}\right\}}$

(4)

donde [N] es la matriz de funciones de interpolación de desplazamientos A partir de los desplazamientos nodales se conocen también las deformaciones en el elemento,

${\textstyle \left\{\epsilon \right\}=}$$\left[L\right]\left[N\right]\left\{{a}^{i}\right\}=\left[B\right]\left\{{a}^{i}\right\}$ (5)

L] es la matriz de operadores diferenciales y [B] es la matriz de forma. En este método las tensiones y los desplazamientos se interpolan de forma independiente [3]:

${\textstyle \left\{\sigma \right\}=}$$\left[Y\right]\left\{{F}^{j}\right\}$ (6)

En el que [Y] es la matriz de funciones de interpolación de fuerzas que satisfacen las ecuaciones de equilibrio. El vector de tensiones se expresa en función del vector deformaciones mediante la matriz de coeficientes de flexibilidad del material [S]:

${\textstyle \left\{\sigma \right\}=}$$\left[S\right]\left\{\epsilon \right\}$ (7)

Sustituyendo las expresiones (5) y (6) en el principio del trabajo virtual (2) se obtiene la ecuación de equilibrio del elemento:

${\textstyle \left(\int _{V}^{}{\left[B\right]}^{T}\left[Y\right]dv\right)\left\{{F}^{j}\right\}=}$$\left\{{P}^{i}\right\}$ (8)

Donde la integral de volumen es la matriz de equilibrio del elemento ${\textstyle \left[E\right]=}$$\int _{v}^{}{\left[B\right]}^{T}\left[Y\right]dV$ . Así la ecuación de equilibrio (8)se escribe como:

${\textstyle \left[E\right]\left\{{F}^{j}\right\}=}$$\left\{{P}^{i}\right\}$ (9)

Ahora, sustituyendo las expresiones (6), (7) y (9) en el principio del trabajo virtual complementario (3), éste se escribe como:

${\textstyle \left(\int _{V}^{}{\left[Y\right]}^{T}\left[S\right]\left[Y\right]dv\right)\left\{{F}^{j}\right\}=}$$\left[{E}^{e}\right]\left\{{a}^{i}\right\}$

(10)

Donde la integral de volumen es la matriz de flexibilidad ${\textstyle \left[{G}^{e}\right]=}$$\int _{v}^{}{\left[Y\right]}^{T}\left[S\right]\left[Y\right]dV$ del elemento [4].

${\textstyle \left[{G}^{e}\right]\left\{{F}^{j}\right\}=}$$\left[{E}^{e}\right]\left\{{a}^{i}\right\}$ (11)

De la ecuación (11) se tiene el vector de fuerzas independientes del elemento {F^j}:

${\textstyle \left\{{F}^{j}\right\}={\left[{G}^{e}\right]}^{-1}\left[{E}^{e}\right]\left\{{a}^{i}\right\}}$

(12)

Multiplicando en ambos lados de dicha expresión por la matriz [E^e] de equilibrio, el primer término de la ecuación corresponde al vector de fuerzas externas aplicadas en los nodos, {Pⁱ},

${\textstyle \left\{{P}^{i}\right\}=\left[E\right]{\left[G\right]}^{-1}{\left[E\right]}^{T}\left\{{a}^{i}\right\}}$

(13)

Así tenemos la ecuación de gobierno IFMD [5] y la matriz de rigidez del elemento:

${\textstyle \left\{{P}^{i}\right\}=\left[K\right]\left\{{a}^{i}\right\}}$

(14)

${\textstyle \left[K\right]=\left[E\right]{\left[G\right]}^{-1}{\left[E\right]}^{T}}$

(15)

3 Elemento cuadrilátero axisimétrico

3.1 Formulación SFM

A continuación se va a aplicar el método SFM para formular el elemento cuadrilátero axisimétrico de 4 nodos de la figura 1.

Figura 1. Elemento cuadrilátero axisimétrico

Como se ha visto en el apartado anterior, para obtener la matriz de rigidez del elemento es necesario calcular previamente las matrices de equilibrio y de flexibilidad.

En la formulación de este elemento se van a utilizar las mismas funciones de interpolación lineales que en el cuadrilátero general [N_i] [6]:

${\textstyle \left[{N}_{i\,}\right]={\frac {1}{4}}\left(1+\xi {\xi }_{i}\right)\left(1+\right.}$$\left.\eta {\eta }_{i}\right)$

(16)

En el caso de un elemento axisimétrico la matriz de operadores diferenciales [L] es la siguiente:

$\left[L\right]=\left[{\begin{array}{cccc}{\frac {\partial }{\partial r}}&{\frac {1}{r}}&0&{\frac {\partial }{\partial z}}\\0&0&{\frac {\partial }{\partial z}}&{\frac {\partial }{\partial r}}\end{array}}\right]$

(17)

La matriz de deformación del elemento [B] se obtiene derivando la matriz de funciones de interpolación es la siguiente:

${\textstyle \left[B\right]={\frac {1}{4}}\left[{\begin{array}{cccccccc}-\left(1-\eta \right)&0&\left(1-\eta \right)&0&\left(1+\eta \right)&0&-\left(1+\eta \right)&0\\{\frac {\left(1-\xi \right)\left(1-\eta \right)}{r_{c}}}&0&{\frac {\left(1+\xi \right)\left(1-\eta \right)}{r_{c}}}&0&{\frac {\left(1+\xi \right)\left(1+\eta \right)}{r_{c}}}&0&{\frac {\left(1-\xi \right)\left(1+\eta \right)}{r_{c}}}&0\\0&-\left(1-\xi \right)&0&-\left(1+\xi \right)&0&\left(1+\xi \right)&0&\left(1-\xi \right)\\-\left(1-\xi \right)&-\left(1-\eta \right)&-\left(1+\xi \right)&\left(1-\eta \right)&\left(1+\xi \right)&\left(1+\eta \right)&\left(1-\xi \right)&-\left(1+\eta \right)\end{array}}\right]}$

(18)

3.3 Matriz de Equilibrio [E^e]

Para el caso de un cuadrilátero general axisimétrico la matriz de equilibrio se calcula mediante la siguiente expresión:

${\textstyle \left[E\right]=2\pi \int _{V}^{}{\left[B\right]}^{T}\,\left[Y\right]{r}_{C}\left|J\right|dV}$

(19)

Donde r_c es la coordenada de cada punto de integración y ${\textstyle \left|J\right|}$ el determinante del Jacobiano. En el caso de un cuadrilátero general la matriz de forma [B] varía con las coordenadas de cada punto, y como matriz de funciones de interpolación de fuerzas [Y] se ha tomado la siguiente expresión [7,8]:

${\textstyle \left[Y\right]=\left[{\begin{array}{ccccccc}1&0&0&J_{11}^{2}\eta &J_{21}^{2}\xi &0&0\\0&1&0&J_{12}^{2}\eta &J_{22}^{2}\xi &0&0\\0&0&0&0&0&1&J_{12}^{}\xi +J_{22}^{}\eta \\0&0&1&J_{11}J_{12}\eta &J_{21}J_{22}\xi &0&0\end{array}}\right]}$

(20)

3.4 Matriz de Flexibilidad [Ge]

La matriz de flexibilidad del elemento se calcula, en un cuadrilátero axisimétrico mediante la siguiente expresión:

${\textstyle \left[G\right]=2\pi \int _{V}^{}{\left[Y\right]}^{T}\,\left[S\right]\left[Y\right]{r}_{C}\left|J\right|dV}$

(21)

Donde [S] es la matriz de flexibilidad del material:

$\left[S\right]=\left[{\begin{array}{cccc}{\frac {1}{E}}&-{\frac {\nu }{E}}&-{\frac {\nu }{E}}&0\\-{\frac {\nu }{E}}&{\frac {1}{E}}&-{\frac {\nu }{E}}&0\\-{\frac {\nu }{E}}&-{\frac {\nu }{E}}&{\frac {1}{E}}&0\\0&0&0&{\frac {2\left(1+\nu \right)}{E}}\end{array}}\right]$

(22)

3.5 Matriz de Rigidez del elemento [K^e]

Conocidas las matrices de equilibrio y de flexibilidad del elemento, lse obtiene la matriz de rigidez del elemento:

${\textstyle \left[{K}^{e}\right]=\left[{E}^{e}\right]{\left[{G}^{e}\right]}^{-1}{\left[{E}^{e}\right]}^{T}}$

(23)

Una vez calculadas las matrices de rigidez de los elementos, la matriz de rigidez del sistema se obtiene ensamblando las matrices de rigidez de los elementos.

4 Análisis mediante SFM del ensayo pull-out de fibra única

El ensayo de tracción de fibra única es una técnica experimental muy utilizada para la caracterización de las propiedades de la interfase entre la fibra y la matriz en materiales compuestos [9]. Existen en la bibliografía aproximaciones analíticas [10,11], basadas en funciones de tensión que cumplen las condiciones de equilibrio, y aproximaciones numéricas [12,13] que, con diferentes tipos de elementos y altos grados de discretización, proporcionan la distribución de tensiones en la fibra y en la matriz.

En este trabajo, se ha analizado por elementos finitos mediante el método SFM el modelo axisimétrico representado en la figura 2 del ensayo pull-out de fibra única.

Figura 2. Ensayo pull-out de fibra única

Para el análisis mediante el método SFM se ha realizado una subrutina de usuario UEL de un elemento cuadrilátero de cuatro nodos axisimétrico, y se ha implementado en el programa comercial Abaqus.

El estado de cargas y las condiciones de contorno del modelo axisimétrico analizado se indican en la figura 3, que en un primer paso se correponden con el estado intacto o previo al despeque.

Figura 3. Modelo axiximétrico y condiciones de contorno.

En las figuras 4 y 5 se muestra la distribución de tensiones axiales σ^f_z y de cortadura σ^f_rz de la fibra en la unión fibra-matriz.

Figura 4. Distribución de tensiones σ_zen la unión fibra matriz.

Figura 5. Distribución de tensiones τ_rz en la unión fibra matriz.

5 Conclusiones

Se ha formulado mediante el método SFM un cuadrilátero general axisimétrico y se ha implementado en una subrutina de usuario (UEL) de Abaqus.

6 Agradecimientos

Los autores agradecen a la Universidad del País Vasco (UPV/EHU) la financiación del Grupo de Investigación Mecánica de Materiales GIU 16/51 en la convocatoria de 2016.

Referencias

[1] S.N. Patnaik et al. Integrated Force Method Solution to Indeterminate Structural Mechanics Problems. National Aeronautics and Space Administration (2004).

[2] E. Oñate. Cálculo de Estructuras por el Método de Elementos Finitos” Centro Internacional de Métodos Numéricos en Ingeniería (1995).

[3] Kaljevic, S.N. Patnaik, D.A. Hopkins. Three dimensional structural analysis by the integrated force method. Computers and Structures 59, 4,pag 691-706 (1996).

[4] J. Robinson, Integrated theory of finite element methods. Wiley and sons (1973).

[5] S.N.Patnaik, M.S. Nagraj, Analysis of continuum by the integrated force method. Computers and Structures, 26: 899-905. (1987).

[6] E. Oñate. Cálculo de Estructuras por el Método de Elementos Finitos. Centro Internacional de Métodos Numéricos en Ingeniería, Barcelona (1995).

[7] THH Pian, K Sumihara. Rational approach for assumed stress finite elements, Int J Numer Methods Eng. 20: 1685-1695 (1984).

[8] C. Zhang, D. Wang, J. Zhang, W. Feng, Q. Huang. On the equivalence of various hybrid finite elements and a new orthogonalization method for explicit element stiffness formulation, Finite Elements Anal.Des. 43: 321-332 (2007).

[9] E. Graciani, V. Mantič, F. París, J. Varnab. Numerical analysis of debond propagation in the single fibre fragmentation test. Composites Science and Technology. 69,15–16: 2514–2520 (2009).

[10] M.Y.Quek, C.Y. Yue. Axisymmetric stress distribution in the single filament pull-out test. Materials Science and Engineering, A198: 105-116 (1994).

[11] Y. Wang, L. Zhou, Z. Wang, H. Huang, L. Ye. Stress distributions in single shape memory alloy fiber composites. Material and Design:32: 3783-3789 (2011).

[12] W. Sun, F. Lin. Computer modeling and FEA simulation for composite single fiber pull-out. Journal of thermoplastic composite materials. 14: 327-343 (2001).

[13] V. Bheemreddy, K. Chandrashekhara, L.R. Dharani, G.E. Hilmas. Modeling of fiber pull-out in continuous fiber reinforced ceramic composites using finite element method and artificial neural networks. Computacional Materials Science. 79: 663-673 (2013)