R/DMMprocess.R

In TPMplt: Tool-Kit for Dynamic Materials Model and Thermal Processing Maps

#' 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)
#' print(DMM)
#'
#' \dontrun{
#' # Calculating for constitutive equation but
#' # Without plots printout.
#' DMM <- DMMprocess(epstable, ConsFunc=TRUE)
#' print(DMM)
#'
#' # Calculating for constitutive equation and
#' # required fitting plots printout.
#' 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(class(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){
      cat("Show fitting plot of log_flow_stress vs. log_strain_rate?")
      ptr <- select.list(c("Yes", "No"))
      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)
            par(new=FALSE)
          } else {
            par(new=TRUE)
          }
        }
      }
      cat("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){
      cat("Show fitting plot of flow_stress vs. log_strain_rate?")
      ptr <- select.list(c("Yes", "No"))
      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)
            par(new=FALSE)
          } else {
            par(new=TRUE)
          }
        }
        cat("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){
      cat("Show fitting plot of log[sinh(alpha*flow_stress)] vs. log_strain_rate?")
      ptr <- select.list(c("Yes", "No"))
      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)
            par(new=FALSE)
          } else {
            par(new=TRUE)
          }
        }
        cat("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){
      cat("Show fitting plot of log[sinh(alpha*flow_stress)] vs. 1/T(K^-1)?")
      ptr <- select.list(c("Yes", "No"))
      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)
            par(new=FALSE)
          } else {
            par(new=TRUE)
          }
        }
        cat("Mean value of K calculated from log[sinh(alpha*flow_stress)] vs. 1/T(K^-1) is", K,"\n")
        cat("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) {
      cat("Show fitting plot of log_Z vs. log[sinh(alpha*flow_stress)]?")
      ptr <- select.list(c("Yes", "No"))
      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)

        cat("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) {
      cat("Show fitting plot of Z vs. [sinh(alpha*flow_stress)]^n?")
      ptr <- select.list(c("Yes", "No"))
      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)

        cat("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){
      cat("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){
    cat("Show 3 order polynomial fitting plot of log_flow_stress vs. log_strain_rate?")
    ptr <- select.list(c("Yes", "No"))
    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)
          par(new=FALSE)
        } else {
          par(new=TRUE)
        }
      }
    }
    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)
}