k (interfacial resistance) is now parameter. Add plotting function.

eff_Drop.jl

@@ -4,6 +4,7 @@ using Plots
 using DelimitedFiles
 using Interpolations
 using LinearAlgebra
+using Revise
 
 """
     analytical(x, x0, Nb, Tf, Di, Do, P)
@@ -99,7 +100,7 @@ end
 
 """
     make_A(geometry,B,α_i, α_o, Δxi, Δxo, x, D_i, D_o, Nb, P, 
-        BEt)
+        BEt, k)
 
     Create matrix A in u = A  b.
 
@@ -117,9 +118,10 @@ end
     - `D_o::Float`: diffusion coefficient outside.
     - `N::Int`: size of matrix A.
     - `Nb::Int`: index of left ghost point.
-    - `P::Float`: Partition coefficient."""
+    - `P::Float`: Partition coefficient.
+    - `k::Float`: Interfacial resistance"""
 function make_A(geometry,B,α_i, α_o, Δxi, Δxo, x, D_i, D_o, N, 
-        Nb, P, BEt)
+        Nb, P, BEt, k)
     A = zeros(N,N)
 
     if geometry=="3D"
@@ -142,12 +144,13 @@ function make_A(geometry,B,α_i, α_o, Δxi, Δxo, x, D_i, D_o, N,
         # A[2,1:3]=[-2*α_i,1+4*α_i,-2*α_i]
         A[2,1:3]=[-3*α_i,1+6*α_i,-3*α_i]
     end
-
     # Boundary Conditions:
     if B=="N"
-        A[Nb, Nb-1:Nb+2] = [1, 0, 0, -P]
-        # k = 0.01;
-        # A[Nb, Nb-2:Nb+2] =   [-D_i, k, D_i, 0,    -k*P]
+        if k=="none"
+            A[Nb, Nb-1:Nb+2] = [1, 0, 0, -P]
+        else
+            A[Nb, Nb-2:Nb+2] =   [-D_i, k, D_i, 0,    -k*P]
+        end
         A[Nb+1, Nb-2:Nb+3] = [D_i/2 , 0,-D_i/2,-D_o/2,0,   D_o/2]
         A[end, end-2:end] = [-α_o,0, α_o]
     elseif B=="E"
@@ -202,7 +205,7 @@ end
     - `perc::Float`: percentage of recovery
     """
 function solve_D(geometry,IC,AG,B,u0E,P,D_i,D_o,Δx,Δt,
-        L,R,BP,x0,T,perc,BE,tE)
+        L,R,BP,x0,T,perc,BE,tE, k)
     # https://www.uni-muenster.de/imperia/md/content/physik_tp/
     # lectures/ws2016-2017/num_methods_i/heat.pdf
     #Defining x, Nb,  α_i,  α_o
@@ -284,10 +287,10 @@ function solve_D(geometry,IC,AG,B,u0E,P,D_i,D_o,Δx,Δt,
         Δxi, Δxo, x, Nb)
     if B=="N"
         A = make_A(geometry,B,α_i(D_i,Δt,Δxi), α_o(D_o,Δt,Δxo), 
-            Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, 1)
+            Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, 1, k)
     elseif B=="E"
         A = make_A(geometry,B,α_i(D_i,Δt,Δxi), α_o(D_o,Δt,Δxo), 
-            Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, BE_t(Δt))
+            Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, BE_t(Δt), k)
     end
 
     # compute u_tot_INSIDE at equilibrium
@@ -343,7 +346,7 @@ function solve_D(geometry,IC,AG,B,u0E,P,D_i,D_o,Δx,Δt,
         if B=="E"
             A=make_A(geometry,B,α_i(D_i,Δt,Δxi), α_i(D_i,Δt,Δxo), 
                 Δxi, Δxo, x, D_i, D_o, length(x), Nb, P, 
-                BE_t(round(Δt*(i-1),digits=length(string(Δt))-2)))
+                BE_t(round(Δt*(i-1),digits=length(string(Δt))-2)), k)
         end
 
 
@@ -363,7 +366,7 @@ function solve_D(geometry,IC,AG,B,u0E,P,D_i,D_o,Δx,Δt,
         #     if B=="E"
         #         A=make_A(geometry,B,α_i(D_i,Δt,Δxi), 
         #           α_i(D_i,Δt,Δxo), Δxi, Δxo, x, D_i, D_o, 
-        #           length(x), Nb, P, BE_t(Δt*(i)))
+        #           length(x), Nb, P, BE_t(Δt*(i)), k)
         #     end
         #     u_next=A\b
         #     # println(u[Id-1])
@@ -420,4 +423,20 @@ function solve_D(geometry,IC,AG,B,u0E,P,D_i,D_o,Δx,Δt,
 
     return A, u_mat, x, b_mat, Nb, t_perc, α_i(D_i,Δt,Δxi), 
             α_o(D_o,Δt,Δxo), Mm, t_vec, Ni, No, m_mat, u0
-end
\ No newline at end of file
+end
+
+function plot_model(x, u, Nb, T, Δt)
+    u = u[:, vcat(1:Nb-1, Nb+2:size(u, 2))]
+    x = vcat(collect(x[1:Nb-1]), collect(x[Nb+2:end]))
+    plt=Plots.plot(xlabel= "x", ylabel = "u", 
+        foreground_color_legend = nothing, palette=:seaborn_colorblind)
+
+    plot!(plt,x, u[1,:], label="Numerics", lw=3, color=1)
+    for i in 101:500:(Int(T/Δt))
+        plot!(plt, x, u[i, :], label=false, color=1,
+            lw=3, xlims=(0, 2))
+    end
+    display(plt)
+    # # Plots.savefig(plt,"/Users/mestroni/Figures/3D_analytical.svg")
+    # Plots.plot!(u[:, Nb-1]./u[:, Nb])
+end