Solve a complex differential equations system

I need to write a script that can solve the system of differential equation written in this article.

Thank you to everyone who will try to help me:)

[Torsten]

And what is your specific question ?

Answer 1

[Torsten]

The main problem with your equations is that they involve a term that has to be evaluated at a boundary point. Otherwise, you could have tried to use "pdepe". You could also look into "pde1m" under

https://uk.mathworks.com/matlabcentral/fileexchange/97437-pde1dm?tab=reviews

I'm not sure if you have access to the values in all grid points with this code.

If it cannot solve your equations, you will have to discretize the spatial derivatives, include the boundary conditions and use "ode15s" to solve the resulting system of ordinary differential equations. Look up "method of lines" for more details.

[Diego Rampi]

Hi, thanks for your reply, sorry but I hadn't seen it. I have tried to implement the resolution script in various ways over time but I never succeeded (I haven't thought about it for several months). I'll show you some examples, all of which are non-functional.

%%Prova 4



clear
clc

% Parametri fisici e costanti per il Decano
H = 15; % Altezza del serbatoio in metri
diametro = H; % Diametro del serbatoio in metri
V_serbatoio = pi * (diametro / 2)^2 * H; % Volume del serbatoio
V_90 = 0.9 * V_serbatoio; % 90% del volume del serbatoio
P = 1e5; % Pressione in Pa
M = 0.14; % Massa molare del Decano in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)

% Proprietà fisiche del Decano (in aria)
k = 0.026; % Conducibilità termica dell'aria in W/m·K
D = 1.1e-5; % Diffusività in m²/s
cp = 2250; % Calore specifico in J/kg·K
d_h = diametro; % Diametro idraulico del serbatoio
mu = 1.8e-5; % Viscosità dinamica dell'aria (Pa·s)

% Costanti di Antoine per il calcolo della tensione di vapore del Decano
A_ant = 0.21021;
B_ant = 440.616;
C_ant = -156.896;

% Costante di Henry per il Decano in aria (mol/(m^3*Pa))
H_c = 2e-6;

% Numero di nodi di discretizzazione 
N = 150; % numero di nodi
dz = H / (N - 1); % Passo spaziale
z = linspace(0, H, N)'; % Vettore spaziale

% Tempo massimo (30 giorni) in secondi
t_max = 30 * 24 * 3600; % 30 giorni in secondi

% Definizione dei tempi per le fasi di riempimento, riposo, e svuotamento
t_riempimento = 24 * 3600; % 24 ore in secondi
t_riposo = 5 * 24 * 3600; % 5 giorni in secondi
t_svuotamento = 24 * 3600; % 24 ore in secondi

% Velocità di riempimento e svuotamento
v_riempimento = V_90 / (pi * (diametro / 2)^2 * t_riempimento); % Velocità di riempimento in m/s
v_svuotamento = -V_90 / (pi * (diametro / 2)^2 * t_svuotamento); % Velocità di svuotamento in m/s

% Numero di cicli
cycles = 3;

% Inizializzazione delle condizioni iniziali
T_init = 300 * ones(N, 1); % Temperatura iniziale in K
C_init = linspace(0.02, 0.01, N)'; % Gradiente iniziale della concentrazione

% Vettore delle condizioni iniziali combinate
y0 = [T_init; C_init];

% Tempo di transizione per passaggi morbidi tra le fasi
transition_time =4*3600; % 1 ora di transizione

% Funzione che definisce la velocità media del liquido nel tempo con transizioni morbide
v_j_func = @(t) v_riempimento * (1 + tanh((t_riempimento - mod(t, t_riempimento + t_riposo + t_svuotamento)) / transition_time)) / 2 + ...
                v_svuotamento * (1 + tanh((mod(t, t_riempimento + t_riposo + t_svuotamento) - (t_riempimento + t_riposo)) / transition_time)) / 2;
% Visualizzazione dell'andamento della velocità
t_total = cycles * (t_riempimento + t_riposo + t_svuotamento); % Tempo totale in secondi
t_vector = linspace(0, t_total, 1000); % Vettore del tempo per la visualizzazione
v_j_values = arrayfun(v_j_func, t_vector); % Calcolo della velocità v_j(t)

figure;
plot(t_vector / 3600, v_j_values, '-o'); % Tempo in ore
xlabel('Tempo (ore)');
ylabel('Velocità (m/s)');
%title('Andamento della velocità di riempimento e svuotamento');
grid on;

% Definizione delle condizioni al contorno variabili nel tempo
T_top_func = 295; % Temperatura costante al tetto
T_ILG = 300; % Inizialmente assumiamo T_ILG costante, ma può variare

% Risoluzione del sistema con ode15s con tolleranze modificate
options = odeset('RelTol',1e-6,'AbsTol',1e-7, 'MaxStep', 5000);

[t, y] = ode15s(@(t, y) odefun(t, y, z, dz, N, R_gas, P, M, T_top_func, T_ILG, v_j_func, H, ...
                k, D, cp, d_h, mu, A_ant, B_ant, C_ant, H_c), ...
                [0 cycles * (t_riempimento + t_riposo + t_svuotamento)], y0, options);

Warning: Failure at t=7.413934e-05. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.168404e-19) at time t.

% Separazione delle soluzioni
T_sol = y(:, 1:N);
C_sol = y(:, N+1:end);

% Calcolo del rapporto C(top,t)/C_sat
T_top = T_sol(:, end);
P_vap = (10 .^ (A_ant - B_ant ./ (T_top + C_ant))) * 10^5; % Tensione di vapore con equazione di Antoine convertita in Pa
C_sat = P_vap ./ H_c; % Calcolo di Csat con la legge di Henry
C_top = C_sol(:, end);
C_top_over_Csat = C_top ./ C_sat; % Rapporto senza vincoli fisici

% Grafico del rapporto C(top,t)/Csat nel tempo
figure;
plot(t / 3600, C_top_over_Csat, '-o'); % Tempo in ore
xlabel('Tempo (ore)');
ylabel('C(top,t)/C_{sat}');
%title('Andamento di C(top,t)/C_{sat} nel tempo');
grid on;

% Funzione ODE per la risoluzione del sistema
function dydt = odefun(t, y, z, dz, N, R_gas, P, M, T_top_func, T_ILG, v_j_func, H, ...
                        k, D, cp, d_h, mu, A_ant, B_ant, C_ant, H_c)
    % Estrazione delle variabili
    T = y(1:N); % Temperatura
    C = y(N+1:end); % Concentrazione
    
    % Calcolo di rho come funzione di T
    rho = P * M ./ (R_gas * T); % Densità variabile

    % Calcolo di v* (velocità del vapore)
    v_j = v_j_func(t); % Velocità media del liquido
    v_bar = v_j; % Assumendo per semplicità che la velocità media sia costante
    dCdz = gradient(C, dz);
    dTdz = gradient(T, dz);
    v_star = max(0, v_j - D ./ rho .* dCdz); % Calcolo di v* limitato a valori positivi

    % Numeri di Reynolds, Schmidt e Prandtl
    nu = mu ./ rho; % Viscosità cinematica
    Sc = nu ./ D; % Numero di Schmidt
    Pr = nu ./ (k ./ (rho .* cp)); % Numero di Prandtl
    Re = v_star .* d_h .* rho ./ nu; % Numero di Reynolds

    % Calcolo di E, E_t e alpha come funzioni di rho
    E = v_star .* d_h .* (1./(Re.*Sc) + (Re.*Sc)./192); % Coefficiente di dispersione assiale
    E_t = k * (1 + 0.011 * v_bar * d_h * Pr); % Coefficiente di dispersione termica
    alpha = E_t ./ (rho .* cp); % Calcolo di alpha, variabile nel tempo e nello spazio

    % Calcolo della tensione di vapore e C_sat per la condizione al contorno
    P_vap = (10 .^ (A_ant - B_ant ./ (T(1) + C_ant))) * 10^5; % Tensione di vapore alla base
    C_sat = P_vap ./ H_c; % Concentrazione di saturazione

    % Condizioni al contorno per temperatura
    T(1) = T_ILG;  % Temperatura all'interfaccia
    T(end) = T_top_func;  % Temperatura al tetto

    % Condizioni al contorno per concentrazione
    C(1) = C_sat;  % Concentrazione di saturazione alla base
    
    % Condizione al contorno per la derivata di C durante exhale
    if v_j < 0 % Fase di exhale
        % Assegniamo al valore dell'ultimo nodo la derivata calcolata
        dCdz_end = (C(end) * v_star(end) * (H - v_j * t)) / E(end); 
        % Ora correggiamo l'ultimo valore usando la derivata calcolata
        C(end) = C(end-1) + dCdz_end * dz;
    end

    % Termini comuni per le equazioni
    common_factor = 1 / (v_bar * t - H);

    % Equazione della temperatura (Eq. 13)
    d2Tdz2 = del2(T, dz);
    term_T = common_factor^2 .* (T .* d2Tdz2 - dTdz.^2 + dTdz(1) .* dTdz);
    dTdt = (z .* v_bar + (v_star - v_bar)) .* common_factor .* dTdz + term_T .* (alpha ./ T_ILG);

    % Equazione della concentrazione (Eq. 14)
    d2Cdz2 = del2(C, dz); % Calcolo della derivata seconda di C
    term_C_T = common_factor^2 .* ((dTdz - dTdz(1)) .* dCdz + C .* d2Tdz2);
    dCdt = (z .* v_bar + (v_star - v_bar)) .* common_factor .* dCdz + common_factor^2 .* E .* d2Cdz2 - term_C_T .* (alpha ./ T_ILG);

    % Unione dei risultati
    dydt = [dTdt; dCdt];
end

[Diego Rampi]

clear;
clc;

% Parametri fisici e costanti per il Decano
H = 14; % Altezza del serbatoio in metri
diametro = H; % Diametro del serbatoio in metri
V_serbatoio = pi * (diametro / 2)^2 * H; % Volume del serbatoio
V_90 = 0.9 * V_serbatoio; % 90% del volume del serbatoio
P = 1e5; % Pressione in Pa
M = 0.14; % Massa molare del Decano in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)
k = 0.026; % Conducibilità termica dell'aria in W/m·K
D = 1.1e-5; % Diffusività in m²/s
cp = 2250; % Calore specifico in J/kg·K
mu = 1.8e-5; % Viscosità dinamica dell'aria (Pa·s)

% Costanti di Antoine per il calcolo della tensione di vapore del Decano
A_ant = 0.21021;
B_ant = 440.616;
C_ant = -156.896;

% Costante di Henry per il Decano in aria (mol/(m^3*Pa))
H_c = 2e-6;

% Densità dell'aria
rho_aria = 1.225; % kg/m^3

% Parametri di discretizzazione
N = 150; % Numero di nodi spaziali
dz = H / (N - 1); % Passo spaziale
dt = 3600; % Passo temporale di 1 ora
t_max = 30 * 24 * 3600; % 30 giorni in secondi
z = linspace(0, H, N)'; % Discretizzazione spaziale
time_steps = t_max / dt; % Numero di passi temporali
t_vec = linspace(0, t_max, time_steps); % Vettore temporale

% Definizione dei tempi per le fasi di riempimento, riposo e svuotamento
t_riempimento = 24 * 3600; % 24 ore per riempimento
t_riposo = 5 * 24 * 3600; % 5 giorni di riposo
t_svuotamento = 24 * 3600; % 24 ore per svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento; % Durata di un ciclo completo

% Condizioni iniziali
C_top = zeros(time_steps, 1); % Concentrazione al tetto
T = 300 * ones(N, 1); % Temperatura iniziale uniforme (K)
C = linspace(0.02, 0.01, N)'; % Gradiente iniziale della concentrazione
C_sat = zeros(time_steps, 1); % Concentrazione di saturazione

% Parametri del serbatoio
T_top_func = 295; % Temperatura costante al tetto

% Parametri per il calcolo di T_ILG
h_L = 10; % Coefficiente di scambio termico del liquido (W/m^2·K)
DeltaH_evap = 33000; % Calore di evaporazione del Decano (J/mol)
m_L = V_90 * 700; % Massa del liquido (kg) con densità stimata a 700 kg/m^3

% Funzione per la velocità del serbatoio (riempimento, riposo e svuotamento)
transition_time = 3600; % Tempo di transizione in secondi
v_avg_func = @(t) v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);

% Metodo di Eulero esplicito per la concentrazione e temperatura nel tempo
for t_idx = 1:time_steps
    t = t_vec(t_idx);
    
    % Metodo di punto fisso per risolvere le non linearità (iterazioni sul sistema)
    tolerance = 1e-5;
    max_iter = 50;
    
    % Iterazioni per trovare la soluzione stabile
    for iter = 1:max_iter
        % Calcolo di rho in funzione di T
        rho = P * M ./ (R_gas * T);
        
        % Calcolo della velocità media (v_avg) e della velocità del vapore v*
        dCdz = gradient(C, dz);
        v_avg = v_avg_func(t);
        v_star = v_avg - (D ./ rho) .* dCdz(end);
        
        % Approssimazione dell'indice per il gradiente di temperatura (per evitare errore)
        index = max(1, min(N, round(v_avg * t / dz)));
        
        % Calcolo della velocità v in funzione di v* e del gradiente di temperatura
        dTdz = gradient(T, dz);
        v = v_star + (k / (T(1) * rho)) .* (dTdz - dTdz(index));
        
        % Numero di Reynolds, Schmidt e Prandtl
        nu = mu ./ rho;
        Sc = nu ./ D;
        Pr = nu ./ (k ./ (rho .* cp));
        Re = v .* diametro .* rho ./ nu;
        
        % Dispersione assiale e termica
        E = v .* diametro .* (1./(Re .* Sc) + (Re .* Sc)./192);
        E_t = k * (1 + 0.011 .* v .* diametro .* Pr);
        
        % Calcolo dei gradienti per le equazioni differenziali
        d2Tdz2 = del2(T, dz);
        d2Cdz2 = del2(C, dz);
        
        % Calcolo di C_sat come massa di Decano per kg di aria
        P_vap = (10^(A_ant - B_ant / (T(1) + C_ant))) * 10^5; % Pa
        C_sat(t_idx) = P_vap * H_c; % mol(dec)/m^3(aria)
        
        % Calcolo di rho_ILG all'interfaccia
        rho_ILG = C_sat(t_idx) * M; % kg(dec)/m^3(aria)
        
        % Calcolo di alpha come E_t / (rho_ILG * cp)
        alpha = E_t ./ (rho_ILG * cp);
        
        % **Calcolo della temperatura T_ILG dal bilancio energetico**
        T_ILG = (k * dTdz(index) + D * dCdz(index) * DeltaH_evap + h_L * T_top_func) ./...
                (h_L + E_t / (rho_ILG * cp));

        % Equazione della Temperatura (Eulero esplicito)
        T_new = T;
        term_T = (v_avg * t - H)^(-2) * (T .* d2Tdz2 - dTdz.^2 + dTdz(1) * dTdz);
        dTdt = (z .* v_avg + (v_star - v_avg)) .* (v_avg * t - H)^(-1) .* dTdz + term_T .* (alpha ./ T_ILG);
        
        % **Aggiorna la temperatura con Eulero esplicito**
        T_new(2:end-1) = T(2:end-1) + dt * dTdt(2:end-1);
        
        % Equazione della Concentrazione (Eulero esplicito)
        C_new = C;
        term_C_T = (v_avg * t - H)^(-2) .* ((dTdz - dTdz(1)) .* dCdz + C .* d2Tdz2);
        dCdt = (z .* v_avg + (v_star - v_avg)) .* (v_avg * t - H)^(-1) .* dCdz + (v_avg * t - H)^(-2) .* E .* d2Cdz2 - term_C_T .* (alpha ./ T_ILG);
        
        % **Aggiorna la concentrazione con Eulero esplicito**
        C_new(2:end-1) = C(2:end-1) + dt * dCdt(2:end-1);
        
        % Condizioni al contorno
        T_new(1) = T_ILG;  % Temperatura alla base
        T_new(end) = T_top_func;  % Temperatura al tetto
        C_new(1) = C_sat(t_idx); % Concentrazione alla base
        
        % Condizione al contorno sulla derivata prima di C
        if v_avg > 0
            % Exhale: la derivata di C è nulla in y=1
            dCdz_end = 0;
        else
            % Inhale: la derivata di C segue la condizione fornita
            dCdz_end = (C(end) * v(end) * (H - v_avg * t)) / E(end);
        end
        C_new(end) = C(end-1) + dCdz_end * dz;
        
        % Controllo di convergenza
        if all(abs(C_new - C) < tolerance) && all(abs(T_new - T) < tolerance)
            break;
        end
        
        % Aggiorna le variabili
        T = T_new;
        C = C_new;
    end
    
    % Salva C(top,t)/Csat
    C_top(t_idx) = C(end) / C_sat(t_idx);
end

Unable to perform assignment because the left and right sides have a different number of elements.

% Visualizzazione della concentrazione C(top,t)/Csat nel tempo
figure;
plot(t_vec / 3600, C_top, '-o');
xlabel('Tempo (ore)');
ylabel('C(top,t) / C_{sat}');
title('Andamento di C(top,t) / C_{sat} nel tempo');
grid on;

% Visualizzazione dell'andamento della velocità
figure;
plot(t_vec / 3600, v_avg_func(t_vec), '-o');
xlabel('Tempo (ore)');
ylabel('Velocità (m/s)');
title('Andamento della velocità del serbatoio');
grid on;

%% Funzione per calcolare la velocità durante i cicli di riempimento-riposo-svuotamento
function v = v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time)
    v_riempimento = V_90 / (pi * (diametro / 2)^2 * t_riempimento);
    v_svuotamento = -V_90 / (pi * (diametro / 2)^2 * t_svuotamento);
    
    cycle_time = mod(t, cycle_duration);
    
    if cycle_time <= t_riempimento
        v = v_riempimento * (1 + tanh((t_riempimento - cycle_time) / transition_time)) / 2;
    elseif cycle_time <= (t_riempimento + t_riposo)
        v = 0;
    else
        v = v_svuotamento * (1 + tanh((cycle_time - (t_riempimento + t_riposo)) / transition_time)) / 2;
    end
end

[Diego Rampi]

This is the system of PDEs:

[Diego Rampi]

Same issues using pde function:

clear;
clc;

% Parametri del serbatoio
H = 14; % Altezza del serbatoio in metri
diametro = H; % Diametro del serbatoio in metri
V_serbatoio = pi * (diametro / 2)^2 * H; % Volume del serbatoio
V_90 = 0.9 * V_serbatoio; % 90% del volume del serbatoio
P = 1e5; % Pressione in Pa

% Parametri fisici e costanti del Decano
M = 0.14; % Massa molare del Decano in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)
k = 0.026; % Conducibilità termica dell'aria in W/m·K
D = 1.1e-5; % Diffusività in m²/s
cp = 2250; % Calore specifico in J/kg·K
mu = 1.8e-5; % Viscosità dinamica dell'aria (Pa·s)
cpL = 2000; % Calore specifico del Decano liquido in J/kg·K
h_L = 10; % Coefficiente di scambio termico del liquido W/m²·K
rho_liq = 740; % Densità del Decano liquido a 298 K in kg/m³
DeltaH_evap = 40000; % Calore latente di evaporazione in J/mol

% Costanti di Antoine per il calcolo della tensione di vapore del Decano
A_ant = 0.21021;
B_ant = 440.616;
C_ant = -156.896;

% Parametri di discretizzazione
N = 150; % Numero di nodi spaziali
dz = H / (N - 1); % Passo spaziale
t_max = 30 * 24 * 3600; % 30 giorni in secondi
z = linspace(0, H, N); % Discretizzazione spaziale
t_vec = linspace(0, t_max, 200); % Aumento numero di punti temporali

% Inizializzazione delle condizioni iniziali
T_init = 300 * ones(N, 1); % Temperatura iniziale uniforme (K)
C_init = linspace(0.02, 0.01, N)'; % Gradiente iniziale della concentrazione

% Vettore iniziale combinato per T e C
y0 = [T_init; C_init];

% Velocità media del liquido
transition_time = 3600; % 1 ora di transizione
t_riempimento = 24 * 3600; % 24 ore di riempimento
t_riposo = 5 * 24 * 3600; % 5 giorni di riposo
t_svuotamento = 24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento;

% Funzione che definisce la velocità media del liquido nel tempo
v_avg_func = @(t) v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);

% Opzioni del solver pdepe con tolleranze modificate
%options = odeset('RelTol',1e-4,'AbsTol',1e-5, 'MaxStep', 1e3);

% Risolvi il sistema di PDEs usando pdepe
m = 0; % Geometria 1D
sol = pdepe(m, @pde_fun, @ic_fun, @(xl, ul, xr, ur, t) bc_fun(xl, ul, xr, ur, t, V_90, diametro, H, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time, P, M, R_gas, mu, D, k, cp), z, t_vec);%, options);

Warning: Failure at t=8.640000e+04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.328306e-10) at time t.

Warning: Time integration has failed. Solution is available at requested time points up to t=7.815075e+04.

% Visualizzazione della soluzione
T_sol = sol(:, :, 1); % Estrae la soluzione di T
C_sol = sol(:, :, 2); % Estrae la soluzione di C

% Calcolo di C_top / C_sat
C_top_over_Csat = C_sol(:, end) ./ C_sol(:, 1);

% Grafico C(top,t)/C_sat nel tempo
figure;
plot(t_vec / 3600, C_top_over_Csat, '-o'); % Grafico nel tempo (ore)

Error using plot
Specify the coordinates as vectors or matrices of the same size, or as a vector and a matrix that share the same length in at least one dimension.

xlabel('Tempo (ore)');
ylabel('C(top,t) / C_{sat}');
title('Andamento di C(top,t) / C_{sat} nel tempo');
grid on;

%% Funzione PDE per il sistema accoppiato (T e C)
function [c, f, s] = pde_fun(x, t, u, dudx)
    % Decomponi la soluzione in temperatura e concentrazione
    T = u(1); % Temperatura
    C = u(2); % Concentrazione
    
    % Gradienti spaziali (dT/dx e dC/dx)
    dTdz = dudx(1);
    dCdz = dudx(2);
    
    % Parametri fisici
    P = 1e5; % Pressione
    M = 0.14; % Massa molare del Decano
    R_gas = 8.314; % Costante dei gas
    k = 0.026; % Cond. termica
    D = 1.1e-5; % Diffusività
    
    rho = P * M / (R_gas * T); % Densità variabile
    alpha = k / (rho * 2250); % Diffusività termica
    
    % Definisci i coefficienti della PDE
    c = [1; 1]; % Termini temporali per T e C
    f = [alpha * dTdz; D * dCdz]; % Termini di diffusione
    s = [0; 0]; % Termini sorgente
    
    % Accoppiamenti tra T e C
    s(1) = -D / rho * dCdz * dTdz;
end

%% Calcolo di T_ILG basato sull'equazione fornita
function T_ILG = calc_T_ILG(T_top_func, dTdz, dCdz, t, h_L, cpL, rho_liq, DeltaH_evap, V_liq)
    % Calcolo della massa del liquido
    m_liq = rho_liq * V_liq; % kg
    
    % Calcola T_ILG con il bilancio termico
    T_ILG = ((k * dTdz + D * dCdz * DeltaH_evap) + h_L * T_top_func) / (h_L + m_liq * cpL);
end

%% Condizioni iniziali
function u0 = ic_fun(x)
    T0 = 300; % Temperatura iniziale (K)
    C0 = 0.02; % Concentrazione iniziale
    u0 = [T0; C0];
end

%% Condizioni al contorno
function [pl, ql, pr, qr] = bc_fun(xl, ul, xr, ur, t, V_90, diametro, H, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time, P, M, R_gas, mu, D, k, cp)
    % Condizioni al contorno per T e C
    % Condizioni alla base (y=0)
    pl = [ul(1) - 300; ul(2) - 0.02];  % T(0)=300, C(0)=C_sat
    ql = [0; 0];
    
    % Calcola la velocità media (v_avg)
    v_avg = v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);
    
    % Condizioni al contorno in y=1 (z=H)
    if v_avg > 0  % Exhale: derivata nulla della concentrazione
        pr = [ur(1) - 295; 0];  % T(H)=295, dC/dy(H) = 0
        qr = [0; 1];  % La derivata di C è zero per exhale
    else  % Inhale: derivata non nulla per C
        rho = P * M / (R_gas * ur(1)); % Densità variabile
        nu = mu / rho;
        Sc = nu / D;
        Re = v_avg * diametro * rho / nu;
        E = v_avg * diametro * (1 / (Re * Sc) + (Re * Sc) / 192);
        
        % Derivata non nulla per C durante inhale
        pr = [ur(1) - 295; (H - v_avg * t) / E];  % T(H)=295, dC/dy(H) = (H - v_avg * t) / E
        qr = [0; 1];  % dC/dy non nullo
    end
end

%% Funzione per calcolare la velocità durante i cicli di riempimento-riposo-svuotamento
function v = v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time)
    % Calcolo della velocità di riempimento e svuotamento con transizioni
    v_riempimento = V_90 / (pi * (diametro / 2)^2 * t_riempimento);
    v_svuotamento = -V_90 / (pi * (diametro / 2)^2 * t_svuotamento);
    
    % Determinazione della fase del ciclo (riempimento, riposo, svuotamento)
    cycle_time = mod(t, cycle_duration);
    
    if cycle_time <= t_riempimento
        % Fase di riempimento
        v = v_riempimento * (1 + tanh((t_riempimento - cycle_time) / transition_time)) / 2;
    elseif cycle_time <= (t_riempimento + t_riposo)
        % Fase di riposo (velocità zero)
        v = 0;
    else
        % Fase di svuotamento
        v = v_svuotamento * (1 + tanh((cycle_time - (t_riempimento + t_riposo)) / transition_time)) / 2;
    end
end

[Torsten]

As you can see, "pdepe" has problems when coefficients in your PDE system change discontinuously. So integrate up to 86400 (24*3600), change parameters and restart with the last solution. Do the same if there are similar time instances following.

But I'm not sure it makes sense starting with "pdepe" since your equations don't permit a solution with this solver. I'd suggest using "pde1dm" right from the beginning. Or the code in which you use ode15s and the method of lines.

[Diego Rampi]

the "pde1dm" solver seems to be to rigid for my problem because allows the solution of a PDE structured in the way reported at page 2 of the manual.

I don't now if is possible to write the script in a way that can solve my system, reported below:

\left\{\begin{array}{l}

\frac{\partial \mathrm{T}}{\partial \mathrm{t}}-\frac{\mathrm{y} \bar{v}+\left(v^{*}-\bar{v}\right)}{\bar{v} \mathrm{t}-\mathrm{H}} \frac{\partial \mathrm{T}}{\partial \mathrm{y}}-\frac{\alpha}{\mathrm{T}_{\text {ILG }}} \frac{1}{(\bar{v} \mathrm{t}-\mathrm{H})^{2}}\left(\mathrm{~T} \frac{\partial^{2} \mathrm{~T}}{\partial \mathrm{y}^{2}}-\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right)^{2}+\left.\frac{\partial \mathrm{T}}{\partial \mathrm{y}} \frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right|_{\mathrm{y}=0}\right)=0 \tag{13}\\

\frac{\partial \mathrm{C}}{\partial \mathrm{t}}-\frac{\mathrm{y} \bar{v}+\left(v^{*}-\bar{v}\right)}{\bar{v} \mathrm{t}-\mathrm{H}} \frac{\partial \mathrm{C}}{\partial \mathrm{y}}-\frac{\mathrm{E}}{(\bar{v} \mathrm{t}-\mathrm{H})^{2}} \frac{\partial^{2} \mathrm{C}}{\partial \mathrm{y}^{2}}+\frac{\alpha}{\mathrm{T}_{\text {ILG }}} \frac{1}{(\bar{v} \mathrm{t}-\mathrm{H})^{2}}\left(\frac{\partial \mathrm{C}}{\partial \mathrm{y}}\left(\frac{\partial \mathrm{T}}{\partial \mathrm{y}}-\left.\frac{\partial \mathrm{T}}{\partial \mathrm{y}}\right|_{\mathrm{y}=0}\right)+\mathrm{C} \frac{\partial^{2} \mathrm{~T}}{\partial \mathrm{y}^{2}}\right)=0

\end{array}\right.

Anyone as an idea?

I really appreciate your help, thank you!

[Diego Rampi]

This is my code for now:

clear;
clc;

% Parametri del serbatoio
H = 14; % Altezza del serbatoio in metri
diametro = H; % Diametro del serbatoio in metri
V_serbatoio = pi * (diametro / 2)^2 * H; % Volume del serbatoio
V_90 = 0.9 * V_serbatoio; % 90% del volume del serbatoio
P = 1e5; % Pressione in Pa

% Parametri fisici e costanti del Decano
M = 0.14; % Massa molare del Decano in kg/mol
R_gas = 8.314; % Costante dei gas perfetti in J/(mol*K)
k = 0.026; % Conducibilità termica dell'aria in W/m·K
D = 1.1e-5; % Diffusività in m²/s
cp = 2250; % Calore specifico in J/kg·K
mu = 1.8e-5; % Viscosità dinamica dell'aria (Pa·s)
cpL = 2000; % Calore specifico del Decano liquido in J/kg·K
h_L = 10; % Coefficiente di scambio termico del liquido W/m²·K
rho_liq = 740; % Densità del Decano liquido a 298 K in kg/m³
DeltaH_evap = 40000; % Calore latente di evaporazione in J/mol

% Costanti di Antoine per il calcolo della tensione di vapore del Decano
A_ant = 0.21021;
B_ant = 440.616;
C_ant = -156.896;
H_c = 2e-6; % Costante di Henry per il Decano in aria (mol/(m^3*Pa))
TL = 298; % Temperatura del liquido costante a 298 K

% Parametri di discretizzazione
N = 150; % Numero di nodi spaziali
dz = H / (N - 1); % Passo spaziale
t_max = 30 * 24 * 3600; % 30 giorni in secondi
z = linspace(0, H, N); % Discretizzazione spaziale
t_vec = linspace(0, t_max, 200); % Punti temporali

% Condizioni iniziali
T_init = 300 * ones(N, 1); % Temperatura iniziale uniforme (K)
C_init = linspace(0.02, 0.01, N)'; % Gradiente iniziale della concentrazione
y0 = [T_init; C_init];

% Velocità media del liquido
transition_time = 3600; % 1 ora di transizione
t_riempimento = 24 * 3600; % 24 ore di riempimento
t_riposo = 5 * 24 * 3600; % 5 giorni di riposo
t_svuotamento = 24 * 3600; % 24 ore di svuotamento
cycle_duration = t_riempimento + t_riposo + t_svuotamento;

% Funzione per la velocità media nel tempo
v_avg_func = @(t) v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time);

% Risolvi il sistema di PDE usando pde1dm (inseriamo m=0 per geometria 1D)
m = 0; % Geometria 1D cartesiana
sol = pde1dm(m, @(x,t,u,dudx) pde_fun(x,t,u,dudx, diametro, v_avg_func, D, mu, cp), @ic_fun, @(xl,ul,xr,ur,t) bc_fun(xl,ul,xr,ur,t,v_avg_func, V_90, diametro, H, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time, P, M, R_gas, mu, D, k, cp, DeltaH_evap, A_ant, B_ant, C_ant, H_c, TL, h_L, cpL, rho_liq), z, t_vec); % Chiamata pde1dm

Warning: Failure at t=8.640000e+04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.069545e-10) at time t.

% Visualizzazione della soluzione
T_sol = sol(:, :, 1); % Estrae la soluzione di T
C_sol = sol(:, :, 2); % Estrae la soluzione di C

% Calcolo di C_top / C_sat
C_top_over_Csat = C_sol(:, end) ./ C_sol(:, 1);

% Grafico di C(top,t)/C_sat nel tempo
figure;
plot(t_vec / 3600, C_top_over_Csat, '-o'); % Tempo in ore

Error using plot
Specify the coordinates as vectors or matrices of the same size, or as a vector and a matrix that share the same length in at least one dimension.

xlabel('Tempo (ore)');
ylabel('C(top,t) / C_{sat}');
title('Andamento di C(top,t) / C_{sat} nel tempo');
grid on;

%% Funzione PDE per il sistema accoppiato (T e C)
function [c, f, s] = pde_fun(x, t, u, dudx, diametro, v_avg_func, D, mu, cp)
    T = u(1); % Temperatura
    C = u(2); % Concentrazione
    dTdz = dudx(1);
    dCdz = dudx(2);
    
    P = 1e5; % Pressione
    M = 0.14; % Massa molare
    R_gas = 8.314; % Costante dei gas
    k = 0.026; % Cond. termica
    
    % Calcolo della densità
    rho = P * M / (R_gas * T); % Densità variabile
    
    % Calcolo della velocità media e v_star
    v_avg = v_avg_func(t);
    v_star = v_avg - (D / rho) * dCdz; % Calcolo di v_star
    
    % Calcolo dei numeri adimensionali
    nu = mu / rho; % Viscosità cinematica
    Sc = nu / D; % Numero di Schmidt
    Pr = nu / (k / (cp * rho)); % Numero di Prandtl
    Re = v_avg * diametro * rho / mu; % Numero di Reynolds
    
    % Coefficiente di dispersione assiale E
    E = v_avg * diametro * (1 / (Re * Sc) + (Re * Sc) / 192);
    
    % Coefficiente di dispersione termica E_t
    E_t = v_avg * diametro * (1 / (Re * Pr) + (Re * Pr) / 192);
    
    % Diffusività termica
    alpha = E_t / (rho * cp); % Diffusività termica
    
    % Velocità del vapore v
    v = v_star + alpha / T * (dTdz - dTdz(1)); % Velocità del vapore
    
    % Derivate temporali
    c = [1; 1];
    
    % Derivate spaziali (prime e seconde)
    f = [alpha * dTdz; D * dCdz]; % Usa E e D per calcolare i contributi spaziali
    
    % Termini sorgente accoppiati
    s = [0; -D / rho * dCdz * dTdz]; % Accoppiamento tra T e C
end

%% Funzione per le condizioni iniziali
function u0 = ic_fun(x)
    T0 = 300;
    C0 = 0.02;
    u0 = [T0; C0];
end

%% Funzione per le condizioni al contorno
function [pl, ql, pr, qr] = bc_fun(xl, ul, xr, ur, t, v_avg_func, V_90, diametro, H, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time, P, M, R_gas, mu, D, k, cp, DeltaH_evap, A_ant, B_ant, C_ant, H_c, TL, h_L, cpL, rho_liq)
    % Calcola la velocità media
    v_avg = v_avg_func(t);
    
    % Calcolo di dTdz e dCdz in y=0
    dTdz_0 = (ul(1) - ur(1)) / (H / (length(ul) - 1)); % Approssimazione della derivata prima per T
    dCdz_0 = (ul(2) - ur(2)) / (H / (length(ul) - 1)); % Approssimazione della derivata prima per C

    % Calcolo dinamico di T_ILG e C_sat tramite Antoine e legge di Henry
    T_ILG = calc_T_ILG(TL, dTdz_0, dCdz_0, h_L, k, cpL, rho_liq, DeltaH_evap, v_avg, diametro, D);
    P_vap = (10^(A_ant - B_ant / (T_ILG + C_ant))) * 1e5; % Pa
    C_sat = P_vap * H_c; % mol/m³ (dalla legge di Henry)

    % Condizioni al contorno alla base (y=0): T = T_ILG, C = C_sat
    pl = [ul(1) - T_ILG; ul(2) - C_sat];  % T(0)=T_ILG, C(0)=C_sat
    ql = [0; 0]; % Condizione di Dirichlet

    % Condizioni al contorno al tetto (y=1)
    if v_avg > 0  % Exhale: derivata nulla della concentrazione
        pr = [ur(1) - 295; 0];  % T(H)=295, dC/dy(H) = 0
        qr = [0; 1];  % Condizione Neumann per la concentrazione
    else  % Inhale: derivata non nulla per C
        rho = P * M / (R_gas * ur(1)); % Densità variabile
        nu = mu / rho;
        Sc = nu / D;
        Re = v_avg * diametro * rho / nu;
        E = v_avg * diametro * (1 / (Re * Sc) + (Re * Sc) / 192);

        % Condizione per inhale: derivata non nulla della concentrazione
        pr = [ur(1) - 295; (H - v_avg * t) / E];
        qr = [0; 1];
    end
end

%% Funzione per calcolare T_ILG basato sul bilancio termico
function T_ILG = calc_T_ILG(TL, dTdz, dCdz, h_L, k, cpL, rho_liq, DeltaH_evap, v_avg, diametro, D)
    % Calcolo della massa del liquido in base alla velocità media
    V_liq = v_avg * (pi * (diametro / 2)^2); % Volume del liquido
    m_liq = rho_liq * V_liq; % Massa del liquido

    % Calcolo di T_ILG tramite bilancio termico
    T_ILG = (h_L * (TL - 300) + m_liq * cpL * (TL - 300) + k * dTdz + D * dCdz * DeltaH_evap) / (h_L + m_liq * cpL);
end

%% Funzione per calcolare la velocità durante i cicli
function v_avg = v_cicli(t, V_90, diametro, t_riempimento, t_riposo, t_svuotamento, cycle_duration, transition_time)
    v_riempimento = V_90 / (pi * (diametro / 2)^2 * t_riempimento);
    v_svuotamento = -V_90 / (pi * (diametro / 2)^2 * t_svuotamento);
    
    cycle_time = mod(t, cycle_duration);
    
    if cycle_time <= t_riempimento
        v_avg = v_riempimento * (1 + tanh((t_riempimento - cycle_time) / transition_time)) / 2;
    elseif cycle_time <= (t_riempimento + t_riposo)
        v_avg = 0;
    else
        v_avg = v_svuotamento * (1 + tanh((cycle_time - (t_riempimento + t_riposo)) / transition_time)) / 2;
    end
end