Skip to content
GitLab
Commit b9c50b34 authored by Kerstin CRAMER's avatar Kerstin CRAMER
Browse files

include continue script & jump condition at membrane

parent 8613bea7
Expand all Hide whitespace changes
Inline Side-by-side
MATLAB/Model/sim_mem.asv
View file @ b9c50b34
......@@ -100,7 +100,7 @@ mem.j = zeros(nx(AIR), 1);
% time-stepping parameters
dt = 0.0001; % time step
ndt = 2000; % number of time steps
ndt = 20000; % number of time steps
% create unknowns; names, positions and sizes
for c=H2O:FIL
......@@ -123,13 +123,16 @@ p(FIL).tot = zeros(nx(FIL), ny(FIL));
% specify initial conditions
t(H2O).val(1:nx(H2O), 1:ny(H2O)) = 80;
t(FIL).val(1:nx(FIL), 1:ny(FIL)) = 20;
t(AIR).val(1:nx(AIR), 1:ny(AIR)) = 50;
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) = 45;
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);
M.val = M_AIR;
a(AIR).val(1:nx(AIR), 1:ny(AIR)) = p_v_sat(min(min((t(AIR).val))*1E-5*M_H2O./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;
a(AIR).val(1:nx(AIR), 1:ny(AIR)) = p_v_sat(t(AIR).val)*1E-5*M_H2O./M.val;
M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
% domain specific boundary conditions
% note that north of air domain and south of water domain are neumann for
......@@ -141,11 +144,13 @@ 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)) = 70.0;
t(H2O).bnd(S).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(S).val(1:nx(H2O),1) = 70.0;
t(H2O).bnd(N).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(N).val(1:nx(H2O),1) = 70.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) = 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(S).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(S).val(1:nx(AIR),1) = 20.0;
......@@ -159,13 +164,13 @@ t(FIL) = adjust_n_boundaries(t(FIL));
p_v.bnd(S).type(1:nx(AIR),1) = 'd'; p_v.bnd(S).val(1:nx(AIR),1) = p_v_sat(t(FIL).bnd(N).val);
p_v.bnd(N).type(1:nx(AIR),1) = 'd'; p_v.bnd(N).val(1:nx(AIR),1) = p_v_sat(t(H2O).bnd(S).val);
M.bnd(S).type(1:nx(AIR),1) = 'd'; M.bnd(S).val(1:nx(AIR),1) = M_AIR;
M.bnd(N).type(1:nx(AIR),1) = 'd'; M.bnd(N).val(1:nx(AIR),1) = M_AIR;
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(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
a(AIR).bnd(S).val = p_v.bnd(S).val./(p(AIR).tot(:,1)+1E5)*M_H2O./M.bnd(S).val;
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
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 = [];
......@@ -193,18 +198,16 @@ for n=1:ndt
% AIR domain values
M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
p_v.val = a(AIR).val .* M.val/M_H2O .* (p(AIR).tot + 1E5);
p_v.val = a(AIR).val .* M.val/M_H2O .* (p(AIR).tot + 1E+5);
% AIR N bnd values
p_v.bnd(N).val = p_v_sat(t(H2O).bnd(S).val); % rather T_interface
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)+1E5)*M_H2O./M.bnd(N).val; % p_v/p_tot * M_H2O/M_mix
% liquid film & AIR S bnd values
M.bnd(S).val = M_AIR*(1-a(AIR).bnd(S).val) + M_H2O*a(AIR).bnd(S).val;
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E5);
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E+5);
M.bnd(S).val = M_AIR*(1-a(AIR).bnd(S).val) + M_H2O*a(AIR).bnd(S).val;
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E+5);
t_int_film = t_sat(p_v.bnd(S).val);
disp(sprintf('mean t_int_film = %5.3f', mean(t_int_film)));
t(AIR).bnd(S).val = t_int_film;
t(FIL).bnd(N).val = t_int_film;
......@@ -214,28 +217,56 @@ for n=1:ndt
m_out = ( 2*prop(AIR).lambda(:,1)./dy(AIR).*(t_int_film - t(AIR).val(:,1))...
+ 2*prop(FIL).lambda(:,ny(FIL))./dy(FIL).*(t_int_film - t(FIL).val(:,ny(FIL)))) ...
* dx(AIR) / h_d; % kg/(m*s)
m_out = m_out;
%------------------------
% membrane
%------------------------
% membrane values
mem.p = (p(H2O).tot(:,1)+p(AIR).tot(:,ny(AIR)))/2.0 + 1E5;
mem.t = (t(H2O).val(:,1)+t(AIR).val(:,ny(AIR)))/2.0 + 273.15;
mem.a = (a(H2O).val(:,1)+a(AIR).val(:,ny(AIR)))/2.0;
mem.pv = (p_v.bnd(N).val +p_v.val(:,ny(AIR))) /2.0;
mem.M = (M.bnd(N).val +M.val(:,ny(AIR))) /2.0;
mem.p = (p(H2O).tot(:,1) + p(AIR).tot(:,ny(AIR)))/2.0 + 1E5;
mem.t = (t(H2O).bnd(S).val + t(AIR).val(:,ny(AIR)))/2.0 + 273.15;
mem.a = (a(AIR).bnd(N).val + a(AIR).val(:,ny(AIR)))/2.0;
mem.pv = (p_v.bnd(N).val + p_v.val(:,ny(AIR))) /2.0;
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_T = 1.0./(1.0./C_K + 1.0./C_M);
% mem.j = C_T .* dx(AIR) .* (p_v.bnd(N).val - p_v.val(:,ny(AIR))); % kg/(m*s)
mem.j(:) = 1E-7; % kg/(m*s)
% jump condition at membrane
lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(:,ny(AIR))*mem.eps;
lhs_lin_mem = ( 2*prop(H2O).lambda(:,1)./dy(H2O) ...
+ 1./(dy(AIR)./(2*prop(AIR).lambda(:,ny(AIR))) + mem.d./lambda_mem)) ...
* mem.eps * dx(AIR) / h_d;
lhs_fun_mem = C_T*dx(AIR);
rhs_mem = C_T*dx(AIR).*p_v.val(:,ny(AIR)) + mem.eps * dx(AIR)/h_d ...
* ( 2*prop(H2O).lambda(:,1)./dy(H2O).*t(H2O).val(:,1) ...
+ 1./(dy(AIR)./(2*prop(AIR).lambda(:,ny(AIR)))+mem.d./lambda_mem) ...
.* 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);
t_int_mem(i) = fzero(jump_cond_mem,t(H2O).val(i,1));
end
disp(sprintf('mean t_int_mem = %5.3f', mean(t_int_mem)));
% 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);
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
mem.j = C_T .* dx(AIR) .* (p_v.bnd(N).val - p_v.val(:,ny(AIR))); % kg/(m*s)
% mem.j(:) = 1E-7; % kg/(m*s)
mem.flowdir = zeros(nx(AIR),1);
mem.flowdir = 1.0*(mem.j >= 0.0);
%------------------------
% concentration
%------------------------
......@@ -355,10 +386,10 @@ for n=1:ndt
%<<<<<<<<<<<<<<<<<<<<<<<<<<
% condensation at liquid film
j = 1;
[c_air, w_air, e_air, s_air, n_air] = loc(i,j,nx(AIR),ny(AIR));
j = ny(FIL);
[c_fil, w_fil, e_fil, s_fil, n_fil] = loc(i,j,nx(FIL),ny(FIL));
A_t(c_air, c_air) = A_t(c_air, c_air) - h_d*m_out(i)/t(AIR).val(i,j);
A_t_fil(c_fil, c_fil) = A_t_fil(c_fil, c_fil) - h_d*m_out(i)/t(FIL).val(i,j);
end
......@@ -494,8 +525,6 @@ for n=1:ndt
disp(sprintf('volume imbalance after correction in phase %1d = %12.5e', c, max(max(abs(b_p)))));
end
mean(t_int_film)
%--------------------
% plot them together
......@@ -519,7 +548,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,80.01,10));
linspace(19.99,60.01,10));
axis equal
title(sprintf('Temperature at time step %d out of %d\n', n, ndt));
% velocity magnitude
......
MATLAB/Model/sim_mem.m
View file @ b9c50b34
......@@ -98,9 +98,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 = 2000; % number of time steps
% create unknowns; names, positions and sizes
for c=H2O:FIL
......@@ -121,15 +118,18 @@ p(FIL).tot = zeros(nx(FIL), ny(FIL));
M.name = 'molar mass'; M.pos = 'c'; M.val = zeros(nx(AIR), ny(AIR));
% specify initial conditions
t(H2O).val(1:nx(H2O), 1:ny(H2O)) = 60;
t(H2O).val(1:nx(H2O), 1:ny(H2O)) = 80;
t(FIL).val(1:nx(FIL), 1:ny(FIL)) = 20;
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)) = 50;
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);
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;
a(AIR).val(1:nx(AIR), 1:ny(AIR)) = p_v_sat(t(AIR).val)*1E-5*M_H2O./M.val;
M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
% domain specific boundary conditions
% note that north of air domain and south of water domain are neumann for
......@@ -140,12 +140,14 @@ u(H2O).bnd(E).type(1,1:ny(H2O)) = 'o'; u(H2O).bnd(E).val(1,1:ny(H2O)) = 4.*0.1.*
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)) = 60.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) = 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(S).type(1:nx(AIR),1) = 'd'; t(AIR).bnd(S).val(1:nx(AIR),1) = 20.0;
......@@ -159,23 +161,30 @@ t(FIL) = adjust_n_boundaries(t(FIL));
p_v.bnd(S).type(1:nx(AIR),1) = 'd'; p_v.bnd(S).val(1:nx(AIR),1) = p_v_sat(t(FIL).bnd(N).val);
p_v.bnd(N).type(1:nx(AIR),1) = 'd'; p_v.bnd(N).val(1:nx(AIR),1) = p_v_sat(t(H2O).bnd(S).val);
M.bnd(S).type(1:nx(AIR),1) = 'd'; M.bnd(S).val(1:nx(AIR),1) = M_AIR;
M.bnd(N).type(1:nx(AIR),1) = 'd'; M.bnd(N).val(1:nx(AIR),1) = M_AIR;
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(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
a(AIR).bnd(S).val = p_v.bnd(S).val./(p(AIR).tot(:,1)+1E5)*M_H2O./M.bnd(S).val;
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
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.0001; % time step
ndt = 200000; % number of time steps
% n_start = 20001;
%-----------
%
% time loop
%
%-----------
for n=1:ndt
for n=n_start:ndt
disp_time_step(n);
......@@ -193,18 +202,16 @@ for n=1:ndt
% AIR domain values
M.val = M_AIR*(1-a(AIR).val) + M_H2O*a(AIR).val;
p_v.val = a(AIR).val .* M.val/M_H2O .* (p(AIR).tot + 1E5);
p_v.val = a(AIR).val .* M.val/M_H2O .* (p(AIR).tot + 1E+5);
% AIR N bnd values
p_v.bnd(N).val = p_v_sat(t(H2O).bnd(S).val); % rather T_interface
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)+1E5)*M_H2O./M.bnd(N).val; % p_v/p_tot * M_H2O/M_mix
% liquid film & AIR S bnd values
M.bnd(S).val = M_AIR*(1-a(AIR).bnd(S).val) + M_H2O*a(AIR).bnd(S).val;
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E5);
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E+5);
M.bnd(S).val = M_AIR*(1-a(AIR).bnd(S).val) + M_H2O*a(AIR).bnd(S).val;
p_v.bnd(S).val = a(AIR).bnd(S).val .* M.bnd(S).val/M_H2O .* (p(AIR).tot(:,1) + 1E+5);
t_int_film = t_sat(p_v.bnd(S).val);
disp(sprintf('mean t_int_film = %5.3f', mean(t_int_film)));
t(AIR).bnd(S).val = t_int_film;
t(FIL).bnd(N).val = t_int_film;
......@@ -214,28 +221,56 @@ for n=1:ndt
m_out = ( 2*prop(AIR).lambda(:,1)./dy(AIR).*(t_int_film - t(AIR).val(:,1))...
+ 2*prop(FIL).lambda(:,ny(FIL))./dy(FIL).*(t_int_film - t(FIL).val(:,ny(FIL)))) ...
* dx(AIR) / h_d; % kg/(m*s)
m_out = m_out;
%------------------------
% membrane
%------------------------
% membrane values
mem.p = (p(H2O).tot(:,1)+p(AIR).tot(:,ny(AIR)))/2.0 + 1E5;
mem.t = (t(H2O).val(:,1)+t(AIR).val(:,ny(AIR)))/2.0 + 273.15;
mem.a = (a(H2O).val(:,1)+a(AIR).val(:,ny(AIR)))/2.0;
mem.pv = (p_v.bnd(N).val +p_v.val(:,ny(AIR))) /2.0;
mem.M = (M.bnd(N).val +M.val(:,ny(AIR))) /2.0;
mem.p = (p(H2O).tot(:,1) + p(AIR).tot(:,ny(AIR)))/2.0 + 1E5;
mem.t = (t(H2O).bnd(S).val + t(AIR).val(:,ny(AIR)))/2.0 + 273.15;
mem.a = (a(AIR).bnd(N).val + a(AIR).val(:,ny(AIR)))/2.0;
mem.pv = (p_v.bnd(N).val + p_v.val(:,ny(AIR))) /2.0;
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_T = 1.0./(1.0./C_K + 1.0./C_M);
% mem.j = C_T .* dx(AIR) .* (p_v.bnd(N).val - p_v.val(:,ny(AIR))); % kg/(m*s)
mem.j(:) = 1E-8; % kg/(m*s)
% jump condition at membrane
lambda_mem = mem.lambda*(1-mem.eps) + prop(AIR).lambda(:,ny(AIR))*mem.eps;
lhs_lin_mem = ( 2*prop(H2O).lambda(:,1)./dy(H2O) ...
+ 1./(dy(AIR)./(2*prop(AIR).lambda(:,ny(AIR))) + mem.d./lambda_mem)) ...
* mem.eps * dx(AIR) / h_d;
lhs_fun_mem = C_T*dx(AIR);
rhs_mem = C_T*dx(AIR).*p_v.val(:,ny(AIR)) + mem.eps * dx(AIR)/h_d ...
* ( 2*prop(H2O).lambda(:,1)./dy(H2O).*t(H2O).val(:,1) ...
+ 1./(dy(AIR)./(2*prop(AIR).lambda(:,ny(AIR)))+mem.d./lambda_mem) ...
.* 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);
t_int_mem(i) = fzero(jump_cond_mem,t(H2O).val(i,1));
end
disp(sprintf('mean t_int_mem = %5.3f', mean(t_int_mem)));
% 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);
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
mem.j = C_T .* dx(AIR) .* (p_v.bnd(N).val - p_v.val(:,ny(AIR))); % kg/(m*s)
% mem.j(:) = 1E-7; % kg/(m*s)
mem.flowdir = zeros(nx(AIR),1);
mem.flowdir = 1.0*(mem.j >= 0.0);
%------------------------
% concentration
%------------------------
......@@ -355,10 +390,10 @@ for n=1:ndt
%<<<<<<<<<<<<<<<<<<<<<<<<<<
% condensation at liquid film
j = 1;
[c_air, w_air, e_air, s_air, n_air] = loc(i,j,nx(AIR),ny(AIR));
j = ny(FIL);
[c_fil, w_fil, e_fil, s_fil, n_fil] = loc(i,j,nx(FIL),ny(FIL));
A_t(c_air, c_air) = A_t(c_air, c_air) - h_d*m_out(i)/t(AIR).val(i,j);
A_t_fil(c_fil, c_fil) = A_t_fil(c_fil, c_fil) - h_d*m_out(i)/t(FIL).val(i,j);
end
......@@ -494,8 +529,6 @@ for n=1:ndt
disp(sprintf('volume imbalance after correction in phase %1d = %12.5e', c, max(max(abs(b_p)))));
end
mean(t_int_film)
%--------------------
% plot them together
......
MATLAB/Model/sim_mem_continue.m 0 → 100644
View file @ b9c50b34
This diff is collapsed.
MATLAB/Model/sim_mem_t_film_connected.asv → MATLAB/Model/sim_mem_wo_jump_mem.m
View file @ b9c50b34
......@@ -16,162 +16,168 @@
% thesis).
%
%==========================================================================
clear
clc;
%------------------
% phase indicators
%------------------
H2O = 1;
AIR = 2;
FIL = 3;
%------------------------------
% set path to common functions
%------------------------------
path(path, '../Util');
path(path, '../Core');
path(path, '../IO');
path(path, '../Model');
path(path, '../Global');
%----------------
% define problem
%----------------
% set domain dimensions and grid resolution
lx(H2O)= 0.1; nx(H2O)=150;
ly(H2O)= 0.02; ny(H2O)= 30;
lx(AIR)= 0.1; nx(AIR)=150;
ly(AIR)= 0.01; ny(AIR)= 15;
lx(FIL)= 0.1; nx(FIL)=150;
ly(FIL)= 0.001; ny(FIL)= 2;
for c=H2O:FIL
dx(c)=lx(c)/nx(c); nxm(c)=nx(c)-1;
dy(c)=ly(c)/ny(c); nym(c)=ny(c)-1;
end
u_max = 0.12;
% set physical properties
[prop(H2O).rho, ...
prop(H2O).mu, ...
prop(H2O).cap, ...
prop(H2O).lambda] = water_properties(1:nx(H2O), 1:ny(H2O));
[prop(AIR).rho, ...
prop(AIR).mu, ...
prop(AIR).cap, ...
prop(AIR).lambda] = air_properties(1:nx(AIR), 1:ny(AIR));
[prop(FIL).rho, ...
prop(FIL).mu, ...
prop(FIL).cap, ...
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;
h_d = 2257E3; % J/kg
M_H2O = 18E-3; % kg/mol
M_AIR = 28E-3; % kg/mol
R = 8.314;
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
% membrane values
mem.p = zeros(nx(AIR), 1);
mem.t = zeros(nx(AIR), 1);
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 = 200; % 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');
v(c) = create_unknown('v-velocity', 'y', [nx(c) nym(c)], 'd');
p(c) = create_unknown('pressure', 'c', [nx(c) ny(c) ], 'n');
t(c) = create_unknown('temperature', 'c', [nx(c) ny(c) ], 'n');
a(c) = create_unknown('concentration','c', [nx(c) ny(c) ], 'n');
p(c) = adjust_n_boundaries(p(c));
end
p(AIR).tot = zeros(nx(AIR), ny(AIR));
p(H2O).tot = zeros(nx(H2O), ny(H2O));
p(FIL).tot = zeros(nx(FIL), ny(FIL));
% just for air:
p_v.name = 'vapor pressure'; p_v.pos = 'c'; p_v.val = zeros(nx(AIR), ny(AIR));
M.name = 'molar mass'; M.pos = 'c'; M.val = zeros(nx(AIR), ny(AIR));
% specify initial conditions
t(H2O).val(1:nx(H2O), 1:ny(H2O)) = 80;
t(FIL).val(1:nx(FIL), 1:ny(FIL)) = 20;
t(AIR).val(1:nx(AIR), 1:ny(AIR)) = 50;
t(AIR).val(1:nx(AIR), 1:ny(AIR)/3) = 20;
t(AIR).val(1:nx(AIR), ny(AIR)/3:2*ny(AIR)/3) = 45;
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);
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;
% domain specific boundary conditions
% note that north of air domain and south of water domain are neumann for
% 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;
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)) = 70.0;
t(H2O).bnd(S).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(S).val(1:nx(H2O),1) = 70.0;
t(H2O).bnd(N).type(1:nx(H2O),1) = 'n'; t(H2O).bnd(N).val(1:nx(H2O),1) = 70.0;
t(H2O) = adjust_n_boundaries(t(H2O));
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(S).type(1:nx(AIR),1) = 'n'; 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) = 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;
t(FIL) = adjust_n_boundaries(t(FIL));
p_v.bnd(S).type(1:nx(AIR),1) = 'd'; p_v.bnd(S).val(1:nx(AIR),1) = p_v_sat(t(FIL).bnd(N).val);
p_v.bnd(N).type(1:nx(AIR),1) = 'd'; p_v.bnd(N).val(1:nx(AIR),1) = p_v_sat(t(H2O).bnd(S).val);
M.bnd(S).type(1:nx(AIR),1) = 'd'; M.bnd(S).val(1:nx(AIR),1) = M_AIR;
M.bnd(N).type(1:nx(AIR),1) = 'd'; M.bnd(N).val(1:nx(AIR),1) = M_AIR;
a(AIR).bnd(S).type(1:nx(AIR),1) = 'd';
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
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 = [];
% clear
% clc;
%
% %------------------
% % phase indicators
% %------------------
% H2O = 1;
% AIR = 2;
% FIL = 3;
%
% %------------------------------
% % set path to common functions
% %------------------------------
% path(path, '../Util');
% path(path, '../Core');
% path(path, '../IO');
% path(path, '../Model');
% path(path, '../Global');
%
% %----------------
% % define problem
% %----------------
%
% % set domain dimensions and grid resolution
% lx(H2O)= 0.1; nx(H2O)=150;
% ly(H2O)= 0.02; ny(H2O)= 30;
%
% lx(AIR)= 0.1; nx(AIR)=150;
% ly(AIR)= 0.01; ny(AIR)= 15;
%
% lx(FIL)= 0.1; nx(FIL)=150;
% ly(FIL)= 0.001; ny(FIL)= 2;
%
% for c=H2O:FIL
% dx(c)=lx(c)/nx(c); nxm(c)=nx(c)-1;
% dy(c)=ly(c)/ny(c); nym(c)=ny(c)-1;
% end
%
% u_max = 0.12;
%
% % set physical properties
% [prop(H2O).rho, ...
% prop(H2O).mu, ...
% prop(H2O).cap, ...
% prop(H2O).lambda] = water_properties(1:nx(H2O), 1:ny(H2O));
% [prop(AIR).rho, ...
% prop(AIR).mu, ...
% prop(AIR).cap, ...
% prop(AIR).lambda] = air_properties(1:nx(AIR), 1:ny(AIR));
% [prop(FIL).rho, ...
% prop(FIL).mu, ...
% prop(FIL).cap, ...
% 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(FIL).lambda(:,ny(FIL)) = 1.0;
% % prop(AIR).lambda(:,1) = 0.5;
%
% h_d = 2257E3; % J/kg
%
% M_H2O = 18E-3; % kg/mol
% M_AIR = 28E-3; % kg/mol
%
% R = 8.314;
% pi= 3.142;
%