Hello all, I have a system of (1D) equations I would like to solve:
Note that the first equation is actually only evaluated at the top of the grid (indicated in the code by subscript 'fs'). From a previous suggestion I went with the method of lines, discretising in space first and then solving the set of ODEs using ode15s. However, I am running into a problem with respect to the time step size. 'Unable to meet integration tolerances without reducing the step size below the smallest value allowed'. For reference, this is the ode function file:
function dydt = ode_system(t,y,nx,T_fs,T_w,rho_fr_bc1,rho_fr_bc2,rhoa_fun,drhoadT_fun,dQdT_fun,ddQdTT_fun,k_fr_fun,cp_fr_fun,D_eff,h_m,W_a,W_fs,L_sub,k_ice)
x_fs = y(1); dx = x_fs/nx; T_fr = y(2:nx-1); rho_fr = y((nx-1)+1:end);
k_fr = k_fr_fun(rho_fr); rho_a = rhoa_fun(T_fr);
dTfsdx = gradient([T_fr ; T_fs],dx); dTfsdx = dTfsdx(end);
dxfsdt = 1./(rho_fr_bc2)*(h_m*(W_a - W_fs) - D_eff*rhoa_fun(T_fs)*dQdT_fun(T_fs)*dTfsdx);
dTdx = zeros(size(T_fr)); ddTdxx = dTdx; drhodx = dTdx; dkdx = dTdx; drhoadx = dTdx; dQdx = dTdx;
drhoadT = drhoadT_fun(T_fr); dQdT = dQdT_fun(T_fr); ddQdTT = ddQdTT_fun(T_fr);
dTdx(1) = (T_fr(3) - T_w)/(2*dx); dTdx(nx-2) = (T_fs - T_fr(nx-3))/(2*dx); ddTdxx(1) = (T_fr(3) - 2*T_fr(2) + T_w)/(dx^2); ddTdxx(nx-2) = (T_fs - 2*T_fr(nx-2) + T_fr(nx-3))/(dx^2); drhodx(1) = (rho_fr(3) - rho_fr_bc1)/(2*dx); drhodx(nx-2) = (rho_fr_bc2 - rho_fr(nx-3))/(2*dx); dkdx(1) = (k_fr(3) - k_ice)/(2*dx); dkdx(nx-2) = (k_fr_fun(rho_fr_bc2) - k_fr(nx-3))/(2*dx); drhoadx(1) = drhoadT(1)*dTdx(1); drhoadx(nx-2) = drhoadT(nx-2)*dTdx(nx-2); dQdx(1) = dQdT(1)*dTdx(1); dQdx(nx-2) = dQdT(nx-2)*dTdx(nx-2); ddQdxx(1) = ddQdTT(1)*ddTdxx(1); ddQdxx(nx-2) = ddQdTT(nx-2)*ddTdxx(nx-2);
for i = 3:nx-3 dTdx(i-1) = (T_fr(i+1) - T_fr(i-1))/(2*dx); ddTdxx(i-1) = (T_fr(i+1) - 2*T_fr(i) + T_fr(i-1))/(dx^2); drhodx(i-1) = (rho_fr(i+1) - rho_fr(i-1))/(2*dx); dkdx(i-1) = (k_fr(i+1) - k_fr(i-1))/(2*dx); drhoadx(i-1) = drhoadT(i-1)*( dTdx(i-1) ); dQdx(i-1) = dQdT(i-1)*dTdx(i-1); ddQdxx(i-1) = ddQdTT(i-1)*ddTdxx(i-1); end
dTdt = ( 1./(cp_fr_fun(T_fr).*rho_fr) ).*( dkdx.*dTdx + k_fr.*ddTdxx + L_sub.*drhodx ); drhodt = -D_eff.*( drhoadx.*dQdx + rho_a.*ddQdxx);
dydt = [dxfsdt ; dTdt ; drhodt];
Now I am not sure what exactly the problem is. I strongly suspect it has to do with the anonymous functions used by the ode function, which are these:
syms rhoa_fun(T_fun) rhov_fun(T_fun) Q_fun(T_fun) k_fr_fun(rho_fun)
rhoa_fun = @(T_fun) ((p_0 - (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15)))))./(R_s.*T_fun))+((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun));
rhov_fun = @(T_fun) (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun);
k_fr_fun = @(rho_fun) (rho_fun - rho_ice)./(rho_air - rho_ice)*(k_air - k_ice) + k_ice;
cp_fr_fun = @(rho_fun) (rho_fun - rho_ice)./(rho_air - rho_ice)*(c_p_air - c_p_ice) + c_p_ice;
drhoadT_fun_sym(T_fun) = diff(rhoa_fun,T_fun);
drhoadT_fun = matlabFunction(drhoadT_fun_sym);
Q_fun = @(T_fun) ((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun))./(((p_0 - (611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15)))))./(R_s.*T_fun))+((611.21 .* exp((22.587.*(T_fun-273.15))./(273.86+(T_fun-273.15))))./(R_sv.*T_fun)));
dQdT_fun_sym(T_fun) = diff(Q_fun,T_fun);
dQdT_fun = matlabFunction(dQdT_fun_sym); ddQdTT_fun_sym(T_fun) = diff(Q_fun,T_fun,2);
ddQdTT_fun = matlabFunction(ddQdTT_fun_sym);
Although it might also relate to my usage of the gradient function? I use this as a quick fix as I lack information beyond the grid. After ignoring dxfsdt in the ode system, however, the step size warning still pops up.
If anyone has experience with this problem, I would appreciate the help!
