Análisis Cluster para Ensayo de Tensión Indirecta
Juan David Oliveros
2023-08-28
##Introducción ¿Qué es el ensayo de tensión indirecta?
El ensayo de tensión indirecta, también llamado “Ensayo brasilero”, es un ensayo de caracterización del concreto, en el que se utiliza un cilindro de concreto con radio R y altura 2R, y se comprime en posición horizontal, como se muestra en el siguiente esquema:
Esquema de un ensayo de tensión indirecta.
Fuente: Elaboración Propia
Este ensayo se denomina de tensión indirecta, pues aunque el ensayo consiste en una compresión del elemento, el modo de falla que presenta es una tensión, donde las partículas del cilindro intentarán desprenderse una de las otras, generando la falla. Un esquema del modo de falla se muestra a continuación:
Modo de falla en un ensayo de tensión indirecta.
Fuente: Elaboración Propia
Los datos que se capturan de este ensayo son desplazamiento del elemento superior que comprime el cilindro, en milímetros, y la fuerza requerida para generar ese desplazamiento, en kilonewtons.
De esta forma, se analizarán 3 ensayos de cilindros de concreto reforzados con fibras metálicas, para conocer su comportamiento, en el rango plástico.
Previo al presente análisis, se realizó un tratamiento de datos, separando el rango elástico de los ensayos realizados, pues al presentar la primera fisura, se presenta un salto en la curva de desplazamiento contra fuerza, salto que no es de relevancia para el presente análisis.
library(ggplot2)
library(growthmodels)
library(Metrics)
library(minpack.lm)#Definición del vector de desplazamientos, en mm
x=c(0,0,0,0.067,0.067,0.067,0.133,0.133,0.133,0.2,0.2,0.2,0.267,0.267,0.267,0.333,0.333,0.333,0.4,0.4,0.4,0.467,0.467,0.467,0.533,0.533,0.533,0.6,0.6,0.6,0.667,0.667,0.667,0.733,0.733,0.733,0.8,0.8,0.8,0.867,0.867,0.867,0.933,0.933,0.933,1,1,1,1.067,1.067,1.067,1.133,1.133,1.133,1.2,1.2,1.2,1.267,1.267,1.267,1.333,1.333,1.333,1.4,1.4,1.4,1.467,1.467,1.467,1.533,1.533,1.533,1.6,1.6,1.6,1.667,1.667,1.667,1.733,1.733,1.733)
#Definición del vector de carga, en kN
load=c(0,0,0,3594.2,4011.12,4862.52,7666.83,9675.42,9913.5,11840.1,15111.82,14921.01,15958.48,20478.58,19416.74,19750.24,25541.42,23610.51,22989.5,30155.95,27849.22,25899.64,34272.91,31632.74,28270.81,37674.92,35288.43,30393.99,41042.01,38604.61,32209.43,44379.57,41611.57,34015.35,47250.66,44547.5,35298.5,49742.99,47144.8,36499.05,51750.7,49478.8,37381.55,53508.46,51436.38,38685.21,55155.95,52844.41,40056.55,56862.2,53553.77,41074.43,58463.49,54327.82,42323.73,59883.09,54978.27,43561.56,61146.62,55491.21,45077.01,61815.37,55881.05,46720.77,62407.8,56404.4,48258.59,62763.21,56782.18,49404.3,62660.11,57330.76,50523.42,62466.46,57826.36,51627.64,62703.57,58387.79,52448.54,62812.75,59000.22)
plot(load~x,pch=16)
#Se redefinen los vectores porque el 1er ensayo consultado está generando distorsiones en el análisis. Se recomienda realizar más de 3 ensayos para tener una muestra representativa.
#Re-definición del vector de desplazamientos, en mm
x=c(0,0,0.067,0.067,0.133,0.133,0.2,0.2,0.267,0.267,0.333,0.333,0.4,0.4,0.467,0.467,0.533,0.533,0.6,0.6,0.667,0.667,0.733,0.733,0.8,0.8,0.867,0.867,0.933,0.933,1,1,1.067,1.067,1.133,1.133,1.2,1.2,1.267,1.267,1.333,1.333,1.4,1.4,1.467,1.467,1.533,1.533,1.6,1.6,1.667,1.667,1.733,1.733)
#Re-definición del vector de carga, en kN
load=c(0,0,4011.12,4862.52,9675.42,9913.5,15111.82,14921.01,20478.58,19416.74,25541.42,23610.51,30155.95,27849.22,34272.91,31632.74,37674.92,35288.43,41042.01,38604.61,44379.57,41611.57,47250.66,44547.5,49742.99,47144.8,51750.7,49478.8,53508.46,51436.38,55155.95,52844.41,56862.2,53553.77,58463.49,54327.82,59883.09,54978.27,61146.62,55491.21,61815.37,55881.05,62407.8,56404.4,62763.21,56782.18,62660.11,57330.76,62466.46,57826.36,62703.57,58387.79,62812.75,59000.22)
plot(load~x,pch=16)
#Se calculan las medias de los ensayos realizados:
medias = tapply(load,x,mean)
plot(load~x,pch=16)
points(medias~unique(x),pch=16,col="red",cex=1)
abline(v=x,col="grey90",lty=2)
medias_diff= diff(medias)
#Se divide entre 0.0666 pues todas las mediciones están hechas para cada 0.0666 mm
TAC=medias_diff/0.0666
TAC## 0.067 0.133 0.2 0.267 0.333 0.4 0.467 0.533
## 66618.919 80445.045 78407.733 74042.718 69494.069 66465.766 59312.913 52985.736
## 0.6 0.667 0.733 0.8 0.867 0.933 1 1.067
## 50174.700 47631.532 43596.246 38210.435 32595.420 27892.943 22939.339 18135.210
## 1.133 1.2 1.267 1.333 1.4 1.467 1.533 1.6
## 17832.883 15540.916 13336.862 7947.372 8376.727 5504.429 3344.444 2266.892
## 1.667 1.733
## 5995.045 5417.492#Se usa la librería nlsLM() que, al parecer, es una opción más robusta que nls. Ambas funciones (nls() y nlsLM()) llegan a la misma solución, pero hay que "iterar" más con nls() para llegar a la solución.
#Ajuste del modelo de crecimiento de Brody
model_1=nlsLM(load~a*(1-exp(-k*x)),start = list(a=1000,k=1))
summary(model_1)##
## Formula: load ~ a * (1 - exp(-k * x))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 6.828e+04 1.268e+03 53.85 <2e-16 ***
## k 1.458e+00 6.360e-02 22.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2278 on 52 degrees of freedom
##
## Number of iterations to convergence: 14
## Achieved convergence tolerance: 1.49e-08#Ajuste del modelo de crecimiento de Weibull
model_2=nlsLM(load~a-(b*exp(-k*x)),start = list(a=1000,b=1000,k=1))
summary(model_2)##
## Formula: load ~ a - (b * exp(-k * x))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 6.684e+04 1.180e+03 56.66 <2e-16 ***
## b 6.956e+04 1.156e+03 60.15 <2e-16 ***
## k 1.593e+00 8.009e-02 19.89 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2148 on 51 degrees of freedom
##
## Number of iterations to convergence: 13
## Achieved convergence tolerance: 1.49e-08#Ajuste del modelo de crecimiento loglogistic
model_3=nlsLM(load~a/(1+exp(-k*log(x))),start = list(a=1000,k=1))
summary(model_3)##
## Formula: load ~ a/(1 + exp(-k * log(x)))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 1.028e+05 8.409e+02 122.31 <2e-16 ***
## k 9.877e-01 2.995e-02 32.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2801 on 52 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.49e-08#Con los parámetros encontrados anteriormente, se realiza la gráfica de los modelos utilizados
plot(load~x,pch=16)
curve((coef(model_1)[1]*(1-exp(-coef(model_1)[2]*x))),add=T,col="orange", lwd=2)
curve((coef(model_2)[1]-(coef(model_2)[2]*exp(-coef(model_2)[3]*x))),add=T,col="palevioletred3", lwd=2)
curve((coef(model_3)[1]/(1+exp(-coef(model_3)[2]*log(x)))),add=T,col="aquamarine3", lwd=2)
legend(x=1, y=20000, legend=c("Modelo Brody", "Modelo Weibull", "Modelo loglogistic"), col=c("orange", "palevioletred3","aquamarine3"), lty=1, cex=0.8, lwd=2)
#Haciendo predicciones de cada modelo:
pred_1=fitted.values(model_1)
pred_2=fitted.values(model_2)
pred_3=fitted.values(model_3)
#RSME compara la raíz de los menores cuadrados
Metrics::rmse(actual=load,predicted=pred_1)## [1] 2235.262Metrics::rmse(actual=load,predicted=pred_2)## [1] 2087.852Metrics::rmse(actual=load,predicted=pred_3)## [1] 2748.542#El segundo modelo genera un error menor. Se toma el modelo de Weibull.#Se define la ecuación de crecimiento adoptada y su derivada
f=expression(a-(b*exp(-k*x)))
D(f,"x")## b * (exp(-k * x) * k)#Se define el TAC_i, en función de la derivada calculada anteriormente
TAC_i=function(x){
y=(coef(model_2)[2]*exp(-coef(model_2)[3]*x)*coef(model_2)[3])
return(y)
}#Se definen 1000 puntos entre 0.5mm y 1.2mm, que es donde se evidencia que hay un comportamiento con cambios de pendientes importantes.
#No se define desde 0mm porque, de 0mm a 0.5mm, la curva se comporta aproximadamente lineal.
x=seq(0.5,1.2,length.out=1000)
y=TAC_i(x)
hist(y);abline(v=mean(y),col="red",lwd=2);abline(v=quantile(y,0.25),col="blue",lwd=2)
legend(x=37500, y=200, legend=c("Media", "Q25"), col=c("red", "blue"), lty=1, cex=0.8, lwd=2)
#Análisis de cluster jerárquico
y_mat = dist(y,method="euclidean")
Hierar_cl=hclust(y_mat, method="ward.D2")
#Crea las etiquetas, en qué datos están en qué cluster
fit = cutree(Hierar_cl,k=5)
#Se ingresan los datos del cluster en un dataframe
df_c=data.frame(y, cluster=fit)
#Se grafica un diagrama de cajas, con los 5 clusters definidos, indicando cómo varía la media de cada cluster y la media del rango
boxplot(df_c$y~df_c$cluster)
points(x=1:5,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)
legend(x=4, y=45000, legend=c("Media del cluster", "Media del rango"), col=c("red", "blue"), lty=c(NA,2), pch=c(16,NA), cex=0.8, lwd=2)
#Se calcula la media del cluster
mean(df_c$y)## [1] 30117.82#Se grafican las derivadas del rango definido, mostrando cada cluster con un color diferente
plot(x=x,y=y,col=as.factor(df_c$cluster))
#Se calcula la media ponderada
loads=table(df_c$cluster)/1000
medias=tapply(df_c$y,df_c$cluster,mean)
mp=t(loads)%*%medias
mp## [,1]
## [1,] 30117.82boxplot(df_c$y~df_c$cluster)
#Se grafica de nuevo comparando que la media ponderada del rango es igual a la media del rango:
points(x=1:5,y=tapply(df_c$y,df_c$cluster,mean),pch=16,col="red")
abline(h=mean(df_c$y),lwd=2,col="blue",lty=2)
abline(h=mp,lwd=2,col="orange",lty=2)
legend(x=4, y=45000, legend=c("Media del cluster", "Media del rango", "Media ponderada"), col=c("red", "blue", "orange"), lty=c(NA,2,2), pch=c(16,NA,NA), cex=0.8, lwd=2)
De esta forma se concluye que la media en este rango, calculada por el
modelo de Weibull, puede dar información importante, al ser la misma que
la media ponderada. En este rango, se define entonces que la media de
resistencia es aproximadamente 30kN.