Rでムーディー線図

library(ggplot2)

# 描画の条件
re_max <- 1e8 # 右側の縦軸の基準にするレイノルズ数
rr_values <- c(
  1e-5, 2e-5, 5e-5, 1e-4, 2e-4, 2.5e-4, 5e-4, 1e-3, 2e-3, 2.5e-3,
  5e-3, 7.5e-3, 1e-2, 1.5e-2, 2e-2, 3e-2, 4e-2, 5e-2
) # 相対粗度 epsilon / D_H
re_turb <- 10^seq(log10(4000), log10(re_max), length.out = 200) # 乱流域のレイノルズ数
re_lam_solid <- 10^seq(log10(600), log10(2300), length.out = 50) # 層流域(実線)のレイノルズ数
re_lam_dash <- 10^seq(log10(2300), log10(4000), length.out = 50) # 臨界域(破線)のレイノルズ数

# 相対粗度ごとの曲線(D_H = 1 とすれば epsilon が相対粗度になる)
curves <- do.call(rbind, lapply(rr_values, function(rr) {
  data.frame(
    Re = re_turb,
    lambda = colebrook_white(epsilon = rr, D_H = 1, Re = re_turb),
    rr = rr
  )
}))

# 滑面管の曲線(epsilon = 0)と、層流の直線(lambda = 64 / Re)
smooth <- data.frame(Re = re_turb, lambda = colebrook_white(epsilon = 0, D_H = 1, Re = re_turb))
laminar_solid <- data.frame(Re = re_lam_solid, lambda = 64 / re_lam_solid)
laminar_dash <- data.frame(Re = re_lam_dash, lambda = 64 / re_lam_dash)

# 右側の縦軸用の変換: Re = re_max での lambda から epsilon / D_H を求める
lambda_to_rr <- function(lambda) {
  suppressWarnings(colebrook_white(lambda = lambda, D_H = 1, Re = re_max))
}

p <- ggplot() +
  # 臨界域(層流から乱流への遷移域)
  annotate("rect", xmin = 2300, xmax = 4000, ymin = 0.006, ymax = 0.11, fill = "grey80", alpha = 0.5) +
  annotate("text", x = sqrt(2300 * 4000), y = 0.0068, label = "臨界域", size = 3.5) +
  annotate("text", x = sqrt(2300 * 4000), y = 0.065, label = "Re = 2300〜4000", size = 3, angle = 90) +
  # 相対粗度ごとの曲線
  geom_line(data = curves, aes(x = Re, y = lambda, group = rr), colour = "steelblue", linewidth = 0.5) +
  # 滑面管
  geom_line(data = smooth, aes(x = Re, y = lambda), colour = "black", linetype = "dashed", linewidth = 0.6) +
  annotate("text", x = 3e6, y = 0.0072, label = "滑面管", size = 3.5, hjust = 0) +
  # 層流
  geom_line(data = laminar_solid, aes(x = Re, y = lambda), colour = "firebrick", linewidth = 0.8) +
  geom_line(data = laminar_dash, aes(x = Re, y = lambda), colour = "firebrick", linewidth = 0.8, linetype = "dashed") +
  annotate("text", x = 640, y = 0.02, label = "層流\nλ=64/Re", size = 3.2, hjust = 0, colour = "firebrick") +
  scale_x_log10(
    breaks = 10^(3:8), # 目盛りの文字は 10^n の位置だけにする
    minor_breaks = as.vector(outer(1:9, 10^(2:7))), # 1〜9 x 10^n の位置に補助線を引く
    labels = function(x) parse(text = paste0("10^", round(log10(x)))), # 指数表記(10^3, 10^4, ...)
    limits = c(600, re_max),
    expand = c(0, 0)
  ) +
  scale_y_log10(
    breaks = c(0.006, 0.008, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1),
    labels = function(x) format(x, scientific = FALSE, trim = TRUE),
    limits = c(0.006, 0.11),
    oob = scales::oob_keep,
    expand = c(0, 0),
    sec.axis = sec_axis(
      lambda_to_rr,
      name = expression(paste("相対粗度 ", epsilon / D[H])), # 数式表記(添え字の H を下付きにする)
      breaks = rr_values,
      labels = function(x) format(x, scientific = FALSE, trim = TRUE, drop0trailing = TRUE)
    )
  ) +
  labs(
    title = "ムーディー線図(コールブルック・ホワイトの式)",
    x = "レイノルズ数 Re",
    y = "管摩擦係数 λ"
  ) +
  theme_bw() +
  theme(
    panel.grid.minor.x = element_line(linewidth = 0.2, colour = "grey90"),
    panel.grid.minor.y = element_blank(),
    axis.text.y.right = element_text(size = 8)
  )

print(p)

# 計算する条件(水が流れる円管の例)
cases <- data.frame(
  name = c("A: 鋼管", "B: 塩ビ管", "C: コンクリート管"),
  epsilon = c(0.045e-3, 0.0015e-3, 1e-3), # 絶対粗度 [m]
  D_H = c(0.1, 0.05, 0.5), # 水力直径(満水の円管では内径に等しい) [m]
  v = c(2, 1, 3), # 平均流速 [m/s]
  nu = 1.0e-6 # 水の動粘度 [m^2/s]
)

# レイノルズ数、相対粗度、管摩擦係数を求める
cases$Re <- cases$v * cases$D_H / cases$nu
cases$rr <- cases$epsilon / cases$D_H
cases$lambda <- colebrook_white(epsilon = cases$epsilon, D_H = cases$D_H, Re = cases$Re)
cases$label <- sprintf(
  "%s\nRe = %s\nλ = %.4f",
  cases$name, format(cases$Re, big.mark = ",", scientific = FALSE, trim = TRUE), cases$lambda
)
print(cases[, c("name", "epsilon", "D_H", "v", "Re", "rr", "lambda")])

# ラベルの位置(A, B, C の順)。重ならないように、1つずつ指定する
# A, B は曲線の下側の空いている場所に、各点の縦の補助線に沿わせて置く
cases$label_x <- cases$Re * c(1.1, 0.9, 1.2) # ラベルの横位置(Re の 1.1 倍など)
cases$label_y <- c(0.0145, 0.0165, cases$lambda[3] * 1.25) # ラベルの縦位置
cases$label_hjust <- c(0, 1, 0) # 0: ラベルの左端を合わせる、1: 右端を合わせる

# 各条件の相対粗度に対応する曲線(ムーディー線図の曲線の間を補う)
curves_cases <- do.call(rbind, lapply(seq_len(nrow(cases)), function(i) {
  data.frame(
    Re = re_turb,
    lambda = colebrook_white(epsilon = cases$rr[i], D_H = 1, Re = re_turb),
    name = cases$name[i]
  )
}))

p2 <- p +
  # 各条件の曲線
  geom_line(data = curves_cases, aes(x = Re, y = lambda, group = name), colour = "darkorange", linewidth = 0.5) +
  # 点から左側の縦軸と横軸へ、読み取り用の補助線
  geom_segment(
    data = cases, aes(x = 600, xend = Re, y = lambda, yend = lambda),
    colour = "grey40", linetype = "dotted"
  ) +
  geom_segment(
    data = cases, aes(x = Re, xend = Re, y = 0.006, yend = lambda),
    colour = "grey40", linetype = "dotted"
  ) +
  # 計算した点とラベル
  geom_point(data = cases, aes(x = Re, y = lambda), colour = "red", size = 3) +
  geom_label(
    data = cases, aes(x = label_x, y = label_y, label = label, hjust = label_hjust),
    size = 3, fill = "white"
  )

print(p2)