Nothing
#' Dynamic material modeling from strain rate temperature table
#'
#' @description Dynamic material modeling based on strain rate-temperature table returned from
#' the function \code{\link[TPMplt:epsExtract]{epsExtract}}. Material constants as well as power
#' dissipation efficiency factors and rheological stability coefficients in current conditions
#' will be returned.
#' @param x A strain rate-temperature table, returned from \code{\link[TPMplt:epsExtract]{epsExtract}}.
#' @param lgbase A numeric value to specify the base of the logarithm calculations for processing map.
#' The default value uses 10.
#' @param InteractMode A boolean value to control figures' output and the printout of related constants
#' during calculations. Default value FALSE means all fitting plots will not be outputed. If these outputs
#' are necessary, set this parameter as TRUE then follow the prompt messages.
#' @param ConsFunc A boolean value to determine whether calculating for constructive equation. The
#' default value uses FALSE.
#' @param legendcex A numeric value to determine the legend scale. It is activated only when the parameter
#' InteractMode is TRUE. The default value is 0.65.
#' @param legendloc A character object to determine the location of legend. It is activated only when the
#' parameter InteractMode is TRUE. The defualt value is "bottomright".
#'
#' @return Serial material constants, constructive function, eta table and xi table through
#' dynamic material model developed by Gegel and Prasad.
#' @import VBTree stats utils RColorBrewer
#' @export DMMprocess
#' @seealso \code{\link[VBTree:VBTree-package]{VBTree}}, \code{\link[TPMplt:epsExtract]{epsExtract}}
#'
#' @examples
#' require(VBTree)
#' dl2vbt(chrvec2dl(colnames(TPMdata)))
#' epstable <- epsExtract(TPMdata, 0.7, 2, 3)
#'
#' # Without calculation for constitutive equation
#' DMM <- DMMprocess(epstable)
#' message(DMM)
#'
#' # Calculating for constitutive equation but
#' # Without plots printout.
#' DMM <- DMMprocess(epstable, ConsFunc=TRUE)
#' message(DMM)
#'
#' # Calculating for constitutive equation and
#' # required fitting plots printout. (message and selection in prompt)
#' \dontrun{
#' DMMprocess(epstable, InteractMode=TRUE, ConsFunc=TRUE)
#' }
#'
#' @keywords DMMprocess DMMresult epsExtract
DMMprocess <- function(x, lgbase=10, InteractMode=FALSE, ConsFunc=FALSE, legendcex=0.65, legendloc="bottomright"){
# check the input data
if(!inherits(x, "SR-T.table")){
stop("the input data should be SR-T Chart generated by the function epsinprt()", call. = FALSE)
}
# # test section
# base <- exp(1)
# x <- epstable
eps <- x$epsilon
x <- x$SRT.table
if(ConsFunc==TRUE){
logbase <- exp(1)
} else{
logbase <- lgbase
}
dims <- dim(x)
SR <- as.numeric(rownames(x))
Tem1 <- rev(1/(as.numeric(colnames(x)) + 273.15)) # Kelvin
lgSR <- log(SR, base = logbase)
clrvec <- brewer.pal(9, "Set1")
if(ConsFunc == TRUE){
# Constitutive equation and related plotss
i <- 1
scatter_mat <- matrix(NA, nrow = dims[2], ncol = dims[1])
scatter_mat1 <- scatter_mat
lmlist <- list(list())[rep(1,dims[2])]
lmlist1 <- lmlist
for (i in 1:dims[2]) {
scatter_mat[i,] <- as.numeric(log(x[,i], base = logbase))
scatter_mat1[i,] <- as.numeric(x[,i])
a <- lgSR
lmlist[[i]] <- lm(scatter_mat[i,]~I(a^1))
lmlist1[[i]] <- lm(scatter_mat1[i,]~I(a^1))
}
temp_vec <- c()
temp_vec1 <- temp_vec
for (i in 1:dims[2]) {
a <- as.numeric(lmlist[[i]][[1]][2]) # Extract slope
a1 <- as.numeric(lmlist1[[i]][[1]][2])
temp_vec <- append(temp_vec, a)
temp_vec1 <- append(temp_vec1, a1)
}
n.StressInd <- mean(1/temp_vec)
beta.StressInd <- mean(1/temp_vec1)
alpha.StressInd <- beta.StressInd/n.StressInd # unit: MPa^(-1)
#Plot1 output section
if(InteractMode == TRUE){
message("Show fitting plot of log_flow_stress vs. log_strain_rate?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes"){
legendvec <- colnames(x)
yscale <- c(min(unlist(scatter_mat)), max(unlist(scatter_mat)))
xscale <- c(min(lgSR), max(lgSR))
labx <- expression(paste("ln(",ring(epsilon), "/s"^-1, ")", sep = ""))
laby <- expression(paste("ln(",sigma, "/MPa)", sep = ""))
for(i in 1:dims[2]){
plot(lgSR, scatter_mat[i,], type = "p", col = clrvec[i], pch=16, xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
lines(x=lgSR, y=predict(lmlist[[i]]), col = clrvec[i], xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
if (i == dims[2]){
legend(legendloc,legend = legendvec, fill = clrvec[1:dims[2]], cex = legendcex, bg = "transparent", box.lty = 0)
cus_par <- par(new=FALSE)
on.exit(par(cus_par))
} else {
cus_par <- par(new=TRUE)
on.exit(par(cus_par))
}
}
}
message("Mean value of n.Stress.Index calculated from log_flow_stress vs. log_strain_rate fitting is", n.StressInd,"\n")
invisible(readline(prompt="Press [enter] to continue"))
}
#Plot2 output section
if(InteractMode == TRUE){
message("Show fitting plot of flow_stress vs. log_strain_rate?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes"){
legendvec <- colnames(x)
yscale <- c(min(unlist(scatter_mat1)), max(unlist(scatter_mat1)))
xscale <- c(min(lgSR), max(lgSR))
labx <- expression(paste("ln(",ring(epsilon), "/s"^-1, ")", sep = ""))
laby <- expression(paste(sigma, "/MPa", sep = ""))
for(i in 1:dims[2]){
plot(lgSR, scatter_mat1[i,], type = "p", col = clrvec[i], pch=16, xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
lines(x=lgSR, y=predict(lmlist1[[i]]), col = clrvec[i], xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
if (i == dims[2]){
legend(legendloc,legend = legendvec, fill = clrvec[1:dims[2]], cex = legendcex, bg = "transparent", box.lty = 0)
cus_par <- par(new=FALSE)
on.exit(par(cus_par))
} else {
cus_par <- par(new=TRUE)
on.exit(par(cus_par))
}
}
message("Mean value of beta.Stress.Index calculated from flow_stress vs. log_strain_rate fitting is", beta.StressInd,"\n")
invisible(readline(prompt="Press [enter] to continue"))
}
}
sinhax <- sinh(alpha.StressInd*x)
lgsinhax <- log(sinhax, base = logbase)
scatter_mat <- matrix(NA, nrow = dims[2], ncol = dims[1])
lmlist <- list(list())[rep(1,dims[2])]
for (i in 1:dims[2]) {
scatter_mat[i,] <- as.numeric(lgsinhax[,i])
a <- lgSR
lmlist[[i]] <- lm(scatter_mat[i,]~I(a^1))
}
temp_vec <- c()
for (i in 1:dims[2]) {
a <- as.numeric(lmlist[[i]][[1]][2])
temp_vec <- append(temp_vec, a)
}
N.StressInd <- mean(1/temp_vec)
#Plot3 output section
if(InteractMode == TRUE){
message("Show fitting plot of log[sinh(alpha*flow_stress)] vs. log_strain_rate?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes"){
legendvec <- colnames(x)
yscale <- c(min(unlist(scatter_mat)), max(unlist(scatter_mat)))
xscale <- c(min(lgSR), max(lgSR))
labx <- expression(paste("ln(",ring(epsilon), "/s"^-1, ")", sep = ""))
laby <- expression(paste("ln[sinh(", alpha, sigma, ")]", sep = ""))
for(i in 1:dims[2]){
plot(lgSR, scatter_mat[i,], type = "p", col = clrvec[i], pch=16, xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
lines(x=lgSR, y=predict(lmlist[[i]]), col = clrvec[i], xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
if (i == dims[2]){
legend(legendloc,legend = legendvec, fill = clrvec[1:dims[2]], cex = legendcex, bg = "transparent", box.lty = 0)
cus_par <- par(new=FALSE)
on.exit(par(cus_par))
} else {
cus_par <- par(new=TRUE)
on.exit(par(cus_par))
}
}
message("Mean value of N.Stress.Index calculated from log[sinh(alpha*flow_stress)] vs. log_strain_rate is", N.StressInd,"\n")
invisible(readline(prompt="Press [enter] to continue"))
}
}
scatter_mat <- matrix(NA, nrow = dims[1], ncol = dims[2])
T_1 <- 1/(as.numeric(colnames(x)) + 273.15) # Convert to Kelvin
lmlist <- list(list())[rep(1,dims[1])]
for (i in 1:dims[1]) {
scatter_mat[i,] <- as.numeric(lgsinhax[i,])
a <- T_1
lmlist[[i]] <- lm(scatter_mat[i,]~I(a^1))
}
temp_vec <- c()
for (i in 1:dims[1]) {
a <- as.numeric(lmlist[[i]][[1]][2])
temp_vec <- append(temp_vec, a)
}
K <- mean(temp_vec)
Q <- K * N.StressInd * 8.314 # R Unit: J/(mol*K)
#Plot4 output section
if(InteractMode == TRUE){
message("Show fitting plot of log[sinh(alpha*flow_stress)] vs. 1/T(K^-1)?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes"){
legendvec <- rownames(x)
yscale <- c(min(unlist(scatter_mat)), max(unlist(scatter_mat)))
xscale <- c(min(T_1), max(T_1))
labx <- "1/T(K^-1)"
laby <- expression(paste("ln[sinh(", alpha, sigma, ")]", sep = ""))
for(i in 1:dims[1]){
plot(T_1, scatter_mat[i,], type = "p", col = clrvec[i], pch=16, xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
lines(x=T_1, y=predict(lmlist[[i]]), col = clrvec[i], xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
if (i == dims[1]){
legend(legendloc,legend = legendvec, fill = clrvec[1:dims[2]], cex = legendcex, bg = "transparent", box.lty = 0)
cus_par <- par(new=FALSE)
on.exit(par(cus_par))
} else {
cus_par <- par(new=TRUE)
on.exit(par(cus_par))
}
}
message("Mean value of K calculated from log[sinh(alpha*flow_stress)] vs. 1/T(K^-1) is", K,"\n")
message("Activation Energy Q calculated based on K and N.Stress.Index is", Q,"J/mol\n")
invisible(readline(prompt="Press [enter] to continue"))
}
}
i <- 1
j <- 1
Zmat <- matrix(NA, nrow = dims[1], ncol = dims[2])
for (i in 1:dims[1]) {
for (j in 1:dims[2]) {
Tcond <- as.numeric(colnames(lgsinhax)[j]) + 273.15
SRcond <- as.numeric(rownames(lgsinhax)[i])
temp <- exp(Q/(8.314*Tcond)) # R=8.314
Zmat[i,j] <- SRcond*temp
}
}
Z <- as.numeric(as.vector(Zmat))
lnZ <- log(Z, base = logbase)
a1 <- as.numeric(as.vector(lgsinhax))
lm1 <- lm(lnZ~I(a1^1))
n.Power <- as.numeric(lm1[[1]][2])
#Plot5 output section
if (InteractMode == TRUE) {
message("Show fitting plot of log_Z vs. log[sinh(alpha*flow_stress)]?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes") {
labx <- expression(paste("ln[sinh(", alpha, sigma, ")]", sep = ""))
laby <- "lnZ"
plot(a1, lnZ, type = "p", col = "black", pch=16, xlab = labx, ylab = laby)
lines(x=a1, y=predict(lm1), col = "red", xlab = labx, ylab = laby)
message("Value of n.Index calculated from log_Z vs. log[sinh(alpha*flow_stress)] is", n.Power,"\n")
invisible(readline(prompt="Press [enter] to continue"))
}
}
a2 <- as.numeric(as.vector(sinhax^n.Power))
lm2 <- lm(Z~I(a2^1))
A <- as.numeric(lm2[[1]][2])
#Plot6 output section
if (InteractMode == TRUE) {
message("Show fitting plot of Z vs. [sinh(alpha*flow_stress)]^n?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes") {
labx <- expression(paste("[sinh(", alpha, sigma, ")]^n", sep = ""))
laby <- "Z"
plot(a2, Z, type = "p", col = "black", pch=16, xlab = labx, ylab = laby)
lines(x=a2, y=predict(lm2), col = "red", xlab = labx, ylab = laby)
message("Value of A.Constant calculated from Z vs. [sinh(alpha*flow_stress)]^n is", A,"\n")
invisible(readline(prompt="Press [enter] to continue"))
}
}
Q.ActivEnerg <- Q*0.001 # Unit: kJ/mol
if(ConsFunc){
message("Constitutive equation in", eps, "strain:\nEpsilon=A{[sinh(Alpha*Sigma)]^n}*exp[-Q/(RT)]\nWhere A =",
A, "s^-1", "\n Alpha =", alpha.StressInd,
"MPa^-1 = ", format(alpha.StressInd, scientific = TRUE), "Pa^-1\n Q =", Q, "J/mol = ", Q.ActivEnerg,
"kJ/mol\n R = ", 8.314,"J/(mol*K)\n Epsilon is strain rate\n Sigma is flow stress\n T is Temperature (K)", "\n\n")
}
paras <- list(A = A, alpha = alpha.StressInd, n = n.Power, Q = Q.ActivEnerg, base = logbase, epsilon.strain = eps, Tcorr.n = n.StressInd) #, Tcorr.Mode = Tcorr.Mode)
}
# recover logarithm base
logbase <- lgbase
scatter_mat <- matrix(NA, nrow = dims[2], ncol = dims[1])
lmlist <- list(list())[rep(1,dims[2])]
for (i in 1:dims[2]) {
scatter_mat[i,] <- as.numeric(log(x[,i], base = logbase))
a <- lgSR
lmlist[[i]] <- lm(scatter_mat[i,]~I(a^1)+I(a^2)+I(a^3))
}
Coef_mat <- matrix(NA, nrow = dims[2], ncol = 4)
for (i in 1:dims[2]) {
Coef_mat[i,] <- as.numeric(lmlist[[i]][[1]])
}
#Plot7 output section
if(InteractMode == TRUE){
message("Show 3 order polynomial fitting plot of log_flow_stress vs. log_strain_rate?")
ptr <- select.list(c("Yes", "No"))
on.exit(select.list(ptr))
if(ptr == "Yes"){
legendvec <- colnames(x)
yscale <- c(min(unlist(scatter_mat)), max(unlist(scatter_mat)))
xscale <- c(min(lgSR), max(lgSR))
labx <- expression(paste("ln(",ring(epsilon), "/s"^-1, ")", sep = ""))
laby <- expression(paste("ln(",sigma, "/MPa)", sep = ""))
for(i in 1:dims[2]){
Coefs <- Coef_mat[i,]
plot(lgSR, scatter_mat[i,], type = "p", col = clrvec[i], pch=16, xlim = xscale, ylim = yscale, xlab = labx, ylab = laby)
curve(Coefs[1]+Coefs[2]*x+Coefs[3]*x^2+Coefs[4]*x^3, xscale[1], xscale[2], col=clrvec[i], ylim = yscale, add = TRUE)
if (i == dims[2]){
legend(legendloc,legend = legendvec, fill = clrvec[1:dims[2]], cex = legendcex, bg = "transparent", box.lty = 0)
cus_par <- par(new=FALSE)
on.exit(par(cus_par))
} else {
cus_par <- par(new=TRUE)
on.exit(par(cus_par))
}
}
}
invisible(readline(prompt="Press [enter] to continue"))
}
m.table <- matrix(NA, nrow = dims[1], ncol = dims[2])
above.table <- m.table
for (i in 1:dims[1]) { # max i = 4
for (j in 1:dims[2]) { # max j = 7
bcd <- Coef_mat[j,c(2:4)] # d3 vector
lgstr <- lgSR[i]
m.table[i,j] <- bcd[1]+2*bcd[2]*lgstr+3*bcd[3]*(lgstr^2) # Add abs?
above.table[i,j] <- 6*lgstr*bcd[3] + 2*bcd[2] # for xi
}
}
temp_mat <- m.table/(1+m.table)
eta_mat <- as.data.frame(2*temp_mat)
rownames(eta_mat) <- rownames(x)
colnames(eta_mat) <- colnames(x)
below.table <- 1/(m.table*(m.table+1)*log(logbase))
xi_mat <- as.data.frame((above.table/below.table) + m.table)
rownames(xi_mat) <- rownames(x)
colnames(xi_mat) <- colnames(x)
if(ConsFunc == FALSE){
paras <- list(base = logbase, epsilon.strain = eps)
}
# package two tables
tables <- list("SRTtable"=x, "etatable"=eta_mat, "xitable"=xi_mat)
result <- list("MaterialCoefficients"=paras, "tablelist"=tables)
class(result) <- "DMMresult"
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.