Solve ODE with function dependent parameters

Hi,

I consider my self a relative novice with MATLAB, but am trying to use it to solve a coupled heat and mass transfer problem using the method of lines, to model drying processes - I expect that as the solution progresses - the "crust" will have different propoerties to the "core" but this isn't happening.

IWhat I hoped would happed is that the physical property functions (which do vary significantly over the course of the solution), would be continuously re-evaluated during the course of the solution (eg varying as the solution progresses), but some basic diagnostics suggest that this is not the case - can you help me see where I'm going wrong?

Many Thanks

Simon

Extract from the ODE function that I'm solving using ode15s - the code runs d=fine, but the solution is suspect.

for i=n:-1:1                       
% Symmetry Boundary Conditions       
if(i==1)
    Tt(i) = 2*(T(i+1) - T(i))  * dx2;           % HT Boundary condition at x=xl, central symmetry
    CWat(i)= 2*(CWa(i+1) - CWa(i)) * dx2;             % MT boundary condition free water at x=xl, central symmetry

% Surface Boundary Conditions                         
elseif(i==n)        % add sorption isotherm here Xgamma as a function of moisture in solid
    Aw = GAB_solve(CWa(n),T(n));
    PSat = AntoineW(T(n));                       % Calculate water vapour pressure
    %Mass fluxes at the interface
    Jwat = (-kmWat * (18.01 / R))  * ( ( (Aw * PSat) / (T(n) - 273.15) ) - (RH / (Tinf -  273.15)) ) ;      
    % Boundary condition at x=xu, convective heat transfer at the surface

    Tt(i) = ((-(h /(lambda(CWa(i-1),T(i-1)))) * (T(i)- Tinf))  - (lambdaW * Jwat) ) * dx2;    % heat transfer boundary condition              
    CWat(i) = Jwat;   % Boundary condition at x=xu, convective heat transfer t=at the surface           % MT boundary condition at x=xu, convective mass transfer at the surface
 
% Body terms
else
    kpot(i) = lambda(CWa(i-1),T(i-1));  % Thermal conductivity as function of composition and Temperature
    rho(i) = density(CWa(i-1),T(i-1));     % Density as function of composition of composition and temperature
    CP(i) = Cp(CWa(i-1),T(i-1));        % Specific Heat Capacity as function of composition and temperature
  
    alphaT = kpot ./ (rho .* CP);      % Thermal diffusivity as function of composition and temperature
    Tt(i) = alphaT(i) * (T(i+1) - 2 * T(i) + T(i-1)) * dx2;             % Diffusive heat transfer

    D = Diffusivity(CWa(i-1),T(i));            % Moisture diffusivity as function of composition and temperature
    CWat(i) = (D * (CWa(i+1) - 2.0 * CWa(i) + CWa(i-1)) * dx2 );      % Diffusive Mass transfer and starch gelatinisation
end

[darova]

Can you attach something more?

[Simon Lawton]

I'm note sure what you need, I'm making some guesses.

The Ode 15 call from the main routine:

[t,u] = ode15s(@(t,u) pdeHTMT_main_coupled_Wa_rect_explicit(t, u, dx, paramHT), tout, u0, options);  

The functions GAB_Solve (GAB sorption isotherm - gives water activity from water composition), AntoineW (saturation pressure from Antoines equation), Lambda (thermal conductivity), density and CP( specific Heat capacity), calculate temperature and composition dependent functions that should change value at each point in space and time. An example of the property function is here:

function[rhoout] = density(X,T)
%% Constants
load('potato_density.mat');
% Potato comps Carbs protein fiber fat ash

% Condition inputs to model units and ranges
The TC = T - 273.15;
DB = (1-X);        % Moisture on a dry weight basis

% Calculate individual component properties
rho.carbs = thermprop(PropCarbs(1), PropCarbs(2), PropCarbs(3), TC);
rho.fat = thermprop(PropFat(1), PropFat(2), PropFat(3), TC);
rho.water = thermprop(PropWater(1), PropWater(2), PropWater(3), TC);
rho.fiber = thermprop(PropFiber(1), PropFiber(2), PropFiber(3), TC);
rho.protein = thermprop(PropProtein(1), PropProtein(2), PropProtein(3), TC);
rho.ash = thermprop(PropAsh(1), PropAsh(2), PropAsh(3), TC);


% Calculate mass fractions for mixing calc
rhofrac.Water = rho.water * X; 
rhofrac.Carbs  = rho.carbs * DB * CompsPotato(1);
rhofrac.Protein = rho.protein * DB * CompsPotato(2);
rhofrac.Fiber = rho.fiber * DB * CompsPotato(3);
rhofrac.Fat = rho.fat * DB * CompsPotato(4);
rhofrac.Ash = rho.ash * DB * CompsPotato(5);

rhoout = rhofrac.Water + rhofrac.Carbs + rhofrac.Protein + rhofrac.Fiber ...
    + rhofrac.Fat + rhofrac.Ash;
end

All of the invividual property functions have been tested independently and return correct values vs literature.

Thanks

Simon

[darova]

Can you show original formulas?

[Simon Lawton]

[darova]

Where is boundary condition (initial)

[Simon Lawton]

INitial conditions are:

T = 293 (Temperature, kelvin scale)

CWa = 0.8 (Moisture concentration)

Answer 1

[darova]

Try this simple example

a = 1/pi^2;
x = linspace(0,1,50)';
T0 = 1-(1/2-x).^2;          % initial condition t=0 (start time)
dx = x(2)-x(1);
% left boundary T=0
% right boundary T=0
f = @(t,T) [0; a*diff(T(1:50),2)/dx^2; 0];
[t,T] = ode45(f,[0 3], T0);
[xx,tt] = meshgrid(x,t);
surf(tt,xx,T,'edgecolor','none')
axis vis3d