Artículos de Investigación

SIMULACIÓN NUMÉRICA DEL COMPORTAMIENTO ESTRUCTURAL DE PERFILES TUBULARES SOMETIDOS A FATIGA DE ULTRA-BAJO CICLAJE

NUMERICAL SIMULATION OF THE STRUCTURAL BEHAVIOR OF TUBULAR PROFILES SUBMITTED TO ULTRA-LOW CYCLES FATIGUE

Leonardo Lopez
Universidad Centroccidental Lisandro Alvarado, República Bolivariana de Venezuela
María Eugenia Marante
Universidad Centroccidental Lisandro Alvarado, República Bolivariana de Venezuela
Ricardo Picón
Universidad Católica de Temuco, República Bolivariana de Venezuela

SIMULACIÓN NUMÉRICA DEL COMPORTAMIENTO ESTRUCTURAL DE PERFILES TUBULARES SOMETIDOS A FATIGA DE ULTRA-BAJO CICLAJE

Gaceta Técnica, vol. 19, núm. 2, pp. 19-35, 2018

Universidad Centroccidental Lisandro Alvarado

La Revista Gaceta Técnica es de acceso libre, gratuito y abierto. Para la reproducción parcial o total de los trabajos o contenidos publicados, se exige reconocer los derechos intelectuales de los autores y además, hacer referencia a esta revista.

Recepción: 14 Diciembre 2017

Aprobación: 03 Mayo 2018

Resumen: En el presente trabajo de investigación se estudia el comportamiento numérico de perfiles estructurales de acero de pared delgada, de producción nacional venezolana Conduven ECO, de sección transversal cuadrada, construidos en forma de viga tipo cantiléver, bajo los efectos de Fatiga de Ultra-Bajo Ciclaje. Se realizaron ensayos experimentales donde se aplicaron historias de desplazamientos cíclicos, con muy baja velocidad y actuando en un rango que va desde el esfuerzo cedente hasta el esfuerzo último del elemento, obteniendo una degradación de la capacidad resistente de los perfiles, estos resultados experimentales fueron comparados con las simulaciones numéricas realizadas mediante un programa comercial de elementos finitos, obteniendo resultados satisfactorios.

Palabras clave: fatiga de ultra-bajo ciclaje, elementos tubulares de acero, falla por pandeo local.

Abstract: In the present work of investigation the numerical behavior of structural profiles of steel of thin wall is studied, of Venezuelan national production Conduven ECO, of square cross section, constructed in the form of beam type cantilever, under the effects of Fatigue of Ultra-Low Cycling. Experimental tests were carried out where histories of cyclic displacements were applied, with very low speed and acting in a range that goes from the yielding effort to the last effort of the element, obtaining a degradation of the resistant capacity of the profiles, these experimental results were compared with the numerical simulations made through a commercial program of finite elements, obtaining satisfactory results.

Keywords: ultra-low cycling fatigue, tubular steel elements, fails due to local buckling.

1. INTRODUCCIÓN

El fenómeno de Fatiga de Ultra-Bajo Ciclaje (FUBC) observado en ensayos experimentales es el causante de la degradación de la resistencia en los elementos cuando cargas cíclicas, inducidas mediante historias de desplazamientos controladas, aplicadas en un rango de esfuerzos que varían entre el esfuerzo cedente y el esfuerzo último sin superarlo y de manera repetitiva [1]. La aplicación y modelado de estos ensayos puede verse de manera extensa en [2]. Para simular el comportamiento observado en los ensayos planteados en la anterior referencia, se utilizó el modelo matemático elaborado por [3], en el que se llevaron a cabo experimentales en elementos estructurales de acero de sección tubular para validar un modelo matemático de elementos finitos para miembros de pórticos de acero. A este modelo se refieren todas las ecuaciones presentadas en este documento, mediante la implementación de una rutina escrita en el lenguaje de programación FORTRAN, que interactúa con el programa de elementos finitos no lineales ABAQUS, se logró ejecutar y sus resultados fueron satisfactorios. Para llevar a cabo dicha tarea se tomó en cuenta la notación mostrada en la Figura 1.

a) Deformaciones generalizadas, b) Esfuerzos generalizados.
Figura 1.
a) Deformaciones generalizadas, b) Esfuerzos generalizados.
Fuente: [3]

Se introduce la matriz generalizada de deformaciones plásticas , donde son las rotaciones plásticas de las rotulas inelásticas del modelo, el ultimo valor de la matriz es 0 porque no se considera elongación del elemento.

También es necesario describir que en el modelo se introduce la variable de daño inducido por momentos positivos ) y daño inducido por momentos negativos , que pueden tomar valores entre 0 y 1. Un elemento estructural que considere daño por momentos negativos y positivos en ambos nodos (i,j) queda representado por el esquema mostrado en la Figura 2.

Modelo de daño considerado.
Figura 2 .
Modelo de daño considerado.
Fuente: [3].

2. LEY DE ESTADO PARA PERFILES DE PARED DELGADA

La Ley de Estado bajo la cual esta formulado el modelo matemático usado en este trabajo de investigación relaciona las deformaciones generalizadas con los esfuerzos generalizados es la siguiente:

En las matrices (2) y (3), I representa la inercia de la sección transversal del elemento y E el módulo de elasticidad del material; en la ecuación (2) el valor de representa las rotación en el nodo y representa las rotaciones permanentes.

2.1. Función de fluencia para perfiles de pared delgada

Para considerar la existencia de dos pandeos diferentes en cada rotula (ver Figura 2), el modelo considera la evolución del daño en cada una de manera independiente mediante el modelo de endurecimiento cinemático e isotrópico no lineal. En cada rótula inelástica se emplea las ecuaciones (4) y (5), el símbolo Max( ) indica que la expresión toma en cuenta los pandeos positivos y negativos seleccionándose aquella que resulte mayor.

Función de endurecimiento cinemático e isotrópico no lineal
Figura 3 .
Función de endurecimiento cinemático e isotrópico no lineal
Fuente: [3]

El conjunto de ecuaciones permite representar el comportamiento del endurecimiento cinemático e isotrópico no lineal mostrado en la Figura 3, que al compararlo con otra ley de endurecimiento como lo es la representada de manera tri-lineal, resulta más ajustado a una curva de momento vs rotación de un ensayo experimental real.

2.2. Función de resistencia al pandeo local

Según la investigación [3], la evolución del daño ocasionado por pandeo local representada en la Figura 4, fue obtenida mediante ensayos experimentales, debiéndose introducir el concepto de rotación plástica conductora del pandeo Pcr, para identificar el valor máximo de rotación plástica a partir del cual la variable de daño se activa y evoluciona.

Comparación de funciones de resistencia al pandeo local.
Figura 4 .
Comparación de funciones de resistencia al pandeo local.
Fuente: [3]

Se considera la representación de la evolución del daño por pandeo local como una función de tendencia logarítmica mostrada en las ecuaciones (8), (9), (10) y (11), donde el parámetro b está relacionada con la velocidad del daño, mientras más alto sea el valor de b el pandeo local evoluciona de manera más rápida. La expresión usada en el modelo matemático que se ajusta al comportamiento descrito anteriormente, se denomina resistencia al pandeo local R(d), representa la rapidez del avance del daño con respecto a la rotación plástica. Las ecuaciones (8) y (9) expresan la resistencia al daño en el nodo i por momentos negativos y positivos, y las ecuaciones (10) y (11) la resistencia al daño en el nodo j.

El otro parámetro a tomar en cuenta es dult, que toma en cuenta la capacidad remanente de las probetas que se aprecia en niveles de daño muy grandes, una vez se inicia la etapa de deterioro la resistencia de las probetas disminuye a medida que avanza el daño, pero tiende en forma asintótica hacia un valor de resistencia remanente que es un porcentaje de la resistencia original de la pieza llamado dult, ilustrado en la Figura 5.

Esquema de zona donde actúa el dult
Figura 5 .
Esquema de zona donde actúa el dult
Fuente: [3]

2.3. Leyes de evolución de daño para perfiles de pared delgada

Según el modelo matemático empleado, la evolución del daño representado en la Figura 4, se activa luego que se acumulan rotaciones plásticas en el elemento hasta que alcanzan un valor mayor a Pcr, la acumulación de deformaciones plásticas se denomina . En caso de ensayos cíclicos, según se plantea en [4], las rotaciones plásticas de signo positivo o negativo tienden a causar pandeo local del mismo signo y a su vez, contribuyen a incrementar la resistencia al pandeo local de signo contrario, a este fenómeno se le denomina contra-pandeo.

En las simulaciones numéricas de los ensayos experimentales realizados con perfiles de acero de pared delgada, se agrega un parámetro llamado h, con el fin de cuantificar este fenómeno se introducen expresiones que relacionan las rotaciones conductoras del pandeo local. Las ecuaciones (12) y (13) representan las funciones de rotación conductora del pandeo local en el nodo i, y las ecuaciones (14) y (15) en el nodo j.

Las Leyes de Evolución del Daño para momentos positivos quedan definidas por la ecuación (16), establece que si las rotaciones conductoras de pandeo local son menores a la resistencia al pandeo local no hay incrementos del daño ( . El daño solo puede incrementarse si las rotaciones conductoras de pandeo local son iguales a la resistencia al pandeo local y el daño solo toma valores entre 0 y 1.

2.4. Determinación de parámetros de calibración del modelo matemático usado en las simulaciones numéricas

Para determinar y calibrar los parámetros del modelo matemático desarrollado por [3] se utilizaron perfiles cuadrados tubulares de pared delgada CONDUVEN ECO (producción venezolana), de sección cuadrada de 10cmx10cm de lado y espesor de pared 3,00 mm para los cuales se tomó en cuenta que todos los efectos inelásticos se producen en la rótula del nodo i mientras que la rótula del nodo j queda intacta (ver Figura 6). La nomenclatura empleada es:

Esquema de idealización de elementos a ensayar
Figura 6 .
Esquema de idealización de elementos a ensayar
Fuente: [3]

Los parámetros de calibración del modelo matemático que involucran las ecuaciones expresadas anteriormente, fueron determinados a partir de los ensayos experimentales descritos en [5] siendo:

E= módulo de elasticidad del material

A= área de la sección transversal del perfil

I= inercia de la sección transversal del perfil

β= parámetro que describe el endurecimiento isotrópico no lineal del material

Pcr= rotación plástica a partir de la cual evoluciona el daño

dult= parámetro que relaciona la tendencia del daño con las rotaciones plásticas

Mu= momento ultimo del perfil ensayado

My= momento cedente del perfil ensayado

b= valor que representa la velocidad del daño

h= factor de contra-pandeo en ensayos cíclicos

Con la realización de un primer ensayo [5], se logró obtener los límites para definir las historias de desplazamientos a ser impuestas en el resto de los ensayos experimentales, estos límites fueron Py=1250kg, δy=0,762cm, Pu=2088kg y δu=2,075cm (ver Figura 7). Con estos datos se ejecutaron dos ensayos de amplitudes ±1,2cm y ±1,6cm, cuyos resultados fueron comparados con los resultados numéricos modelados mediante el elemento finito descrito en la investigación [3]. El parámetro Pcr representa la rotación alcanzada en una curva de momento vs rotación, de un ensayo monotónico como se muestra en la Figura 8, a partir de la cual las pendientes de descarga en cada ciclo comienzan a disminuir debido al incremento y evolución del daño, lo que implica una degradación de la resistencia del perfil.

Curva de comportamiento (δ vs P) del ensayo inicial
Figura 7 .
Curva de comportamiento (δ vs P) del ensayo inicial
Fuente: los autores

Curva de comportamiento (∅ vs M) del ensayo inicial
Figura 8.
Curva de comportamiento (∅ vs M) del ensayo inicial
Fuente: los autores

Tomando en cuenta el modelo matemático elaborado por [3] aplicado en esta investigación para describir el fenómeno de FUBC, el parámetro b incluido en las ecuaciones (8) (9) (10) y (11), tiene un valor de 5,5 para el cual las simulaciones numéricas coinciden en gran medida con los ensayos experimentales. Para calcular el parámetro β, que está relacionado con el endurecimiento isotrópico no lineal en el elemento, se usan las ecuaciones (6) y (7) en puntos consecutivos de la curva característica de momento vs rotación en un ensayo monotónico a flexión, como se muestra en la Figura 9.

Curva de comportamiento de ensayo inicial donde se toman los datos para la ecuación (7)
Figura 9 .
Curva de comportamiento de ensayo inicial donde se toman los datos para la ecuación (7)
Fuente: los autores

Con los 4 puntos resaltados se realizó una tabla (ver Tabla 1) para aplicar la ecuación (7), donde:

Tabla 1.
Valores considerados para el cálculo de β.
Punto
1200180,001404137810,001001
2338890,002405212080,003499
3550970,00590450590,001724
4601560,007628--
Nota: My(kg*cm)=100.000,00 y Mu(kg*cm)=167.000,00
Fuente: los autores

El valor de la variable dult en las ecuaciones (8) (9) (10) y (11) toma en cuenta la capacidad remanente de las probetas, una vez que se inicia la etapa de deterioro la resistencia de las probetas disminuye a medida que avanza el daño, pero tiende en forma asintótica hacia un valor de resistencia remanente que es un porcentaje de la resistencia original de la pieza llamado dult, como se mostró en la Figura 5. El valor h que representa el fenómeno de contra-pandeo según [3] y que involucran las ecuaciones (12) (13) (14) y (15), fue estimado según este autor mediante la aplicación de ensayos experimentales cíclicos y su comparación con ensayos monosigno.

En el perfil estudiado la inercia es I=164,39cm4 y el área de la sección transversal A=10,58cm2, para la determinación del módulo de elasticidad experimental se empleó la zona elástica de la gráfica momento-rotación del ensayo de identificación de Py y Pu que va desde el punto (0;0) hasta el punto (0,00921;100.000), en donde a partir del valor de momento de 100.000 kgf*cm el comportamiento deja de ser netamente elástico. En esta zona se puede aplicar la relación entre momento rotación de una viga empotrada en un extremo y libre en el otro M=3EIϕ⁄L en donde se despeja el valor de EI=LM⁄3ϕ. Para el ensayo M=100.000 kgf*cm, L=80cm, ϕ=0,009952, por lo que EI=268.100.000kgf*cm2 y EA=174.000 kgf.

El valor de la variable dult en las ecuaciones (8) (9) (10) y (11) fue determinado por [3] y el valor de los resultados ajustados a los ensayos experimentales realizados es un valor promedio de 0,8. El valor que representa el fenómeno de contra-pandeo según [3] y que involucran las ecuaciones (12) (13) (14) y (15), es h=-0,8, estimado según este autor mediante la aplicación de ensayos experimentales cíclicos y su comparación con ensayos monosigno. Los parámetros de calibración introducidos en el modelo matemático que se usó para describir el fenómeno de FUBC en el presente trabajo de investigación fueron: EI=268.100.000 kgf*cm², EA=174.000 kgf, β= 180, Pcr=0,01150, dult =0,8, Mu= 167.000kgf*cm, My= 100.000 kgf-cm, b=5,50, h=-0,80.

3. METODOLOGÍA

3.1. Implementación del modelo matemático usado para representar el fenómeno de FUBC

Para realizar la simulación numérica de los ensayos experimentales antes descritos, se aplicó el modelo matemático propuesto y validado en la investigación [3] mediante un programa de elementos finitos desarrollado por [6], el cual puede interactuar con el software de cálculo ABAQUS [7]. Este programa permite al usuario crear un elemento finito e introducirlo al programa por medio de subrutinas de programación (UEL, User’s Elements) [6]; las características y comportamiento del elemento finito son introducidos en el software principal mediante el lenguaje de programación FORTRAN, permitiendo procesar la entrada de datos para así calcular deformaciones, fuerzas internas, esfuerzos, variables internas y todos los elementos necesarios para la solución del problema (ver Figura 10).

La rutina escrita en lenguaje de programación FORTRAN realiza la resolución del problema paso a paso, mediante incrementos a través del tiempo; esta solución se divide en dos problemas: local y global (ver Figura 11).

Ventana de trabajo para rutinas escritas en lenguaje FORTRAN
Figura 10.
Ventana de trabajo para rutinas escritas en lenguaje FORTRAN
Fuente: los autores

El problema local se resuelve mediante el uso de la rutina “mod90a.f” basada en la teoría de pórticos y la Teoría de Daño Concentrado [8], esta rutina resuelve el problema central del análisis no lineal de pórticos degradables. En cambio, el problema global por la rutina “AnalisePorticosFadiga.f90”, encargada de tomar los valores del resultado de la rutina definida por el usuario e incorporarlos a su programación.

Esquema de funcionamiento general de interacción de rutinas escritas en FORTRAN
Figura 11.
Esquema de funcionamiento general de interacción de rutinas escritas en FORTRAN
Fuente: los autores

4. ANÁLISIS Y DISCUSIÓN DE RESULTADOS

4.1. Comparación de ensayos experimentales de diferentes amplitudes con simulaciones numéricas

Todos los ensayos descritos están referidos a la investigación [5], donde se trata a detalle los procedimientos usados para realizarlos. Para la simulación del ensayo de amplitud ±1,6cm se tomó en cuenta el modelo de FUBC aplicando el modelo matemático mencionado anteriormente [3], y se puede observar en la Figura 13 que la máxima rotación plástica que puede ocurrir en este ensayo es 0,004622 y un momento máximo de 155.000 kgf*cm, sin embargo mediante la acumulación de rotaciones plásticas sobrepasa el valor de Pcr=0,01150, y se activa el daño en el perfil. Los parámetros usados para la simulación de este ensayo fueron los visualizados en la Figura 12.

Comparación de las curvas características del ensayo de amplitud ±1,6 cm: a) Ensayo experimental (b) Simulación numérica
Figura 12 .
Comparación de las curvas características del ensayo de amplitud ±1,6 cm: a) Ensayo experimental (b) Simulación numérica
Fuente: los autores

Mediante la comparación de la curva experimental con la simulación numérica, se puede notar que el modelo matemático para FUBC es aplicable y describe con notable precisión el ensayo experimental realizado; la diferencia entre el valor máximo de momento del ensayo realizado y la simulación numérica fue de 1,55%. En la Figura 13 se observa la comparación de la fuerza generada en el ensayo y la simulación numérica.

Comparación de curva de fuerza vs tiempo del ensayo de amplitud ±1,6cm y la simulación numérica
Figura 13 .
Comparación de curva de fuerza vs tiempo del ensayo de amplitud ±1,6cm y la simulación numérica
Fuente: los autores

Se puede notar que a medida que pasa el tiempo la fuerza comienza a disminuir hasta la culminación del ensayo; esto se debe a que la sección transversal del perfil comienza a experimentar mayor degradación de su capacidad resistente a medida que pasa el tiempo, la diferencia entre estos valores numéricos y experimentales varía entre 0,78% y 1,40%.

En las Figuras 14 y 15 se tiene la evolución del daño producido tanto por momentos positivos como por momentos negativos en la simulación numérica del ensayo, donde se refleja que a medida que aumenta el número de ciclos, el daño aumenta y sigue esa tendencia hasta el final de la simulación.

Evolución del daño por momentos positivos. Ensayo de amplitud ±1,6 cm
Figura 14 .
Evolución del daño por momentos positivos. Ensayo de amplitud ±1,6 cm
Fuente: los autores

Evolución del daño por momentos negativos. Ensayo de amplitud ±1,6 cm
Figura 15 .
Evolución del daño por momentos negativos. Ensayo de amplitud ±1,6 cm
Fuente: los autores

En la simulación del ensayo de amplitud ±1,2 cm se observar en la Figura 16 que la rotación plástica máxima que puede ocurrir en este ensayo es 0,002022 y un momento máximo de 129.000 kgf*cm, sin embargo mediante la acumulación de rotaciones plásticas en la función sobrepasa el valor de Pcr =0,01150, y se activa el daño en la simulación numérica. Los parámetros introducidos para obtener la Figura 16 fueron los mismos que para el ensayo de amplitud ±1,2 cm, esto se debe a que se trata del mismo perfil con la misma sección transversal, siendo la única variante entre ensayos es la historia de desplazamientos. La diferencia entre el valor máximo de momento del ensayo realizado y la simulación numérica fue de 0,88%.

Comparación de las curvas característica del ensayo de amplitud ± 1,6cm: (a) Ensayo experimental (b) Simulación numérica
Figura 16 .
Comparación de las curvas característica del ensayo de amplitud ± 1,6cm: (a) Ensayo experimental (b) Simulación numérica
Fuente: los autores

En la Figura 17 se tiene la comparación de la fuerza generada en el ensayo 7 y su correspondiente simulación numérica, en ambas gráficas se nota que a medida que pasa el tiempo la fuerza comienza a disminuir hasta la culminación del ensayo, esto se debe a que la sección transversal del perfil comienza a experimentar mayores daños a medida que transcurre el tiempo, la diferencia entre los valores de fuerza está entre 2,68% y 3,20%.

Comparación de curva de fuerza vs tiempo del ensayo de ±1,2cm (cíclico) y la simulación numérica
Figura 17 .
Comparación de curva de fuerza vs tiempo del ensayo de ±1,2cm (cíclico) y la simulación numérica
Fuente: los autores

En las Figuras 18 y 19 se visualiza la evolución del daño producido tanto por momentos positivos como por momentos negativos, en la simulación numérica del ensayo cíclico de amplitud ±1,2 cm, reflejándose que a medida que aumenta el número de ciclos, el daño se acrecienta y sigue esa tendencia hasta el final de la simulación.

Daño por momentos positivos.  Ensayo de amplitud ±1,2 cm
Figura 18.
Daño por momentos positivos. Ensayo de amplitud ±1,2 cm
Fuente: los autores

Daño por momentos negativos.  Ensayo de amplitud ±1,2 cm
Figura 19.
Daño por momentos negativos. Ensayo de amplitud ±1,2 cm
Fuente: los autores

5. CONCLUSIONES

Aunque en la concepción de esta investigación se invirtió tiempo y esfuerzo para tratar de observar y de alguna manera medir la aparición y consecutiva propagación de grietas, el fenómeno que conduce el fallo de los perfiles cerrados de pared delgada es el fenómeno del pandeo local, y la degradación de la resistencia de los elementos ensayados ocurre cuando las rotaciones plásticas acumuladas superan el valor de Pcr.

Es importante resaltar por otro lado, que el empleo del software de análisis estructural para la simulación numérica puede ser una herramienta muy útil, si pueden introducirse los parámetros que representan condiciones de contorno de los problemas a estudiar. En esta investigación fue necesario la obtención de estos parámetros, como lo son b, β, Mu, My, h, dultI, A y E, que solo fue posible mediante la realización de ensayos experimentales, con los cuales se pudo calibrar el modelo empleado hasta obtener resultados numéricos muy acordes con lo observado en ensayos experimentales.

Por otro lado, el modelo matemático elaborado por [3] es válido para representar el fenómeno de FUBC; este modelo introduce la acumulación de rotaciones plásticas como elemento para activar el daño en el perfil. Cuando la rotación plástica acumulada es mayor al valor de Pcrentonces se activa la evolución del daño como se pudo observar en las gráficas comparativas entre los ensayos experimentales realizados y las simulaciones numéricas.

5. REFERENCIAS

[1]A. Hamidreza, A. Aghakouchak, S.Sharif y M. D. Engelhardt «Finite element simulation of ultra low cycle fatigue cracking in steel structures» Journal of Constructional Steel Research, 89, pp 175-184, 2013

[2]L. Lopez, R.Picon y M. Marante «ANÁLISIS EXPERIMENTAL DE PERFILES TUBULARES DE PARED DELGADA SOMETIDOS A FATIGA DE ULTRA-BAJO CICLAJE» Revista Gaceta Técnica, 18, pp 23-34, 2017

[3]R. Febres, "Modelo de daño para pórticos planos de acero bajo cargas histeréticas" Merida, Venezuerla, 2002

[4]R. Febres, P. Inghlesis y J. Florez-Lopez «Modeling of local buckling in tubular steel frames subjected to cyclic loading» Computers and Structures, 81, pp 32-46, 2003

[5]L.Lopez, R. Picon y M.Marante «ANALISIS EXPERIMENTAL DE PERFILES TUBULARES DE PARED DELGADA SOMETIDOS A FATIGA DE ULTRA-BAJO CICLAJE» Trabajo de Grado, Universidad Centroccidental Lisandro Alvarado, Barquisimeto, Venezuela, 2017

[6]M. Y. Uzcategui «Desarrollo de un programa de elementos finitos tridimensional basado en la WEB» Merida. 2012.

[7]HKS «Abaqus User's Manual 6.2-1» Hibbit, Karlson and Soresen Inc, Providence, EEUU, 2000

[8]P. Inglessis «Modelado del comportamiento de pórticos de acero mediante la teoria del daño concentrado» Merida, Venezuela, 1999

Artículo relacionado

[Artículo corregido , vol. 19 ] https://revistas.ucla.edu.ve/index.php/gt/article/view/1146

HTML generado a partir de XML-JATS4R por