## import Base: ÷ using Plots using ColorSchemes Plots.scalefontsizes(1.25) ## const R::Float64 = 8.31446261815324 const M_Si::Float64 = 28.0855 const M_C::Float64 = 12.011 const M_Fe::Float64 = 55.845 ## a::Real ÷ b::Real = round(Int, a / b) function wt2mol(wt, M) Σxm = sum(wt ./ M) return wt ./ M ./ Σxm end function mol2wt(x::Vector{Float64}, M) Σm = sum(x .* M) return x .* M ./ Σm end function get_mol2wt(y_Si, y_C) ny = length(y_Si) M = [28.0855, 12.0107, 55.845] wt_Si = zeros(ny) wt_C = zeros(ny) for i in 1:ny wt_Si[i], wt_C[i], _ = mol2wt([y_Si[i], y_C[i], (1 - y_Si[i])], M) end return wt_Si, wt_C end ## function update(D_CC, D_CM, D_MC, D_MM, μ_C1, y_C1, y_Si1, T) for i in eachindex(y_C1) y_C = y_C1[i] y_Si = y_Si1[i] G1 = 0 G2 = 51000 - 21.8T G3 = 77207 - 15.877T G4 = -20510 + 38.7T L1 = -125248 + 41.116T - 142708 * (1 - 2y_Si) + 89907 * (1 - 2y_Si)^2 L2 = 143219.9 + 39.31T - 216320.5 * (1 - 2y_Si) L3 = -34671 L4 = 0 D_L1 = 719256y_Si - 74212 D2_L1 = 719256 D_L2 = 432641 D2_L2 = 0 Ω_Fe = 7.0e-5 * exp(-286000 / (R * T)) / (R * T) Ω_Si = 9.05e-3 * exp(-322465 / (R * T)) / (R * T) Ω_C = 4.529e-7 * exp(-(1 - 2.221e-4 * T) * (-72007 * y_C + 147723 * (1 - y_C)) / (R * T)) / (R * T) dGm_dyC = ( -(1 - y_Si) * G1 - y_Si * G2 + (1 - y_Si) * G3 + y_Si * G4 + R * T * (log(y_C) - log(1 - y_C)) - y_Si * (1 - y_Si) * L1 + y_Si * (1 - y_Si) * L2 + (1 - 2y_C) * (1 - y_Si) * L3 + y_Si * (1 - 2y_C) * L4 ) d2Gm_dyC2 = R * T * (1 / y_C + 1 / (1 - y_C)) - 2 * (1 - y_Si) * L3 - 2 * y_Si * L4 d2Gm_dySi2 = ( R * T * (1 / y_Si + 1 / (1 - y_Si)) - 2 * (1 - y_C) * L1 - 2 * y_C * L2 + 2 * (1 - 2y_Si) * (1 - y_C) * D_L1 + 2 * y_C * (1 - 2y_Si) * D_L2 + y_Si * (1 - y_Si) * (1 - y_C) * D2_L1 + y_C * y_Si * (1 - y_Si) * D2_L2 ) # ySi_1ySi_d2Gm_dySi2 = ( # -2 * (1 - y_C) * L1 - 2 * y_C * L2 # + 2 * (1 - 2y_Si) * (1 - y_C) * D_L1 + 2 * y_C * (1 - 2y_Si) * D_L2 # + y_Si * (1 - y_Si) * (1 - y_C) * D2_L1 + y_C * y_Si * (1 - y_Si) * D2_L2 # ) * y_Si * (1 - y_Si) + R * T d2Gm_dySidyC = G1 - G2 - G3 + G4 + (1 - 2y_Si) * (L2 - L1) + (1 - 2y_C) * (L4 - L3) + y_Si * (1 - y_Si) * (D_L2 - D_L1) μ_C = dGm_dyC dμC_dyC = d2Gm_dyC2 dμC_dySi = d2Gm_dySidyC dμSi_dyC = μ_C + (1 - y_Si) * d2Gm_dySidyC dμSi_dySi = (1 - y_Si) * d2Gm_dySi2 dμFe_dyC = μ_C - y_Si * d2Gm_dySidyC dμFe_dySi = -y_Si * d2Gm_dySi2 # ySi_dμSi_dySi = ySi_1ySi_d2Gm_dySi2 # _1ySi_dμFe_dySi = -ySi_1ySi_d2Gm_dySi2 μ_C1[i] = μ_C D_CC[i] = y_C * (1 - y_C) * Ω_C * dμC_dyC D_CM[i] = y_C * (1 - y_C) * Ω_C * dμC_dySi D_MC[i] = (1 - y_Si) * y_Si * Ω_Si * dμSi_dyC - y_Si * (1 - y_Si) * Ω_Fe * dμFe_dyC D_MM[i] = (1 - y_Si) * y_Si * Ω_Si * dμSi_dySi - y_Si * (1 - y_Si) * Ω_Fe * dμFe_dySi # D_MM[i] = (1 - y_Si) * Ω_Si * ySi_dμSi_dySi - y_Si * Ω_Fe * _1ySi_dμFe_dySi end end ## function solve(wt_L, wt_R, T::Float64, L, t_max, dt_save::Float64, dx::Float64, dt::Float64) # Conversion M = [M_Si, M_C, M_Fe] x_Si_L, x_C_L, x_Fe_L = wt2mol(wt_L, M) x_Si_R, x_C_R, x_Fe_R = wt2mol(wt_R, M) # Arrays length nx = (L ÷ dx) + 1 nt_save = (t_max ÷ dt_save) + 1 nt = (t_max ÷ dt) + 1 ny = 2nx - 1 t_range = [0.0] x_range = range(-L, L, ny) # Allocate arrays y_C1 = zeros(ny) y_C2 = zeros(ny) y_Si1 = zeros(ny) y_Si2 = zeros(ny) D_CC = zeros(ny) D_CM = zeros(ny) D_MC = zeros(ny) D_MM = zeros(ny) μ_C1 = zeros(ny) # Initial value # Left y_C1[1:nx] .= x_C_L / (x_Si_L + x_Fe_L) y_Si1[1:nx] .= x_Si_L / (x_Si_L + x_Fe_L) # Right y_C1[nx+1:end] .= x_C_R / (x_Si_R + x_Fe_R) y_Si1[nx+1:end] .= x_Si_R / (x_Si_R + x_Fe_R) # Initial update(D_CC, D_CM, D_MC, D_MM, μ_C1, y_C1, y_Si1, T) y_C = [copy(y_C1)] y_Si = [copy(y_Si1)] μ_C = [copy(μ_C1)] for j in 2:nt update(D_CC, D_CM, D_MC, D_MM, μ_C1, y_C1, y_Si1, T) for i in 2:(ny-1) y_C2[i] = ( y_C1[i] + dt / dx^2 * (sqrt(D_CC[i+1] * D_CC[i]) * (y_C1[i+1] - y_C1[i]) - sqrt(D_CC[i] * D_CC[i-1]) * (y_C1[i] - y_C1[i-1])) + dt / dx^2 * (sqrt(D_CM[i+1] * D_CM[i]) * (y_Si1[i+1] - y_Si1[i]) - sqrt(D_CM[i] * D_CM[i-1]) * (y_Si1[i] - y_Si1[i-1])) ) y_Si2[i] = ( y_Si1[i] + dt / dx^2 * (sqrt(D_MC[i+1] * D_MC[i]) * (y_C1[i+1] - y_C1[i]) - sqrt(D_MC[i] * D_MC[i-1]) * (y_C1[i] - y_C1[i-1])) + dt / dx^2 * (sqrt(D_MM[i+1] * D_MM[i]) * (y_Si1[i+1] - y_Si1[i]) - sqrt(D_MM[i] * D_MM[i-1]) * (y_Si1[i] - y_Si1[i-1])) ) end y_C2[1] = ( y_C1[1] + dt / dx^2 * (sqrt(D_CC[1+1] * D_CC[1]) * (y_C1[1+1] - y_C1[1]) - sqrt(D_CC[1] * D_CC[1+1]) * (y_C1[1] - y_C1[1+1])) + dt / dx^2 * (sqrt(D_CM[1+1] * D_CM[1]) * (y_Si1[1+1] - y_Si1[1]) - sqrt(D_CM[1] * D_CM[1+1]) * (y_Si1[1] - y_Si1[1+1])) ) y_Si2[1] = ( y_Si1[1] + dt / dx^2 * (sqrt(D_MC[1+1] * D_MC[1]) * (y_C1[1+1] - y_C1[1]) - sqrt(D_MC[1] * D_MC[1+1]) * (y_C1[1] - y_C1[1+1])) + dt / dx^2 * (sqrt(D_MM[1+1] * D_MM[1]) * (y_Si1[1+1] - y_Si1[1]) - sqrt(D_MM[1] * D_MM[1+1]) * (y_Si1[1] - y_Si1[1+1])) ) y_C2[end] = ( y_C1[end] + dt / dx^2 * (sqrt(D_CC[end-1] * D_CC[end]) * (y_C1[end-1] - y_C1[end]) - sqrt(D_CC[end] * D_CC[end-1]) * (y_C1[end] - y_C1[end-1])) + dt / dx^2 * (sqrt(D_CM[end-1] * D_CM[end]) * (y_Si1[end-1] - y_Si1[end]) - sqrt(D_CM[end] * D_CM[end-1]) * (y_Si1[end] - y_Si1[end-1])) ) y_Si2[end] = ( y_Si1[end] + dt / dx^2 * (sqrt(D_MC[end-1] * D_MC[end]) * (y_C1[end-1] - y_C1[end]) - sqrt(D_MC[end] * D_MC[end-1]) * (y_C1[end] - y_C1[end-1])) + dt / dx^2 * (sqrt(D_MM[end-1] * D_MM[end]) * (y_Si1[end-1] - y_Si1[end]) - sqrt(D_MM[end] * D_MM[end-1]) * (y_Si1[end] - y_Si1[end-1])) ) # y_C2[1] = y_C1[1] # y_Si2[1] = y_Si1[1] # y_C2[end] = y_C1[end] # y_Si2[end] = y_Si1[end] if ((j - 1) % ((nt - 1) ÷ (nt_save - 1)) == 0) || j == nt update(D_CC, D_CM, D_MC, D_MM, μ_C1, y_C2, y_Si2, T) @info "λ = $(dt/dx^2*maximum([maximum(D_CC), maximum(D_CM), maximum(D_MC), maximum(D_MM)]))" push!(μ_C, copy(μ_C1)) push!(t_range, j * dt) push!(y_C, copy(y_C2)) push!(y_Si, copy(y_Si2)) end for i in 1:ny y_C1[i] = y_C2[i] y_Si1[i] = y_Si2[i] end end return y_C, y_Si, t_range, x_range, μ_C end ## # Left wt_Si1 = 3.8 / 100 wt_C1 = 0.478 / 100 wt_Fe1 = 1 - wt_Si1 - wt_C1 # Right wt_Si2 = 0.000001 / 100 # wt_Si2 = 0 / 100 wt_C2 = 0.441 / 100 wt_Fe2 = 1 - wt_Si2 - wt_C2 # Parameters T = 1323.0 L = 25e-3 t_max = 13 * 3600.0 * 24 dt_save = 13 * 3600.0 * 24 dx = 0.1e-3 dt = 30.0 # @time y_C, y_Si, t_range, x_range, μ_C = solve([wt_Si1, wt_C1, wt_Fe1], [wt_Si2, wt_C2, wt_Fe2], T, L, t_max, dt_save, dx, dt) nothing ## ref_p = [-19.5 -17.5 -15.5 -13.5 -11.0 -9.0 -7.0 -4.6 -3.6 -2.7 -1.5 -0.5 0.5 1.5 2.5 3.5 5.5 7.0 11.0 15.1 19.4 23.3; 0.472 0.45 0.455 0.432 0.432 0.418 0.401 0.379 0.363 0.359 0.333 0.422 0.574 0.565 0.55 0.555 0.534 0.518 0.49 0.472 0.458 0.455] p = plot(ylabel=raw"wt% Carbon", xlabel=raw"Distance $(mm)$") x = x_range * 1e3 idx = findall(x -> (-30 <= x <= 30), x) _t = round(Int, t_range[end] / (3600 * 24)) # wt_Si, wt_C = get_mol2wt(y_Si[end], y_C[end]) y = wt_C .* 100 plot!(x[idx], y[idx], label="Day $_t") scatter!(ref_p[1, :], ref_p[2, :], label="ref") plot!(ylims=(0.25, 0.65), yticks=(0.3:0.1:0.6, string.(0.3:0.1:0.6))) # savefig("output/compare.png") ## p = plot(ylabel="%wt", xlabel=raw"Distance $(mm)$") plot!(title="Initial Condition") wt_Si, wt_C = get_mol2wt(y_Si[1], y_C[1]) plot!(x_range * 1e3, wt_Si * 100, label=raw"$Si$") plot!(x_range * 1e3, wt_C * 100, label=raw"$C$") # savefig("output/initial.png") ## t_max = 100 * 3600.0 * 24 dt_save = 1 * 3600.0 * 24 y_C, y_Si, t_range, x_range, μ_C = solve([wt_Si1, wt_C1, wt_Fe1], [wt_Si2, wt_C2, wt_Fe2], T, L, t_max, dt_save, dx, dt) ## jdx = [1, 5, 10, 20, 50, 100] .+ 1 colori = [get(ColorSchemes.viridis, j) for j in range(0, 1, length(jdx) + 1)] nothing ## p = plot(ylabel=raw"wt% Carbon", xlabel=raw"Distance $(mm)$") len_t = length(t_range) for (i, j) in enumerate(jdx) # x = x_range * 1e3 idx = findall(x -> (-30 <= x <= 30), x) _t = round(Int, t_range[j] / (3600 * 24)) # wt_Si, wt_C = get_mol2wt(y_Si[j], y_C[j]) y = wt_C .* 100 plot!(x[idx], y[idx], label="Day $_t", color=colori[i], line=1.5) end plot!(ylims=(0.25, 0.65), yticks=(0.3:0.1:0.6, string.(0.3:0.1:0.6))) plot!(legend=:topleft) # savefig("output/results.png") ## p = plot(ylabel=raw"wt% $Si$", xlabel=raw"Distance $(mm)$") len_t = length(t_range) for (i, j) in enumerate(jdx) wt_Si, wt_C = get_mol2wt(y_Si[j], y_C[j]) _t = round(Int, t_range[j] / (3600 * 24)) plot!(x_range * 1e3, wt_Si, label="Day $_t", color=colori[i]) end p # savefig("output/Si.png") ## p = plot(ylabel="Chemical Potential", xlabel=raw"Distance $(mm)$") len_t = length(t_range) for (i, j) in enumerate([1; jdx]) _t = round(Int, t_range[j] / (3600 * 24)) plot!(x_range * 1e3, μ_C[j], label="Day $(_t)", color=colori[i], line=1.5) end p # savefig("output/Chemical_Potential.png")