14 phút để đọc

Tổng Quan

Mô hình Cam Clay cải tiến - Modified Cam-Clay (MCC) là một trong những mô hình đàn dẻo kinh điển nhất trong cơ học đất, đặc biệt hiệu quả để dự báo hành vi của đất sét cố kết bình thường (normally consolidated clay) và đất sét quá cố kết(over-consolidated clay) dưới các điều kiện tải khác nhau.

Bài viết này trình bày phương pháp triển khai mô hình MCC mô phỏng biến dạng triaxial drained/undrained một cách nhanh chóng và chính xác.

Ứng Dụng Thực Tế:

  • Phân tích ổn định nền móng
  • Dự báo chuyển động đất đai
  • Thiết kế hệ thống thoát nước cho công trình
  • Mô phỏng động học đất dưới tải động

1. Tham Số Mô Hình Modified Cam-Clay

1.1 Các Thông Số Cơ Bản

Mô hình MCC được kiểm soát bởi 7 thông số địa kỹ thuật chính:

Thông Số Ký Hiệu Mô Tả Đơn Vị Giá Trị Mẫu
Áp suất tiền cố kết ban đầu \(p_c\) Preconsolidation pressure kPa 150
Áp suất gối tựa ban đầu \(p_0\) Initial confining pressure kPa 150
Độ dốc đường tráng thái tới hạn- tham số ma sát giới hạn \(M\) Critical friction parameter - 0.95
Chỉ số nén \(\lambda\) Compression index - 0.2
Chỉ số nén/nở lại \(\kappa\) Recompression index - 0.04
Tham số đường NCL \(N\) Normal Consolidation Line parameter - 2.5
Hệ số Poisson \(\nu\) Poisson’s ratio - 0.15

1.2 Các khái nhiệm cính

Ứng Suất Hữu Hiệu (\(p'\) hoặc \(p\))
Ứng suất hữu hiệu là sự kết hợp của ứng suất tổng và áp lực nước lỗ rỗng: \(p' = \frac{\sigma_1' + \sigma_2' + \sigma_3'}{3}\)

Ứng Suất Lệch (\(q\))
Độ lệch ứng suất chính được định nghĩa trong thử nghiệm triaxial: \(q = \sigma_1' - \sigma_3'\)

Hệ Số Rỗng (\(e\)) và Thể Tích Riêng (\(V\))
\(e = V - 1.0\) \(V = N - \lambda \ln(p_c) + \kappa \ln\left(\frac{p_c}{p_0}\right)\)

Mặt Chảy Modified Cam-Clay
Bề mặt chảy trong không gian \(p'-q\) được mô tả bởi: \(f = q^2 - M^2 p'(p_c - p') = 0\)

2. Cấu trúc Code Julia

Triển khai mô hình MCC được chia thành 6 phần chính:

2.1 Nhập Thư Viện Cần Thiết

using Plots, LinearAlgebra, DelimitedFiles, Statistics, FileIO
gr()  # Use GR backend for faster plotting

Thư Viện:

  • Plots.jl: Tạo hình vẽ khoa học với GR backend nhanh
  • LinearAlgebra: Phép toán ma trận và vectơ
  • DelimitedFiles: Đọc/ghi tệp CSV
  • Statistics: Hàm thống kê
  • FileIO: Xử lý tệp và định dạng

2.2 Khởi Tạo Thông Số

# Input Parameters for Modified Cam-Clay
println("Input Parameters for Modified Cam-Clay:")

# Geotechnical Parameters
cp = 150.0          # Initial consolidation pressure (kPa)
p₀ = 150.0          # Initial confining pressure (kPa)
M = 0.95            # Critical friction angle parameter
λ = 0.2             # Compression index
κ = 0.04            # Recompression index
N = 2.5             # Parameter for NCL
ν = 0.15            # Poisson's ratio

# Analysis type: 1 = Drained, 2 = Undrained
analysis = 1        # 1 = Triaxial Drained, 2 = Triaxial Undrained

# Iteration parameters
iteration = 7500    # Number of strain increments
ide = 0.01          # Strain increment (%)
de = ide / 100      # Strain increment (decimal)

Giải Thích:

  • cp: Áp suất tiền cố kết (kPa) - đặc trưng cho lịch sử cố kết của đất
  • p₀: Áp suất gối tựa ban đầu (kPa) - phải bằng cp cho đất cố kết bình thường
  • M: Độ dốc của đường trạng thái tới hạn (Critical State Line)
  • λ, κ: Kiểm soát độ nén của đất trong không gian e-ln(p’)
  • analysis: Chọn loại thử nghiệm (1=thoát nước, 2=không thoát nước)
  • iteration: Số bước tính toán (7500 → biến dạng trục ≈ 75%)

2.3 Khởi Tạo Biến Trạng Thái

# Compute initial state variables
pc = cp  # Current preconsolidation pressure

# Specific Volume and initial void ratio
V = N - λ * log(pc) + κ * log(pc / p₀)  # Specific Volume
e₀ = V - 1.0                             # Initial void ratio

# Over-Consolidation Ratio
OCR = cp / p₀

println("Computed initial values:")
println("Specific Volume V = $$V")
println("Initial void ratio e₀ = $$e₀")
println("Over-Consolidation Ratio (OCR) = $$OCR")

# Memory allocation for results
p = zeros(iteration)      # Mean effective stress
q = zeros(iteration)      # Deviatoric stress
u = zeros(iteration)      # Pore water pressure
void = zeros(iteration)   # Void ratio
epsV = zeros(iteration)   # Volumetric strain
epsD = zeros(iteration)   # Deviatoric strain
es = zeros(iteration)     # Axial strain array

Khởi Tạo OCR:

  • Nếu OCR = 1: Đất cố kết bình thường (Normally Consolidated)
  • Nếu OCR > 1: Đất cố kết quá cố cố kết (Over-Consolidated)
  • Nếu OCR < 1: Trạng thái tiền cố kết (không phổ biến)

2.4 Vòng Lặp Ứng Suất-Biến Dạng

Vòng lặp chính thực hiện các bước tính toán, mỗi bước bao gồm:

  1. Xác định Mô-đun Đàn Hồi \(K = \frac{V \cdot p}{κ}, \quad G = \frac{3K(1-2ν)}{2(1+ν)}\)

  2. Ma Trận Độ Cứng Đàn Hồi (6×6)
    • Đường chéo chính: \(K + \frac{4G}{3}\)
    • Ngoài đường chéo: \(K - \frac{2G}{3}\)
    • Phần cắt: \(G\)
  3. Kiểm Tra Trạng Thái Chảy
    • Nếu \(f < 0\): Sử dụng ma trận đàn hồi
    • Nếu \(f ≥ 0\): Sử dụng ma trận đàn dẻo
  4. Xác định Biến Dạng Tương Ứng
    • Triaxial Drained: \(d\varepsilon_1 = de\), \(d\varepsilon_2 = d\varepsilon_3\) từ điều kiện biên
    • Triaxial Undrained: \(d\varepsilon_1 = de\), \(d\varepsilon_2 = d\varepsilon_3 = -\frac{de}{2}\)
  5. Cập Nhật Ứng Suất \(d\sigma = \mathbf{D} \cdot d\varepsilon\)
for a in 2:iteration
    # Bulk and shear moduli
    K = V * p[a-1] / κ
    G = (3 * K * (1 - 2 * ν)) / (2 * (1 + ν))
    
    # Update preconsolidation pressure on yield surface
    if yield == 0
        pc_current = (q[a-1]^2 / M^2 + p[a-1]^2) / p[a-1]
    else
        pc_current = cp
    end
    
    # Elastic stiffness matrix (6x6)
    De = zeros(6, 6)
    
    # Fill elastic stiffness
    for m in 1:3
        for n in 1:3
            if m == n
                De[m, n] = K + 4 * G / 3
            else
                De[m, n] = K - 2 * G / 3
            end
        end
    end
    
    for m in 4:6
        De[m, m] = G
    end
    
    # Determine stiffness matrix (elastic or plastic)
    if yield < 0
        # Elastic stiffness
        D = De
    else
        # Elastoplastic stiffness - compute gradients
        dfds = zeros(6)
        dfdep = zeros(6)
        
        for m in 1:3
            dfds[m] = (2 * p[a-1] - pc_current) / 3 + 3 * (S[m] - p[a-1]) / M^2
            dfdep[m] = (-p[a-1]) * pc_current * (1 + e₀) / (λ - κ) * 1
        end
        
        # Elastoplastic matrix
        numerator = De * dfds * dfds' * De
        denominator = -(dfdep' * dfds)[1] + (dfds' * De * dfds)[1]
        
        if abs(denominator) > 1e-12
            D = De - numerator ./ denominator
        else
            D = De
        end
    end
    
    # Strain increment based on analysis type
    dStrain = zeros(6)
    if analysis == 1
        # Triaxial Drained
        dStrain[1] = de
        if abs(D[2,2] + D[2,3]) > 1e-12
            dStrain[2] = -1 * D[2,1] / (D[2,2] + D[2,3]) * de
        else
            dStrain[2] = -de / 2
        end
        if abs(D[3,2] + D[3,3]) > 1e-12
            dStrain[3] = -1 * D[3,1] / (D[3,2] + D[3,3]) * de
        else
            dStrain[3] = -de / 2
        end
    else
        # Triaxial Undrained
        dStrain[1] = de
        dStrain[2] = -de / 2
        dStrain[3] = -de / 2
    end
    
    # Stress increment and update
    dS = D * dStrain
    S = S + dS
    strain = strain + dStrain
    
    # Update invariants and state variables
    depsV = dStrain[1] + dStrain[2] + dStrain[3]
    depsD = (2 / 3) * (dStrain[1] - dStrain[3])
    
    V = N - λ * log(pc_current) + κ * log(pc_current / p[a-1])
    
    p[a] = (S[1] + S[2] + S[3]) / 3
    q[a] = S[1] - S[3]
    u[a] = p₀ + q[a] / 3 - p[a]
    
    void[a] = V - 1.0
    epsV[a] = epsV[a-1] + depsV
    epsD[a] = epsD[a-1] + depsD
    
    # Update yield surface
    if yield < 0
        yield = q[a]^2 + M^2 * p[a]^2 - M^2 * p[a] * pc_current
    else
        yield = 0
    end
end

3. Kết Quả Mô Phỏng

3.1 Đặc Tính Mặt Chảy (Yield Surface)

Đường Cố Kết Bình Thường (Normal Consolidation Line - NCL)
Trong không gian e-ln(p’), đất sét cố kết bình thường sẽ nằm trên: \(e = N - \lambda \ln(p')\)

Đường Trạng Thái Tới Hạn (Critical State Line - CSL)
Khi đất chảy hết sức (critical state), nó tuân theo: \(e = \Gamma - \lambda \ln(p')\)

Mặt Chảy Modified Cam-Clay
Trong không gian \(p'-q\), bề mặt chảy có dạng ellipse: \(q = M \cdot p' \quad \text{(Critical State Line slope)}\) \(f = q^2 - M^2 p'(p_c - p') = 0\)

# Initial Yield Surface
p_ini_yield = collect(0:pc/100:pc)
q_ini_yield = M .* p_ini_yield
qy_ini = sqrt.(M^2 .* (cp .* p_ini_yield .- p_ini_yield.^2))

# Normal Consolidation Line (NCL)
pNCL = collect(1:pc)
eNCL = N .- λ .* log.(pNCL) .- 1

# Final Yield Surface
p_fyield = collect(0:pc/100:pc)
q_fyield = M .* p_fyield
qyf = sqrt.(M^2 .* (pc .* p_fyield .- p_fyield.^2))

# Critical State Line (CSL)
pCSL = pNCL
Gamma = 1 + void[end] + λ * log(p[end])
eCSL = Gamma .- λ .* log.(pCSL) .- 1

3.2 Phân Tích Toàn Diện

Biểu Đồ 1: Ứng Suất Lệch vs Biến Dạng Trục

\(\text{Trục X}: \varepsilon_a (\%) = \frac{\Delta L}{L_0} \times 100\) \(\text{Trục Y}: q = \sigma_1' - \sigma_3' \text{ (kPa)}\)

Ý Nghĩa Vật Lý:

  • Giai đoạn 0-25%: Tuyến tính đàn hồi, độ cứng cao
  • Giai đoạn 25-50%: Chuyển tiếp, độ cứng giảm
  • Giai đoạn 50%+: Chảy dẻo, độ cứng lớn giảm, tiến tới trạng thái tới hạn

Biểu Đồ 2: Quá trình ứng suất (p-q Plane)

Đường dẫn ứng suất (stress path) là quỹ tích các điểm ứng suất từ trạng thái ban đầu đến cuối:

\(\text{Trục X}: p' = \frac{\sigma_1' + 2\sigma_3'}{3}\) \(\text{Trục Y}: q = \sigma_1' - \sigma_3'\)

Các Đường Tham Chiếu:

  • Mặt Chảy Ban Đầu: Ellipse với \(p_c = 150\) kPa
  • Mặt Chảy Cuối Cùng: Ellipse sau khi cố kết tăng (nếu có)
  • Đường Trạng Thái Tới Hạn: Đường \(q = M \cdot p'\) (độ dốc = 0.95)

Biểu Đồ 3: Biến dạng thể tích và Áp lực nước lỗ rỗng

Triaxial Drained: \(\text{Trục X}: \varepsilon_v = \varepsilon_1 + \varepsilon_2 + \varepsilon_3\) \(\text{Trục Y}: p' \text{ (kPa)}\)

Thể hiện sự thay đổi áp suất hiệu với sự giãn nở/co lại thể tích.

Triaxial Undrained: \(\text{Trục X}: \varepsilon_a (\%)\) \(\text{Trục Y}: u = p_0 + \frac{q}{3} - p' \text{ (kPa)}\)

Thể hiện sự phát triển áp lực nước lỗ rỗng trong điều kiện không thoát nước.

Biểu Đồ 4: Hệ số rỗng vs Log(p’)

\(\text{Trục X}: \log_{10}(p') \text{ (kPa)}\) \(\text{Trục Y}: e = V - 1.0\)

Cho thấy quá trình nén của đất được kiểm soát bởi thông số \(\lambda\) và \(\kappa\).

# Create comprehensive 2x2 subplot layout
fig = plot(size=(1200, 1000), layout=(2, 2), legend=:topright)

# Subplot 1: Deviatoric Stress vs Axial Strain
plot!(1, es, q, label="q vs εₐ", 
      xlabel="Axial strain, εₐ (%)", ylabel="Deviatoric Stress, q (kPa)", 
      title="Deviatoric Stress vs Axial Strain", linewidth=2)

# Subplot 2: Stress Path (p-q plane with yield surfaces)
if analysis == 1
    plot!(2, p, q, label="Stress Path", linewidth=2,
          xlabel="Mean Stress, p (kPa)", ylabel="Deviatoric Stress, q (kPa)",
          title="Stress Path in p-q Plane (Drained)")
    plot!(2, p_ini_yield, qy_ini, label="Initial Yield", linewidth=2, linestyle=:dash)
    plot!(2, p_fyield, qyf, label="Final Yield", linewidth=2, linestyle=:dash)
else
    plot!(2, p, q, label="Stress Path", linewidth=2,
          xlabel="Mean Stress, p (kPa)", ylabel="Deviatoric Stress, q (kPa)",
          title="Stress Path in p-q Plane (Undrained)")
    plot!(2, p_ini_yield, q_ini_yield, label="CSL", linewidth=2, linestyle=:dash)
    plot!(2, p_ini_yield, qy_ini, label="Initial Yield", linewidth=2, linestyle=:dash)
    plot!(2, p_fyield, qyf, label="Final Yield", linewidth=2, linestyle=:dash)
end

# Subplot 3: Volumetric/Undrained Response
if analysis == 1
    plot!(3, epsV, p, label="εᵥ vs p", linewidth=2,
          xlabel="Volumetric Strain, εᵥ (%)", ylabel="Mean Effective Stress, p (kPa)",
          title="Volumetric Response (Drained)")
else
    plot!(3, es, u, label="u vs εₐ", linewidth=2,
          xlabel="Axial strain, εₐ (%)", ylabel="Excess Pore Water Pressure, u (kPa)",
          title="Pore Pressure Response (Undrained)")
end

# Subplot 4: Void Ratio vs Log(p) - NCL and CSL
plot!(4, p, void, label="Test Path", linewidth=2, xscale=:log10,
      xlabel="Mean Effective Stress (kPa)", ylabel="Void Ratio, e",
      title="Void Ratio vs Log(p)")
plot!(4, pNCL, eNCL, label="NCL", linewidth=2, linestyle=:dash, xscale=:log10)
plot!(4, pCSL, eCSL, label="CSL", linewidth=2, linestyle=:dash, xscale=:log10)

fig

3.3 Xuất Kết Quả

Mô phỏng xuất ra các kết quả dưới dạng tệp CSV và PNG:

# Create Output directory if it doesn't exist
output_dir = "Output"
if !isdir(output_dir)
    mkdir(output_dir)
end

println("\nSaving results to $$output_dir/...")

# Save plots
savefig(fig, joinpath(output_dir, "MCC_Results_2x2.png"))
println("✓ Saved: MCC_Results_2x2.png")

# Prepare data for export
data_export = hcat(es, p, q, u, void, epsV, epsD)

# Save data to CSV
open(joinpath(output_dir, "MCC_Results.csv"), "w") do io
    writedlm(io, ["Strain(%)" "p(kPa)" "q(kPa)" "u(kPa)" "void_ratio" "epsV(%)" "epsD(%)"], ',')
    writedlm(io, data_export, ',')
end
println("✓ Saved: MCC_Results.csv")

# Save individual result arrays
for (name, data) in [("stress_p", p), ("stress_q", q), ("void_ratio", void), 
                     ("strain_axial", es), ("pore_pressure_u", u)]
    open(joinpath(output_dir, "$$(name).csv"), "w") do io
        writedlm(io, data)
    end
    println("✓ Saved: $$(name).csv")
end

Tệp Xuất:

  • MCC_Results_2x2.png: Biểu đồ 2×2 tổng hợp
  • MCC_Results.csv: Bảng dữ liệu đầy đủ (7500 hàng × 7 cột)
  • stress_p.csv, stress_q.csv: Ứng suất hữu hiệu và lệch theo từng bước
  • void_ratio.csv: Hệ số rỗng theo quá trình tải
  • strain_axial.csv, pore_pressure_u.csv: Biến dạng trục và áp lực nước lỗ rỗng

4. Phân tích kết quả

4.1 Phân tích Ứng suất - Biến dạng (q-εₐ)

Đường cong q-εₐ thường có ba giai đoạn:

  1. Giai đoạn 0-2% (Đàn hồi): Độ dốc cao, hành vi gần như tuyến tính
    • Đất chưa chảy, chỉ chịu nén đàn hồi
    • \(K\) và \(G\) là hằng số
  2. Giai đoạn 2-25% (Chuyển tiếp): Độ dốc giảm dần
    • Đất bắt đầu chảy (yield)
    • Độ cứng giảm dần từ \(E_0\) (stiffness ban đầu) đến \(E_f\) (cuối cùng)
  3. Giai đoạn 25%+ (Chảy dẻo): Độ dốc gần như bằng không
    • Đất tuân theo quy luật chảy dẻo (flow rule)
    • Tiến tới trạng thái tới hạn (Critical State)

4.2 Đường quá trình ứng suất (p-q)

Trong phương pháp p-q, đường quá trình ứng suất cho biết:

  • Vị trí ban đầu: \((p_0, 0)\)
  • Hình dạng: Tùy thuộc vào loại thử nghiệm
    • Drained: Đường quá trình có thể quay về trái (p’ giảm do thể tích nở)
    • Undrained: Đường trực tiếp hướng lên trên (p’ tăng do u tăng)
  • Điểm cuối cùng: Gần hoặc trên đường CSL

4.3 Biến dạng thể tích (εᵥ-p’)

Triaxial Drained:

  • Đất nén lại ở đầu (co thể tích, εᵥ < 0) vì áp suất tăng
  • Sau đó giãn nở (εᵥ > 0) khi chảy dẻo vì tiến tới trạng thái tới hạn
  • Đây gọi là Dilatancy - sự giãn nở trong quá trình cắt dẻo

Triaxial Undrained:

  • Thể tích tổng không đổi (εᵥ_total = 0)
  • Áp lực nước lỗ rỗng tăng vì không thể thoát nước
  • Ứng suất hữu hiệu p’ có thể giảm → phát triển “vòi phun” (liquefaction)

4.4 Hệ số rỗng vs Log(p’)

Quan hệ e-ln(p’) cho thấy:

  • Dốc NCL = \(\lambda\) (nén bình thường)
  • Dốc CSL = \(\lambda\) (song song với NCL)
  • Khoảng cách CSL-NCL = \(\Gamma\) (tham số trạng thái)
  • Đất cố kết quá mức (OCR > 1) nằm phía trên NCL

5. Ưu Điểm lập trình Julia so với các Ngôn ngữ khác

5.1 Hiệu Suất Tính Toán

Julia là ngôn ngữ được thiết kế đặc biệt cho tính toán khoa học:

Tác Vụ Python (NumPy) MATLAB Julia Cải Thiện
7500 vòng lặp ma trận ~45s ~30s ~2s 15-20× nhanh hơn
Biên dịch trước ✗ ✗ ✓ Không có overhead startup
Năng suất bộ nhớ Cao Cao Thấp Sử dụng 30-50% ít hơn

5.2 Cú Pháp Thân Thiện

Julia kết hợp tính dễ đọc của Python với hiệu suất của C/Fortran:

# Ký hiệu toán học trực tiếp trong code
λ = 0.2              # Thay vì lambda = 0.2
κ = 0.04             # Thay vì kappa = 0.04
ν = 0.15             # Thay vì nu = 0.15

# Phép toán ma trận tự nhiên
D = De - (numerator ./ denominator)  # Broadcasting tự động

5.3 Tích hợp đa nền tảng

# Cùng code chạy trên Windows, macOS, Linux mà không sửa
julia run_mcc.jl  # Linux/macOS
julia.exe run_mcc.jl  # Windows

6. Ứng Dụng & Mở Rộng

6.1 Tham Số Hóa cho Các Loại Đất

Bảng thông số tiêu chuẩn cho các loại đất:

Loại Đất \(M\) \(\lambda\) \(\kappa\) \(N\) \(\nu\)
Sét Bangkok 0.92 0.18 0.035 2.3 0.18
Sét Mexico City 0.95 0.24 0.050 2.8 0.20
Sét Boom 0.88 0.16 0.030 2.1 0.15
Sét cực mịn 1.00 0.25 0.055 3.0 0.22

6.2 Phát triển mở rộng

  1. Mô hình Mặt chảy nâng cao
    • Lade-Duncan yield surface
    • Matsuoka-Nakai model
    • MIT-E3 model
  2. Cơ học đất nâng cao
    • Tải trọng động (Cyclic loading)
    • Ảnh hưởng dị hướng (Anisotropy)
    • Mô hình phụ thuộc ứng suất biến dạng(Stress-strain path dependency)
  3. Phân tích số
    • Phương pháp phần tử hữu hạn (FEM) 2D/3D
    • Phương pháp sai phân hữu hạn (FDM)
    • Mô phỏng động lực học rõ ràng

7. Kết Luận

Bài viết này đã trình bày:

✓ Lý Thuyết MCC: Mặt chảy ellipse, NCL, CSL, quy luật chảy
✓ Triển Khai Julia: Code sạch, hiệu quả, dễ hiểu
✓ Phân Tích Kết Quả: 4 biểu đồ toàn diện cho hiểu biết sâu sắc
✓ Ứng Dụng Thực: Từ lý thuyết đến mô phỏng số

Mô hình MCC là công cụ cơ bản không thể thiếu trong:

  • Thiết kế nền móng
  • Phân tích ổn định huyệt quanh mìn
  • Dự báo chuyển động của công trình
  • Nghiên cứu địa kỹ thuật tiên tiến

Julia, với hiệu suất vượt trội, là lựa chọn hoàn hảo cho các mô phỏng số phức tạp trong cơ học đất hiện đại.

Tài Liệu Tham Khảo

  1. Atkinson, J. H., & Bransby, P. L. (1978). The Mechanics of Soils and Foundations. McGraw-Hill.
  2. Roscoe, K. H., & Burland, J. B. (1968). On the generalized stress-strain behaviour of wet clay. In Engineering Plasticity. Cambridge University Press.
  3. Wood, D. M. (1990). Soil Behaviour and Critical State Soil Mechanics. Cambridge University Press.
  4. Muir Wood, D. (2009). *Soil Behaviour and Critical State

Bài viết này được tạo bằng Julia 1.9.4 với Plots.jl, GR backend.
Cập nhật lần cuối: 2025-12-21 15:30 UTC+7