Before you start, have you read the introduction to the Ternary package?

A ternary plot can be shaded and contoured according to a known function:

library("Ternary")
par(mar = rep(0.2, 4))

FunctionToContour <- function(a, b, c) {
  a - c + (4 * a * b) + (27 * a * b * c)
}

values <- TernaryPointValues(FunctionToContour, resolution = 24L)
TernaryPlot(alab = "a", blab = "b", clab = "c",
            # Place an opaque fill behind grid lines:
            panel.first = ColourTernary(values, spectrum = hcl.colors(256)))
TernaryContour(FunctionToContour, resolution = 36L, legend = TRUE, 
               legend... = list(bty = "n"))

If it is not computationally feasible to execute a function at every point, it is possible to interpolate between known values.

Inverse distance weighting is a simple approach:

# If using your own data, set
# abc <- [Three-column matrix containing a, b, c coordinates of points]
# response <- [Vector of values at the points specified in abc]

# We'll use random data for this example.

# First we'll generate some points at which to sample:
set.seed(0)
nPts <- 50
a <- runif(nPts)
b <- runif(nPts) * (1 - a)
c <- 1 - a - b
abc <- rbind(a, b, c)

# Then we'll measure the response variable at each of those points:
response <- FunctionToContour(a, b, c)


# Now we must start a plot, to define the coordinate system
par(mar = rep(0.2, 4))
TernaryPlot(alab = "a", blab = "b", clab = "c")

# Convert measured points to XY
xy <- TernaryToXY(abc)

# Use an inverse distance weighting to interpolate between measured points
Predict <- function(predXY) {
  Distance <- function(a, b) {
    apply(a, 2, function(pt) sqrt(colSums((pt - b) ^ 2)))
  }
  dists <- Distance(xy, predXY)
  id <- 1 / dists
  idw <- id / rowSums(id)

  # Return:
  colSums(response * t(idw))
}

# Predict at triangle centres
tri <- TriangleCentres(resolution = 12L)
# Adjust the resolution to suit your own dataset

# Now we interpolate between our known values to generate a colour for each
# of our tiles
predicted <- Predict(tri[1:2, ])
map <- rbind(x = tri["x", ], y = tri["y", ], z = predicted,
             down = tri["triDown", ])

# Place a semitransparent colour fill over grid lines:
ColourTernary(map)

# Calculate contours
PredictABC <- function(a, b, c) Predict(TernaryToXY(rbind(a, b, c)))
TernaryContour(PredictABC, resolution = 36L, legend = 6, 
               legend... = list(bty = "n"))

# Mark the points at which we took measurements
TernaryPoints(abc, pch = 3, col = "#cc3333")