salt implemented

MATLAB/Model/p_v_sat_salt.m

+%-------------------------------------------------------------------------------
+function p_v = p_v_sat_salt(t, a)
+%-------------------------------------------------------------------------------
+% t      temperature in °C
+% a      salt concentration in kg/kg
+% p_v    saturation vapor pressure in Pa
+% Correlations from Bridgeman and Aldrich, 1964 (NIST)
+%-------------------------------------------------------------------------------
+
+t = t + 273.15; % to convert to K
+
+i = find(t>334 & t<=363.15);
+j = find(t>304 & t<=334);
+k = find(t>273 & t<=304);
+
+p_v = zeros(size(t));
+
+p_v(i) = 1E+5 * 10.^(5.0768  - 1659.793./(t(i)-45.854));
+p_v(j) = 1E+5 * 10.^(5.20389 - 1733.926./(t(j)-39.485));
+p_v(k) = 1E+5 * 10.^(5.40221 - 1838.675./(t(k)-31.737));
+
+% correction for salt
+p_v = p_v .* ( 1 - 0.5.*a - 10.*a.*a ) .* ( 1 - a );
+
+p_v(p_v == 0) = NaN;
+
+end
\ No newline at end of file

MATLAB/Model/sim_mem.m

 %==========================================================================
 %
-% This program demonstrate the use of two independent fluids in two
-% independent domains (one in water, the other in air), for simulation of 
-% membrane for Kerstin.
-%
-% It is important to note that all variables in each of the domains (water 
-% and air) are discretized independently, but temperature is coupled 
-% through the system matrix.
-%
-% A nice thing about this implementation is that it didn't need
-% modifications in any of the supporting functions (for example in "Core")
-%
-% An interesting thing, for the future, that it can serve as a basis for
-% implemetation of two-fluid model (but it has nothing to do with Kerstin's
-% thesis).
-%
 %==========================================================================
 clear
 clc;
@@ -54,7 +39,8 @@ for c=H2O:FIL
   dy(c)=ly(c)/ny(c);   nym(c)=ny(c)-1;
 end  
 
-u_max = 0.12;
+u_max  = 0.12;
+a_salt = 90; % g/l
 
 % set physical properties
 [prop(H2O).rho, ...
@@ -71,7 +57,7 @@ u_max = 0.12;
  prop(FIL).lambda] = water_properties(1:nx(FIL), 1:ny(FIL));
 
 prop(AIR).diff(1:nx(AIR),1:ny(AIR)) = 4.0E-4;    % diffusivity of vapor in air
-prop(H2O).diff(1:nx(H2O),1:ny(H2O)) = 0.0;
+prop(H2O).diff(1:nx(H2O),1:ny(H2O)) = 1.99E-09;
 
 % prop(FIL).lambda(:,ny(FIL)) = 1.0; 
 % prop(AIR).lambda(:,1)       = 0.5;
@@ -86,10 +72,10 @@ pi= 3.142;
 
 % Membrane properties
 mem.d      = 166E-6;  % m thickness
-mem.lambda = 0.2;    % W/mK thermal conductivity
-mem.eps    = 0.687;    % porosity
-mem.tau    = 2;      % tortuosity
-mem.r      = 0.1E-6; % m pore radius
+mem.lambda = 0.2;     % W/mK thermal conductivity
+mem.eps    = 0.687;   % porosity
+mem.tau    = 2;    % tortuosity
+mem.r      = 0.1E-6;  % m pore radius
 
 % membrane values
 mem.p  = zeros(nx(AIR), 1);
@@ -98,7 +84,6 @@ mem.a  = zeros(nx(AIR), 1);
 mem.pv = zeros(nx(AIR), 1);
 mem.j  = zeros(nx(AIR), 1);
 
-
 % create unknowns; names, positions and sizes
 for c=H2O:FIL 
   u(c) = create_unknown('u-velocity',   'x', [nxm(c) ny(c) ], 'd');
@@ -124,7 +109,10 @@ t(AIR).val(1:nx(AIR),  1:ny(AIR)) = 40;
 t(AIR).val(1:nx(AIR),  1:ny(AIR)/3)           = 20;
 t(AIR).val(1:nx(AIR),  ny(AIR)/3:2*ny(AIR)/3) = 40;
 t(AIR).val(1:nx(AIR),  2*ny(AIR)/3:ny(AIR))   = 70;
+
 u(H2O).val(1:nxm(H2O), 1:ny(H2O)) = repmat(4.*0.1.*((1:ny(H2O))).*(ny(H2O)+1-(1:ny(H2O)))/(ny(H2O)+1)^2,nxm(H2O),1);
+
+a(H2O).val(1:nx(H2O),  1:ny(H2O)) = a_salt./prop(H2O).rho(1:nx(H2O), 1:ny(H2O));
 M.val = M_AIR;
 a(AIR).val(1:nx(AIR), 1:ny(AIR)) = p_v_sat(t(AIR).val)*1E-5*M_H2O./M_AIR;
 M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
@@ -136,22 +124,23 @@ M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
 % temperature
 u(H2O).bnd(W).type(1,1:ny(H2O)) = 'd'; u(H2O).bnd(W).val(1,1:ny(H2O)) = 4.*0.1.*((1:ny(H2O))).*(ny(H2O)+1-(1:ny(H2O)))/(ny(H2O)+1)^2;
 u(H2O).bnd(E).type(1,1:ny(H2O)) = 'o'; u(H2O).bnd(E).val(1,1:ny(H2O)) = 4.*0.1.*((1:ny(H2O))).*(ny(H2O)+1-(1:ny(H2O)))/(ny(H2O)+1)^2;
+v(H2O).bnd(S).type(1:nx(H2O),1) = 'd'; v(H2O).bnd(S).val(1:nx(H2O),1) = - mem.j;
 
 u(AIR).bnd(W).type(1,1:ny(AIR)) = 'd'; u(AIR).bnd(W).val(1,1:ny(AIR)) =   0.0;
 u(AIR).bnd(E).type(1,1:ny(AIR)) = 'd'; u(AIR).bnd(E).val(1,1:ny(AIR)) =  +0.0;
 
 t(H2O).bnd(W).type(1,1:ny(H2O)) = 'd'; t(H2O).bnd(W).val(1,1:ny(H2O)) = 80.0;
-t(H2O).bnd(E).type(1,1:ny(H2O)) = 'n'; t(H2O).bnd(E).val(1,1:ny(H2O)) = 60.0;
-t(H2O).bnd(S).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(S).val(1:nx(H2O),1) = 60.0;
-t(H2O).bnd(N).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(N).val(1:nx(H2O),1) = 60.0;
+t(H2O).bnd(E).type(1,1:ny(H2O)) = 'n';
+t(H2O).bnd(S).type(1:nx(H2O),1) = 'd'; t(H2O).bnd(S).val(1:nx(H2O),1) = 80.0;
+t(H2O).bnd(N).type(1:nx(H2O),1) = 'n';
 t(H2O) = adjust_n_boundaries(t(H2O));
 
 t_int_mem = t(H2O).bnd(S).val;
 
-t(AIR).bnd(W).type(1,1:ny(AIR)) = 'n'; t(AIR).bnd(W).val(1,1:ny(AIR)) = 70.0;
-t(AIR).bnd(E).type(1,1:ny(AIR)) = 'n'; t(AIR).bnd(E).val(1,1:ny(AIR)) = 70.0;
+t(AIR).bnd(W).type(1,1:ny(AIR)) = 'n';
+t(AIR).bnd(E).type(1,1:ny(AIR)) = 'n';
 t(AIR).bnd(S).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(S).val(1:nx(AIR),1) = 20.0;
-t(AIR).bnd(N).type(1:nx(AIR),1) = 'n'; t(AIR).bnd(N).val(1:nx(AIR),1) = 70.0;
+t(AIR).bnd(N).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(N).val(1:nx(AIR),1) = 70.0;
 t(AIR) = adjust_n_boundaries(t(AIR));
 
 t(FIL).bnd(S).type(1:nx(FIL),1) = 'd'; t(FIL).bnd(S).val(1:nx(FIL),1) = 20.0;
@@ -164,6 +153,9 @@ p_v.bnd(N).type(1:nx(AIR),1)    = 'd'; p_v.bnd(N).val(1:nx(AIR),1)    = p_v_sat(
 M.bnd(S).type(1:nx(AIR),1)      = 'd'; M.bnd(S).val(1:nx(AIR),1)      = M.val(:,1);
 M.bnd(N).type(1:nx(AIR),1)      = 'd'; M.bnd(N).val(1:nx(AIR),1)      = M.val(:,ny(AIR));
 
+a(H2O).bnd(W).type(1,1:ny(H2O)) = 'd'; a(H2O).bnd(W).val(1,1:ny(H2O)) = a_salt./prop(H2O).rho(1,:);
+a(H2O).bnd(E).type(1,1:ny(H2O)) = 'd'; a(H2O).bnd(E).val(1,1:ny(H2O)) = a(H2O).val(nx(H2O),1:ny(H2O));
+
 a(AIR).bnd(S).type(1:nx(AIR),1) = 'n';
 a(AIR).bnd(N).type(1:nx(AIR),1) = 'd';
 a(AIR).bnd(N).val = p_v.bnd(N).val./(p(H2O).tot(:,1) + 1E5)*M_H2O./M.bnd(N).val; % p_v/p_tot * M_H2O/M_mix
@@ -175,9 +167,9 @@ obst = [];
 % load('workspace_sim_mem_20000_60');
 
 % time-stepping parameters
-dt  =    0.0001;  % time step
-ndt =    200000;  % number of time steps
-% n_start = 20001;
+dt  =    0.001;  % time step
+ndt =    10;  % number of time steps
+n_start = 1;
 
 %-----------
 %
@@ -234,8 +226,8 @@ for n=n_start:ndt
   mem.M  = (M.bnd(N).val      + M.val(:,ny(AIR)))     /2.0;
   
   % Diffusion Coefficients
-  C_K = 2*mem.eps*mem.r/(3*mem.tau*mem.d*dy(AIR)).*(8*M_H2O./(mem.t*R*pi));
-  C_M = mem.eps*mem.p.*prop(AIR).diff(:,ny(AIR))./(mem.tau*R*mem.t.*(mem.p-mem.pv)*dy(AIR));
+  C_K = 2*mem.eps*mem.r/(3*mem.tau*mem.d).*(8*M_H2O./(mem.t*R*pi)).^0.5;
+  C_M = mem.eps*mem.p.*prop(AIR).diff(:,ny(AIR))./(mem.tau*R*mem.t.*(mem.p-mem.pv));
   C_T = 1.0./(1.0./C_K + 1.0./C_M);
   
   % jump condition at membrane
@@ -252,7 +244,7 @@ for n=n_start:ndt
         .* t(AIR).val(:,ny(AIR)) );
       
   for i=1:nx(H2O)
-    jump_cond_mem = @(t) lhs_lin_mem(i) * t + lhs_fun_mem(i) * p_v_sat(t) - rhs_mem(i);
+    jump_cond_mem = @(t) lhs_lin_mem(i) * t + lhs_fun_mem(i) * p_v_sat_salt(t, a(H2O).val(i,1)) - rhs_mem(i);
     t_int_mem(i) = fzero(jump_cond_mem,t(H2O).val(i,1));
   end
 
@@ -261,7 +253,7 @@ for n=n_start:ndt
   % AIR N bnd values 
   t(H2O).bnd(S).val = t_int_mem; 
   t(AIR).bnd(N).val = t_int_mem; 
-  p_v.bnd(N).val = p_v_sat(t_int_mem);
+  p_v.bnd(N).val = p_v_sat_salt(t_int_mem, a(H2O).val(:,1));
   
   M.bnd(N).val = M_AIR*(1-a(AIR).bnd(N).val) + M_H2O*a(AIR).bnd(N).val;
   a(AIR).bnd(N).val = p_v.bnd(N).val./(p(H2O).tot(:,1) + 1E+5)*M_H2O./M.bnd(N).val; % p_v/p_tot * M_H2O/M_mix
@@ -271,6 +263,8 @@ for n=n_start:ndt
   mem.flowdir = zeros(nx(AIR),1);
   mem.flowdir = 1.0*(mem.j >= 0.0);
   
+  v(H2O).bnd(S).val(1:nx(H2O),1) = - mem.j;
+  
   %------------------------
   % concentration
   %------------------------
@@ -286,14 +280,17 @@ for n=n_start:ndt
     
     % set v(AIR).bnd(N) temporarily to zero to avoid convection from
     % membrane into domain; it is accounted for by j_a
-    v_temp = v(AIR).bnd(N).val;
+    v_temp_AIR = v(AIR).bnd(N).val;
+    v_temp_H2O = v(H2O).bnd(S).val;
     v(AIR).bnd(N).val(1:nx(AIR),1) = 0.0; 
+    v(H2O).bnd(S).val(1:nx(H2O),1) = 0.0;
     
     % advection terms for momentum                            
     c_a = advection(prop(c).rho, a(c), u(c), v(c), dx(c), dy(c), dt, 'superbee');
     
     % set v(AIR).bnd(N) back to old value
-    v(AIR).bnd(N).val = v_temp;
+    v(AIR).bnd(N).val = v_temp_AIR;
+    v(H2O).bnd(S).val = v_temp_H2O;
     
     % innertial term for concentration
     i_a = reshape( (a(c).val .*prop(c).rho ./dt *(dx(c)*dy(c))), nx(c)*ny(c),1);   
@@ -327,6 +324,7 @@ for n=n_start:ndt
     % solve concentration
     a(c).val = reshape( A_a\(f_a), nx(c), ny(c) );
     a(c) = adjust_n_boundaries(a(c));
+    a(H2O).bnd(E).val(1,1:ny(H2O)) = a(H2O).val(nx(H2O),1:ny(H2O));
                                      
   end
   
@@ -351,12 +349,7 @@ for n=n_start:ndt
                                      (prop(FIL).rho.*prop(FIL).cap)./dt, ... % innertial term
                                      prop(FIL).lambda,                   ... % diffusive coefficent
                                      dx(FIL), dy(FIL),                   ... % dimensions
-                                     obst);                              ... % obstacles
- 
-  % combine two matrices
-  A_t = [ [ A_t_air ; zeros(nx(H2O)*ny(H2O), nx(AIR)*ny(AIR)) ], ...
-          [ zeros(nx(AIR)*ny(AIR), nx(H2O)*ny(H2O)) ; A_t_h2o ] ];
-  b_t = [ b_t_air; b_t_h2o ];    
+                                     obst);                              ... % obstacles  
   
   for i=1:nx(H2O)
     %>>>>>>>>>>>>>>>>>>>>>>
@@ -367,22 +360,19 @@ for n=n_start:ndt
     j = ny(AIR);
     [c_air, w_air, e_air, s_air, n_air] = loc(i,j,nx(AIR),ny(AIR));
     
+    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,ny(AIR))*mem.eps;
+    cy = -prop(AIR).lambda(i,ny(AIR)) * 2*dx(AIR)/dy(AIR) ...
+      + dx(AIR)/(mem.d/lambda_mem + dy(AIR)/(2*prop(AIR).lambda(i,ny(AIR))));
+    
+    b_t_air(c_air) = b_t_air(c_air) + cy * t(AIR).bnd(N).val(i,1) ;
+    A_t_air(c_air, c_air) = A_t_air(c_air, c_air) + cy ...
+      - h_d*mem.j(i).*(1-mem.flowdir(i))/t(AIR).val(i,j);
+    
     % number of the unknown in the water domain (on top of air)
     j = 1;
     [c_h2o, w_h2o, e_h2o, s_h2o, n_h2o] = loc(i,j,nx(H2O),ny(H2O));
-    c_h2o = c_h2o + nx(AIR)*ny(AIR);
-    s_h2o = s_h2o + nx(AIR)*ny(AIR);
-    
-    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,j)*mem.eps;
-    cy = dx(AIR)/(dy(H2O)/(2*prop(H2O).lambda(i,1)) + mem.d/lambda_mem ...
-      + dy(AIR)/(2*prop(AIR).lambda(i,ny(AIR))));
-    
-    A_t(c_air, n_air) = -cy;
-    A_t(c_air, c_air) = A_t(c_air, c_air) + cy ...
-      - h_d*mem.j(i).*(1-mem.flowdir(i))/t(AIR).val(i,j);
 
-    A_t(c_h2o, s_h2o) = -cy;
-    A_t(c_h2o, c_h2o) = A_t(c_h2o, c_h2o) + cy ...
+    A_t_h2o(c_h2o, c_h2o) = A_t_h2o(c_h2o, c_h2o) ...
       + h_d*mem.j(i).*mem.flowdir(i)/t(H2O).val(i,j);
     
     %<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -401,26 +391,26 @@ for n=n_start:ndt
   c_t_air = advection(prop(AIR).rho.*prop(AIR).cap, t(AIR), u(AIR), v(AIR), dx(AIR), dy(AIR), dt, 'superbee');
   c_t_h2o = advection(prop(H2O).rho.*prop(H2O).cap, t(H2O), u(H2O), v(H2O), dx(H2O), dy(H2O), dt, 'superbee');
   c_t_fil = advection(prop(FIL).rho.*prop(FIL).cap, t(FIL), u(FIL), v(FIL), dx(FIL), dy(FIL), dt, 'superbee');
-  c_t = [ c_t_air; c_t_h2o ];
     
   % innertial term for enthalpy; form separatelly then combine
   i_t_air = reshape( (t(AIR).val.*prop(AIR).rho.*prop(AIR).cap/dt*(dx(AIR)*dy(AIR))), nx(AIR)*ny(AIR),1);
   i_t_h2o = reshape( (t(H2O).val.*prop(H2O).rho.*prop(H2O).cap/dt*(dx(H2O)*dy(H2O))), nx(H2O)*ny(H2O),1);
   i_t_fil = reshape( (t(FIL).val.*prop(FIL).rho.*prop(FIL).cap/dt*(dx(FIL)*dy(FIL))), nx(FIL)*ny(FIL),1);
-  i_t = [ i_t_air; i_t_h2o ];
   
   % the entire source term
-  f_t = b_t - c_t + i_t;
+  f_t_air = b_t_air - c_t_air + i_t_air;
+  f_t_h2o = b_t_h2o - c_t_h2o + i_t_h2o;
   f_t_fil = b_t_fil - c_t_fil + i_t_fil;
   
   % solve temperature
-  t_all      = A_t\(f_t);
+  t(AIR).val = A_t_air\(f_t_air);
+  t(H2O).val = A_t_h2o\(f_t_h2o);
   t(FIL).val = A_t_fil\(f_t_fil);
   
   % split the temperature in air and water domain
-  t(AIR).val = reshape(t_all(                1 : nx(AIR)*ny(AIR)),                 nx(AIR), ny(AIR));
-  t(H2O).val = reshape(t_all(nx(AIR)*ny(AIR)+1 : nx(AIR)*ny(AIR)+nx(H2O)*ny(H2O)), nx(H2O), ny(H2O));
-  t(FIL).val = reshape(t(FIL).val,                                                 nx(FIL), ny(FIL));
+  t(AIR).val = reshape(t(AIR).val, nx(AIR), ny(AIR));
+  t(H2O).val = reshape(t(H2O).val, nx(H2O), ny(H2O));
+  t(FIL).val = reshape(t(FIL).val, nx(FIL), ny(FIL));
   
   t(AIR) = adjust_n_boundaries(t(AIR));
   t(H2O) = adjust_n_boundaries(t(H2O));
@@ -457,6 +447,10 @@ for n=n_start:ndt
     p_st_v = reshape( dif_y(p(c).tot) / dy(c) * (dx(c)*dy(c)), nx(c)*nym(c), 1 );
   
     g_v = reshape( (avg_y(prop(c).rho)*(dx(c)*dy(c))), nx(c)*nym(c),1);
+    
+    b_v = reshape( b_v , nx(c), nym(c));
+    b_v(:,1) = - mem.j;
+    b_v = reshape( b_v , nx(c)*nym(c),1);
   
     % full force terms for momentum equations
     f_u = b_u - c_u + i_u - p_st_u;
@@ -552,7 +546,7 @@ for n=n_start:ndt
     % temperature
     subplot(3,2,1);
       [C,h] = contourf(x, y, [t(FIL).val'; t(AIR).val'; t(H2O).val'], ...
-                       linspace(19.99,60.01,10)); 
+                       linspace(19.99,80.01,10)); 
       axis equal
       title(sprintf('Temperature at time step %d out of %d\n', n, ndt));
     % velocity magnitude
@@ -588,6 +582,10 @@ for n=n_start:ndt
     drawnow;
   end
   
+  if mod(n,500) == 0
+    save('workspace_sim_mem.mat')
+  end
+  
 end
 
 % % plot in tecplot format

MATLAB/Model/sim_mem_wo_jump_mem.m

@@ -340,7 +340,7 @@ for n=5596:ndt
     c_h2o = c_h2o + nx(AIR)*ny(AIR);
     s_h2o = s_h2o + nx(AIR)*ny(AIR);
     
-    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,j)*mem.eps;
+    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,ny(AIR))*mem.eps;
     cy = dx(AIR)/(dy(H2O)/(2*prop(H2O).lambda(i,1)) + mem.d/lambda_mem ...
       + dy(AIR)/(2*prop(AIR).lambda(i,ny(AIR))));
     

MATLAB/Model/sim_mem.asv → MATLAB/Model/sim_mem_wo_salt.m

 %==========================================================================
 %
-% This program demonstrate the use of two independent fluids in two
-% independent domains (one in water, the other in air), for simulation of 
-% membrane for Kerstin.
-%
-% It is important to note that all variables in each of the domains (water 
-% and air) are discretized independently, but temperature is coupled 
-% through the system matrix.
-%
-% A nice thing about this implementation is that it didn't need
-% modifications in any of the supporting functions (for example in "Core")
-%
-% An interesting thing, for the future, that it can serve as a basis for
-% implemetation of two-fluid model (but it has nothing to do with Kerstin's
-% thesis).
-%
 %==========================================================================
 clear
 clc;
@@ -86,10 +71,10 @@ pi= 3.142;
 
 % Membrane properties
 mem.d      = 166E-6;  % m thickness
-mem.lambda = 0.2;    % W/mK thermal conductivity
-mem.eps    = 0.687;    % porosity
-mem.tau    = 2;      % tortuosity
-mem.r      = 0.1E-6; % m pore radius
+mem.lambda = 0.2;     % W/mK thermal conductivity
+mem.eps    = 0.687;   % porosity
+mem.tau    = 2;    % tortuosity
+mem.r      = 0.1E-6;  % m pore radius
 
 % membrane values
 mem.p  = zeros(nx(AIR), 1);
@@ -98,10 +83,6 @@ mem.a  = zeros(nx(AIR), 1);
 mem.pv = zeros(nx(AIR), 1);
 mem.j  = zeros(nx(AIR), 1);
 
-% time-stepping parameters
-dt  =    0.0001;  % time step
-ndt =     20000;  % number of time steps
-
 % create unknowns; names, positions and sizes
 for c=H2O:FIL 
   u(c) = create_unknown('u-velocity',   'x', [nxm(c) ny(c) ], 'd');
@@ -140,21 +121,22 @@ M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
 u(H2O).bnd(W).type(1,1:ny(H2O)) = 'd'; u(H2O).bnd(W).val(1,1:ny(H2O)) = 4.*0.1.*((1:ny(H2O))).*(ny(H2O)+1-(1:ny(H2O)))/(ny(H2O)+1)^2;
 u(H2O).bnd(E).type(1,1:ny(H2O)) = 'o'; u(H2O).bnd(E).val(1,1:ny(H2O)) = 4.*0.1.*((1:ny(H2O))).*(ny(H2O)+1-(1:ny(H2O)))/(ny(H2O)+1)^2;
 
+
 u(AIR).bnd(W).type(1,1:ny(AIR)) = 'd'; u(AIR).bnd(W).val(1,1:ny(AIR)) =   0.0;
 u(AIR).bnd(E).type(1,1:ny(AIR)) = 'd'; u(AIR).bnd(E).val(1,1:ny(AIR)) =  +0.0;
 
 t(H2O).bnd(W).type(1,1:ny(H2O)) = 'd'; t(H2O).bnd(W).val(1,1:ny(H2O)) = 80.0;
-t(H2O).bnd(E).type(1,1:ny(H2O)) = 'n'; t(H2O).bnd(E).val(1,1:ny(H2O)) = 60.0;
-t(H2O).bnd(S).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(S).val(1:nx(H2O),1) = 60.0;
-t(H2O).bnd(N).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(N).val(1:nx(H2O),1) = 60.0;
+t(H2O).bnd(E).type(1,1:ny(H2O)) = 'n';
+t(H2O).bnd(S).type(1:nx(H2O),1) = 'd'; t(H2O).bnd(S).val(1:nx(H2O),1) = 80.0;
+t(H2O).bnd(N).type(1:nx(H2O),1) = 'n';
 t(H2O) = adjust_n_boundaries(t(H2O));
 
 t_int_mem = t(H2O).bnd(S).val;
 
-t(AIR).bnd(W).type(1,1:ny(AIR)) = 'n'; t(AIR).bnd(W).val(1,1:ny(AIR)) = 70.0;
-t(AIR).bnd(E).type(1,1:ny(AIR)) = 'n'; t(AIR).bnd(E).val(1,1:ny(AIR)) = 70.0;
+t(AIR).bnd(W).type(1,1:ny(AIR)) = 'n';
+t(AIR).bnd(E).type(1,1:ny(AIR)) = 'n';
 t(AIR).bnd(S).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(S).val(1:nx(AIR),1) = 20.0;
-t(AIR).bnd(N).type(1:nx(AIR),1) = 'n'; t(AIR).bnd(N).val(1:nx(AIR),1) = 70.0;
+t(AIR).bnd(N).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(N).val(1:nx(AIR),1) = 70.0;
 t(AIR) = adjust_n_boundaries(t(AIR));
 
 t(FIL).bnd(S).type(1:nx(FIL),1) = 'd'; t(FIL).bnd(S).val(1:nx(FIL),1) = 20.0;
@@ -175,12 +157,19 @@ a(AIR).bnd(S).val = p_v.bnd(S).val./(p(AIR).tot(:,1) + 1E5)*M_H2O./M.bnd(S).val;
 % define location of obstacles (e.g. membrane)
 obst = [];
 
+% load('workspace_sim_mem_20000_60');
+
+% time-stepping parameters
+dt  =    0.001;  % time step
+ndt =    22000;  % number of time steps
+n_start = 1;
+
 %-----------
 %
 % time loop 
 %
 %-----------
-for n=1:ndt
+for n=n_start:ndt
 
   disp_time_step(n);
   
@@ -230,8 +219,8 @@ for n=1:ndt
   mem.M  = (M.bnd(N).val      + M.val(:,ny(AIR)))     /2.0;
   
   % Diffusion Coefficients
-  C_K = 2*mem.eps*mem.r/(3*mem.tau*mem.d*dy(AIR)).*(8*M_H2O./(mem.t*R*pi));
-  C_M = mem.eps*mem.p.*prop(AIR).diff(:,ny(AIR))./(mem.tau*R*mem.t.*(mem.p-mem.pv)*dy(AIR));
+  C_K = 2*mem.eps*mem.r/(3*mem.tau*mem.d).*(8*M_H2O./(mem.t*R*pi)).^0.5;
+  C_M = mem.eps*mem.p.*prop(AIR).diff(:,ny(AIR))./(mem.tau*R*mem.t.*(mem.p-mem.pv));
   C_T = 1.0./(1.0./C_K + 1.0./C_M);
   
   % jump condition at membrane
@@ -347,12 +336,7 @@ for n=1:ndt
                                      (prop(FIL).rho.*prop(FIL).cap)./dt, ... % innertial term
                                      prop(FIL).lambda,                   ... % diffusive coefficent
                                      dx(FIL), dy(FIL),                   ... % dimensions
-                                     obst);                              ... % obstacles
- 
-  % combine two matrices
-  A_t = [ [ A_t_air ; zeros(nx(H2O)*ny(H2O), nx(AIR)*ny(AIR)) ], ...
-          [ zeros(nx(AIR)*ny(AIR), nx(H2O)*ny(H2O)) ; A_t_h2o ] ];
-  b_t = [ b_t_air; b_t_h2o ];    
+                                     obst);                              ... % obstacles  
   
   for i=1:nx(H2O)
     %>>>>>>>>>>>>>>>>>>>>>>
@@ -363,22 +347,19 @@ for n=1:ndt
     j = ny(AIR);
     [c_air, w_air, e_air, s_air, n_air] = loc(i,j,nx(AIR),ny(AIR));
     
+    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,ny(AIR))*mem.eps;
+    cy = -prop(AIR).lambda(i,ny(AIR)) * 2*dx(AIR)/dy(AIR) ...
+      + dx(AIR)/(mem.d/lambda_mem + dy(AIR)/(2*prop(AIR).lambda(i,ny(AIR))));
+    
+    b_t_air(c_air) = b_t_air(c_air) + cy * t(AIR).bnd(N).val(i,1) ;
+    A_t_air(c_air, c_air) = A_t_air(c_air, c_air) + cy ...
+      - h_d*mem.j(i).*(1-mem.flowdir(i))/t(AIR).val(i,j);
+    
     % number of the unknown in the water domain (on top of air)
     j = 1;
     [c_h2o, w_h2o, e_h2o, s_h2o, n_h2o] = loc(i,j,nx(H2O),ny(H2O));
-    c_h2o = c_h2o + nx(AIR)*ny(AIR);
-    s_h2o = s_h2o + nx(AIR)*ny(AIR);
-    
-    lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(i,j)*mem.eps;
-    cy = dx(AIR)/(dy(H2O)/(2*prop(H2O).lambda(i,1)) + mem.d/lambda_mem ...
-      + dy(AIR)/(2*prop(AIR).lambda(i,ny(AIR))));
-    
-    A_t(c_air, n_air) = -cy;
-    A_t(c_air, c_air) = A_t(c_air, c_air) + cy ...
-      - h_d*mem.j(i).*(1-mem.flowdir(i))/t(AIR).val(i,j);
 
-    A_t(c_h2o, s_h2o) = -cy;
-    A_t(c_h2o, c_h2o) = A_t(c_h2o, c_h2o) + cy ...
+    A_t_h2o(c_h2o, c_h2o) = A_t_h2o(c_h2o, c_h2o) ...
       + h_d*mem.j(i).*mem.flowdir(i)/t(H2O).val(i,j);
     
     %<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -397,26 +378,26 @@ for n=1:ndt
   c_t_air = advection(prop(AIR).rho.*prop(AIR).cap, t(AIR), u(AIR), v(AIR), dx(AIR), dy(AIR), dt, 'superbee');
   c_t_h2o = advection(prop(H2O).rho.*prop(H2O).cap, t(H2O), u(H2O), v(H2O), dx(H2O), dy(H2O), dt, 'superbee');
   c_t_fil = advection(prop(FIL).rho.*prop(FIL).cap, t(FIL), u(FIL), v(FIL), dx(FIL), dy(FIL), dt, 'superbee');
-  c_t = [ c_t_air; c_t_h2o ];
     
   % innertial term for enthalpy; form separatelly then combine
   i_t_air = reshape( (t(AIR).val.*prop(AIR).rho.*prop(AIR).cap/dt*(dx(AIR)*dy(AIR))), nx(AIR)*ny(AIR),1);
   i_t_h2o = reshape( (t(H2O).val.*prop(H2O).rho.*prop(H2O).cap/dt*(dx(H2O)*dy(H2O))), nx(H2O)*ny(H2O),1);
   i_t_fil = reshape( (t(FIL).val.*prop(FIL).rho.*prop(FIL).cap/dt*(dx(FIL)*dy(FIL))), nx(FIL)*ny(FIL),1);
-  i_t = [ i_t_air; i_t_h2o ];
   
   % the entire source term
-  f_t = b_t - c_t + i_t;
+  f_t_air = b_t_air - c_t_air + i_t_air;
+  f_t_h2o = b_t_h2o - c_t_h2o + i_t_h2o;
   f_t_fil = b_t_fil - c_t_fil + i_t_fil;
   
   % solve temperature
-  t_all      = A_t\(f_t);
+  t(AIR).val = A_t_air\(f_t_air);
+  t(H2O).val = A_t_h2o\(f_t_h2o);
   t(FIL).val = A_t_fil\(f_t_fil);
   
   % split the temperature in air and water domain
-  t(AIR).val = reshape(t_all(                1 : nx(AIR)*ny(AIR)),                 nx(AIR), ny(AIR));
-  t(H2O).val = reshape(t_all(nx(AIR)*ny(AIR)+1 : nx(AIR)*ny(AIR)+nx(H2O)*ny(H2O)), nx(H2O), ny(H2O));
-  t(FIL).val = reshape(t(FIL).val,                                                 nx(FIL), ny(FIL));
+  t(AIR).val = reshape(t(AIR).val, nx(AIR), ny(AIR));
+  t(H2O).val = reshape(t(H2O).val, nx(H2O), ny(H2O));
+  t(FIL).val = reshape(t(FIL).val, nx(FIL), ny(FIL));
   
   t(AIR) = adjust_n_boundaries(t(AIR));
   t(H2O) = adjust_n_boundaries(t(H2O));
@@ -548,7 +529,7 @@ for n=1:ndt
     % temperature
     subplot(3,2,1);
       [C,h] = contourf(x, y, [t(FIL).val'; t(AIR).val'; t(H2O).val'], ...
-                       linspace(19.99,60.01,10)); 
+                       linspace(19.99,80.01,10)); 
       axis equal
       title(sprintf('Temperature at time step %d out of %d\n', n, ndt));
     % velocity magnitude
@@ -584,6 +565,10 @@ for n=1:ndt
     drawnow;
   end
   
+  if mod(n,500) == 0
+    save('workspace_sim_mem.mat')
+  end
+  
 end
 
 % % plot in tecplot format