# Routine for Life-Stress Relationship Plotting Example
# Coded by Reuel Smith 2015-2017
# ========================================================================
# Example Problem 4.7
# Consider an accelerated test where temperature is the accelerating
# variable and 20 units were tested to failure assuming complete failures
# (no censoring). Eight units were tested at 406 K, and six units each at
# 436 K and 466 K, with times to failure tabulated below. Assuming an
# Arrhenius-Weibull life-stress relationship, find the parameters of the
# model, using both the plotting method and MLE method and compare the
# results. Find the expected life at the use level temperature of 353 K
# and compare the results.
# Temperature Time to Failure (hrs)
# =========== =====================================
# 406 K 248 456 528 731 813 965 972 1528
# 436 K 164 176 289 319 386 459
# 466 K 92 105 155 184 219 235
#
# Solution:
# This routine computes the expected life at use level temperature 353 ?K
# and plots the relationship plot. The assumed Life-Stress relationship
# for this problem is the Arrhenius model L(t) = b*exp(a/T) or ln(L(t)) =
# ln(b) + a*(1/T) Eqn(4.40) where T-temperature, a-slope, and ln(b) is the
# intercept. The Life-stress relationship plot may take additional models
# in place of Arrhenius (see textbook).
# ========================================================================
# Known alpha (life) values at 406, 436, and 466 degrees K. These are
# evaluated by the Weibull_Plotting_Multiple_Data_Example_Problem_4_7.R
# R script.
alpha_1 <- 898.698471763768
alpha_2 <- 341.942889990736
alpha_3 <- 187.571102122010
# 1/Temperature
Tempaxis <- c(1/406,1/436,1/466)
# Life axis (hours)
Lifeaxis <- c(alpha_1,alpha_2,alpha_3)
# Fit the X and Y vectors to a linear equation (gets m and b)
p4 <- lm(Lifeaxis ~ poly(Tempaxis, 1, raw=TRUE))
slope <- summary(p4)$coefficients[2,1]
intercept<- summary(p4)$coefficients[1,1]
# Extrapolate best fit relation line for Life and Stress
x1 <- min(Tempaxis)+(((1/300)-min(Tempaxis))/1000)*c(1:1000)
y1 <- intercept + slope*x1
# Value of life at use condition 353 degrees K
Uselife <- intercept + slope*(1/353)
# The Relationship plot
plot(1/353,Uselife, col="black",
xlab="1/Temp (K^-1)", ylab="Characteristic Life (hours)", pch=13,
xlim=c(2e-3,3.4e-3), ylim=c(0,3000), axes=FALSE)
lines(x1, y1, col="black", lty=3)
# STRESS LEVEL 1 PLOT
points(1/406,alpha_1, pch=16, col="blue")
# STRESS LEVEL 2 PLOT
points(1/436,alpha_2, pch=17, col="red")
# STRESS LEVEL 3 PLOT
points(1/466,alpha_3, pch=18, col="green")
axis(1, at=NULL, las=2, cex.axis=0.7)
axis(2, at=NULL, las=2, cex.axis=0.7)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = par("lwd"), equilogs = TRUE)
# The legend heading
# The first two places signify the (X,Y) where the legend will be located. CHANGE the
# legend labels in the LEGEND="" portion. PCH changes the label
# types and LTY changes the line types.
legend(0.002, 3000, legend=c("Stress Level 1 - 406 deg K","Stress Level 2 - 436 deg K","Stress Level 3 - 466 deg K","Characteristic Use life at 353 deg K","Best Fit Line"),
col=c("blue", "red", "green","black","black"), pch=c(16,17,18,13,-1), lty=c(0,0,0,0,3), cex=0.8,
text.font=4, bg='white')
Top