Nothing
## STRESS
## Kyle Elmy and Jim Kaklamanos
## 1 October 2015
########################################################################################
## 1. VERTICAL STRESS AT A POINT
##
## Output: A three-element list containing:
## sigmaZ.eff = effective vertical stress
## sigmaZ.total = total vertical stress
## u = pore water pressure
##
## Input: gamma = vector of unit weights (pcf or kN/m^3)
## thk = vector of layer thicknesses (ft or m)
## depth = vector of layer bottom depths (ft or m)
## zw = depth of groundwater table (ft or m)
## zout = desired depth of output (ft or m): a single value
## gammaW = unit weight of water (default = 62.4 pcf for English units;
## 9.81 kN/m^3 for metric units)
## metric = logical variable: TRUE (for metric units) or FALSE (for English units)
##
## Note: Either layer thicknesses or depths to layer bottoms must be specified.
sigmaZ <- function(gamma, thk = NA, depth = NA, zw, zout, gammaW = NA, metric){
## Calculate layer thicknesses, if not provided
if(all(is.na(thk) == TRUE)){
if(length(depth) == 1){
thk <- depth
} else{
if(length(depth) > 1){
thk <- c(depth[1], diff(depth))
}
}
}
## Define gammaW
if(metric == TRUE){
if(is.na(gammaW) == TRUE){
gammaW <- 9.81
}
} else{
if(metric == FALSE){
if(is.na(gammaW) == TRUE){
gammaW <- 62.4
}
}
}
## Case: zout = 0
if(zout == 0){
sigmaZ.total <- 0
## Case: zout > 0
} else{
## Add ficticious zeroth layer
thk <- c(0, thk)
gamma <- c(0, gamma)
## Depths of layer interfaces
zint <- cumsum(thk)
## Layer above point of interest
Nlayers <- max(which(zint < zout))
## Total stress
sigmaZ.total <- sum(thk[1:Nlayers] * gamma[1:Nlayers]) + (zout - zint[Nlayers]) * gamma[Nlayers+1]
}
## Pore water pressure
u <- gammaW * max(0, zout - zw)
## Effective stress
sigmaZ.eff <- sigmaZ.total - u
return(list(sigmaZ.eff = sigmaZ.eff, sigmaZ.total = sigmaZ.total, u = u))
}
## EXAMPLE:
## sigmaZ(gamma = c(108, 116), depth = c(15, 40), zout = 18, zw = 15, metric = FALSE)
########################################################################################
## 2. VERTICAL STRESS PROFILE
##
## Output: A four-element list containing four vectors:
## depth = depth
## sigmaZ.eff = effective vertical stress
## sigmaZ.total = total vertical stress
## u = pore water pressure
##
## Input: gamma = vector of unit weights (pcf or kN/m^3)
## thk = vector of layer thicknesses (ft or m)
## depth = vector of layer bottom depths (ft or m)
## zw = depth of groundwater table (ft or m)
## zout = desired depths of output (ft or m): defaults to critical locations in
## the profile (the top and bottom of the profile, layer interfaces, and the
## groundwater table) [recommended]
## gammaW = unit weight of water (default = 62.4 pcf for English units;
## 9.81 kN/m^3 for metric units)
## metric = logical variable: TRUE (for metric units) or FALSE (for English units)
sigmaZ.profile <- function(gamma, thk = NA, depth = NA, zw, zout = NA, gammaW = NA, metric){
## Determine critical depths
if(all(is.na(depth) == TRUE)){
depth <- cumsum(thk)
}
if(all(is.na(zout) == TRUE)){
zout <- sort(unique(c(0, depth, zw)))
}
## Calculate stresses
sigmaZ.eff <- vector(length = length(zout))
sigmaZ.total <- vector(length = length(zout))
u <- vector(length = length(zout))
for(i in 1:length(zout)){
temp <- sigmaZ(gamma = gamma, thk = thk, depth = depth, zw = zw,
zout = zout[i], gammaW = gammaW, metric = metric)
sigmaZ.eff[i] <- temp$sigmaZ.eff
sigmaZ.total[i] <- temp$sigmaZ.total
u[i] <- temp$u
}
return(list(depth = zout, sigmaZ.eff = sigmaZ.eff, sigmaZ.total = sigmaZ.total, u = u))
}
## EXAMPLE:
## sigmaZ.profile(gamma = c(108, 116), depth = c(15, 40), zw = 15, metric = FALSE)
########################################################################################
## 3. VERTICAL STRESS PLOT
##
## Output: Plot of vertical stress versus depth
##
## Input: depth = vector of depths (ft or m)
## sigmaZ.eff = vector of effective stresses (psf or kPa)
## sigmaZ.total = vector of total stresses (psf or kPa)
## u = vector of pore water pressures (psf or kPa)
## metric = logical variable: TRUE (for metric units) or FALSE (for English units)
##
## Notes: o If sigmaZ.total and u are left blank, the plot is only constructed for effective stress
## o Once constructed, additional profiles may be added to this plot
## (for example, for induced stress or maximum past pressure)
sigmaZ.plot <- function(depth, sigmaZ.eff, sigmaZ.total = NA, u = NA, metric){
## Expression for axes
if(metric == TRUE){
xLab <- "Stress (kPa)"
yLab <- "Depth, z (m)"
} else{
if(metric == FALSE){
xLab <- "Stress (psf)"
yLab <- "Depth, z (ft)"
}
}
## X axis limit
if(all(is.na(sigmaZ.total)) == FALSE){
xMax <- max(sigmaZ.total)
} else{
xMax <- max(sigmaZ.eff)
}
## Plot
par(mgp = c(2.8, 1, 0), las = 1)
plot(sigmaZ.eff, depth, type = "l", yaxs = "i", xaxt = "n", col = "black", lwd = 2,
xlim = c(0, xMax), ylim = c(max(depth), 0),
xlab = "", ylab = yLab)
axis(3)
mtext(xLab, side = 3, line = 2.5)
if(all(is.na(sigmaZ.total)) == FALSE && all(is.na(u)) == FALSE){
lines(u, depth, col = "blue", lwd = 1)
lines(sigmaZ.total, depth, col = "gray70", lty = 1, lwd = 3)
lines(sigmaZ.eff, depth, col = "black", lwd = 2)
}
## Line for gwt
if(all(is.na(u)) == FALSE){
Zgwt <- depth[max(which(u == 0))]
abline(h = Zgwt, lwd = 1, lty = 3)
points(x = 0.95*xMax, y = Zgwt*0.97, pch = 6, cex = 1.1)
}
## Legend
if(all(is.na(sigmaZ.total)) == FALSE && all(is.na(u)) == FALSE){
legend(x = "topright", bty = "n", lty = c(1, 1, 1), lwd = c(3, 2, 1),
col = c("gray70", "black", "blue"),
legend = c("Total vertical stress", "Effective vertical stress",
"Pore water pressure"))
} else{
legend(x = "topright", bty = "n", lty = c(1), lwd = c(3),
col = c("black"), legend = c("Effective vertical stress"))
}
}
## EXAMPLE:
## Site with constant unit weight = 100 pcf, GWT at 10 ft depth
## temp <- sigmaZ.profile(gamma = rep(100, 3), depth = c(10, 20, 30), zw = 10, metric = FALSE)
## depth <- temp$depth
## sigmaTotal <- temp$sigmaZ.total
## u <- temp$u
## sigmaEff <- temp$sigmaZ.eff
## sigmaZ.plot(depth = depth, sigmaZ.eff = sigmaEff,
## metric = FALSE, sigmaZ.total = sigmaTotal, u = u)
########################################################################################
## 4. INDUCED STRESS DUE TO POINT LOAD
##
## Output: A six-element list containing the induced stresses due to a point load
## following the Boussinesq method:
## o sigmaX = induced horizontal stress in the x direction
## o sigmaY = induced horizontal stress in the y direction
## o sigmaZ = induced vertical stress
## o tauZX = induced shear stress on the ZX plane
## o tauYX = induced shear stress on the YX plane
## o tauYZ = induced shear stress on the YZ plane
##
## Input: P = point load
## x = horizontal distance from load in the X direction
## y = horizontal distance from load in the Y direction
## z = depth of interest
## nu = Poisson's ratio
induced.point <- function(P, x, y, z, nu){
## Compute Pythagorean distances
R <- sqrt((x^2)+(y^2)+(z^2))
r <- sqrt((x^2)+(y^2))
## Factor
f <- P/(2*pi)
## -----------------------------------------------
## Induced horizontal stress in X direction
## First section
S1X <- (3*(x^2)*z)/(R^5)
## Second section
S2X <- ((x^2)-(y^2))/((R*(r^2))*(R+z))
## Third section
S3X <- ((y^2)*z)/((R^3)*(r^2))
## Combine second and third
S4X <- (1-(2*nu))*(S2X+S3X)
## Compute change in stress in the x direction
sigmaX <- f*(S1X-S4X)
## -----------------------------------------------
## Induced horizontal stress in Y direction
## First section
S1Y <- ((3*(y^2)*z)/(R^5))
## Second section
S2Y <-((y^2)-(x^2))/((R*(r^2))*(R+z))
## Third section
S3Y <- ((x^2)*z)/((R^3)*(r^2))
## Combine second and third
S4Y <- (1-(2*nu))*(S2Y+S3Y)
## Compute change in stress in the y direction
sigmaY <- f*(S1Y-S4Y)
## -----------------------------------------------
## Induced vertical stress in Z direction
sigmaZ <- (3*P*(z^3))/(2*pi*(R^5))
## -----------------------------------------------
## Induced shear stress in the ZX direction
tauZX <- (3*P*(z^2)*x)/(2*pi*(R^5))
## -----------------------------------------------
## Induced shear stress in the YX direction
## First section
S1yx <- ((3*x*y*z)/(R^5))
## Second Section
S2yx <- (((2*R)+2)*x*y)/((R+(z^2))*(R^3))
## Combine first and second section
S3yx <- (1-(2*nu))*S2yx
## Compute TauYX
tauYX <- f*(S1yx-S3yx)
## -----------------------------------------------
## Induced shear stress in the YZ direction
tauYZ <- (3*P*(z^2)*y)/(2*pi*(R^5))
## Return elements
return(list(sigmaX = sigmaX, sigmaY = sigmaY, sigmaZ = sigmaZ,
tauZX = tauZX, tauYX = tauYX, tauYZ = tauYZ))
}
## EXAMPLE:
## induced.point(P = 100000, x = 5, y = 2, z = 6, nu = 0.35)
########################################################################################
## 5. INDUCED STRESS DUE TO POINT LOAD: PROFILE
##
## Output: A seven-element list containing seven vectors displaying the depth variation of
## induced stresses due to a point load following the Boussinesq method:
## o depth = vector of depths
## o sigmaX = induced horizontal stress in the x direction
## o sigmaY = induced horizontal stress in the y direction
## o sigmaZ = induced vertical stress
## o tauZX = induced shear stress on the ZX plane
## o tauYX = induced shear stress on the YX plane
## o tauYZ = induced shear stress on the YZ plane
##
## Input: P = point load
## x = horizontal distance from load in the X direction
## y = horizontal distance from load in the Y direction
## z = vector of depths of interest (default: 1-ft or 1-m increments,
## to a maximum depth of 50 ft or 50 m)
## nu = Poisson's ratio
induced.point.profile <- function(P, x, y, z = NA, nu){
## Specify default values for depths
if(all(is.na(z) == TRUE)){
z <- seq(from = 0, to = 50, by = 1)
}
## Calculate stresses
sigmaX <- vector(length = length(z))
sigmaY <- vector(length = length(z))
sigmaZ <- vector(length = length(z))
tauZX <- vector(length = length(z))
tauYX <- vector(length = length(z))
tauYZ <- vector(length = length(z))
for(i in 1:length(z)){
temp <- induced.point(P = P, x = x, y = y, z = z[i], nu = nu)
sigmaX[i] <- temp$sigmaX
sigmaY[i] <- temp$sigmaY
sigmaZ[i] <- temp$sigmaZ
tauZX[i] <- temp$tauZX
tauYX[i] <- temp$tauYX
tauYZ[i] <- temp$tauYZ
}
## Return elements
return(list(depth = z, sigmaX = sigmaX, sigmaY = sigmaY, sigmaZ = sigmaZ,
tauZX = tauZX, tauYX = tauYX, tauYZ = tauYZ))
}
## EXAMPLE:
## induced.point.profile(P = 100000, x = 5, y = 2, nu = 0.35)
########################################################################################
## 6. INDUCED STRESS DUE TO AREA LOAD
##
## Output: Induced vertical stress at the center of the loaded area
##
## Input: z = depth of interest
## q = applied pressure at ground surface
## B = width of loaded area
## L = length of loaded area (rectangular foundations only)
## shape = shape of loaded area (a string containing
## "circle", "square", "strip", or "rectangle")
##
## Notes: This function currently uses the approximate method of Poulos and Davis (1974).
## More advanced formulations are expected in future versions of this package.
induced.area <- function(z, q, B, L = NA, shape){
## Check input for shape
if(shape != "circle" && shape != "square" && shape != "strip" && shape != "rectangle"){
stop("Shape must be either 'circle', 'square', 'strip', or 'rectangle'.")
}
## Compute the influence factor for circular loaded areas
IC <- 1-(1/(1+(B/(2*z))^(2))^(1.5))
## Compute the influence factor for square loaded areas
IS <- 1-(1/(1+(B/(2*z))^(2))^(1.76))
## Compute the influence factor for square loaded areas
Ist <- 1-(1/(1+(B/(2*z))^(1.38))^(2.60))
## Compute the influence factor for rectangular loaded areas of width B and length L
if(is.na(L) == FALSE){
EXPI <- 1.38+((0.62*B)/L)
EXPE <- 2.60-((0.84*B)/L)
IR <- 1-(1/(1+((B/(2*z))^(EXPI)))^(EXPE))
}
## Induced stress
if(shape == "circle"){
I <- IC
} else{
if(shape == "square"){
I <- IS
} else{
if(shape == "strip"){
I <- Ist
} else{
if(shape == "rectangle"){
I <- IR
}
}
}
}
return(I * q)
}
## EXAMPLE:
## induced.area(z = 10, q = 1000, B = 3, shape = "square")
########################################################################################
## 7. INDUCED STRESS DUE TO AREA LOAD: PROFILE
##
## Output: A two-element list:
## depth = vector of depths
## sigmaZ = vector of induced vertical stresses below center of the loaded area
##
## Input: z = vector of depths (default: 1-ft or 1-m increments,
## to a maximum depth of 50 ft or 50 m)
## q = applied pressure at ground surface
## B = width of loaded area
## L = length of loaded area (rectangular foundations only)
## shape = shape of loaded area (a string containing
## "circle", "square", "strip", or "rectangle")
##
## Notes: This function currently uses the approximate method of Poulos and Davis (1974).
## More advanced formulations are expected in future versions of this package.
induced.area.profile <- function(z = NA, q, B, L = NA, shape){
## Specify default values for depths
if(all(is.na(z) == TRUE)){
z <- seq(from = 0, to = 50, by = 1)
}
## Calculate stresses
sigmaZ <- vector(length = length(z))
for(i in 1:length(z)){
sigmaZ[i] <- induced.area(z = z[i], q = q, B = B, L = L, shape = shape)
}
## Return values
return(list(depth = z, sigmaZ = sigmaZ))
}
## EXAMPLE:
## induced.area.profile(q = 1000, B = 3, shape = "square")
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.