Esercitazioni Matlab - esame di Metodi numerici per l'ingegneria navale

Esercitazione 1

clc; clear; close all; % Pulisce la command window, cancella le variabili e chiude tutte le figure
tic % Inizia il cronometro per misurare il tempo di esecuzione

Importazione dei dati dal file Excel

filename = ['Aletta_Sez_Rettangolare.xlsx']; % Nome del file Excel (assicurarsi che si trovi nella stessa cartella)
data = readmatrix(filename, 'Sheet', 2); % Legge i dati numerici dal secondo foglio del file Excel

Estrazione delle colonne dai dati

T = data(:,1); % Colonna 1: Temperatura in K
k = data(:,2); % Colonna 2: Conducibilità termica in W/m·K
cp = data(:,3); % Colonna 3: Calore specifico in J/kg·K
Prandtl_Number = data(:,4);% Colonna 4: Numero di Prandtl (adimensionale)
rho = data(:,5); % Colonna 5: Densità in kg/m³
mu = data(:,6); % Colonna 6: Viscosità dinamica in Pa·s

Grafico del numero di Prandtl

figure; % Crea una nuova finestra grafica
plot(T,Prandtl_Number,'ro--','LineWidth',2) % Plotta il numero di Prandtl in funzione della temperatura (punti rossi e linea tratteggiata)
xlabel('Temperatura (°C)'); ylabel('Prandtl Number'); % Etichette degli assi
title('Prandtl Number'); grid on; % Titolo del grafico e griglia
xticks(0:200:1000); % Imposta i tick dell'asse x da 0 a 1000 con passo 200

Creazione di una finestra con 4 subplot per visualizzare i dati originali

figure; % Nuova finestra grafica

Subplot 1: Conducibilità termica

subplot(2,2,1);
plot(T, k, 'b-', 'LineWidth', 2); % Linea blu continua
xlabel('Temperatura (°C)'); ylabel('Conducibilità termica (W/mK)');
title('Conducibilità termica'); grid on;
xticks(0:200:1000);

Subplot 2: Calore specifico

subplot(2,2,2);
plot(T, cp, 'r-', 'LineWidth', 2); % Linea rossa continua
xlabel('Temperatura (°C)'); ylabel('Calore specifico (J/kgK)');
title('Calore specifico'); grid on;
xticks(0:200:1000);

Subplot 3: Densità

subplot(2,2,3);
plot(T, rho, 'g-', 'LineWidth', 2); % Linea verde continua
xlabel('Temperatura (°C)'); ylabel('Densità (kg/m^3)');
title('Densità'); grid on;
xticks(0:200:1000);

Subplot 4: Viscosità

subplot(2,2,4);
plot(T, mu, 'm-', 'LineWidth', 2); % Linea magenta continua
xlabel('Temperatura (°C)'); ylabel('Viscosità (Pa s)');
title('Viscosità'); grid on;
xticks(0:200:1000);

Approssimazioni con metodi di regressione (minimi quadrati)

figure; % Nuova finestra grafica per i fit
Qui creeremo un nuovo grafico con dati + curve di fit.

Conducibilità (fit lineare)

subplot(2,2,1)
[p_k, k_fit] = fitPolinomiale(T, k, 1, 'Temperatura (°C)', 'Conducibilità termica(W/mK)');
p_k vettore dei coefficienti del polinomio [p1, p2].→- k_fit vettore dei valori stimati della conducibilità termica calcolati con la retta del fit→-

Calore specifico (fit quadratico)

subplot(2,2,2)
[p_cp, cp_fit] = fitPolinomiale(T, cp, 2, 'Temperatura (°C)', 'Calore specifico(J/kgK)');

Densità (modello potenza)

subplot(2,2,3)
[p_rho, rho_fit] = fitPotenza(T, rho, 'Temperatura (°C)', 'Densità (kg/m^3)');

Viscosità (modello potenza)

subplot(2,2,4)
[p_mu, mu_fit] = fitPotenza(T, mu, 'Temperatura (°C)', 'Viscosità (Pa·s)');
toc

Funzione 1 — Fit polinomiale

function [p, y_fit] = fitPolinomiale(x, y, grado, xlab, ylab)
p = polyfit(x, y, grado); % Coefficienti polinomio
y_fit = polyval(p, x); % Valori del fit sui dati

% Grafico dati + fit
plot(x, y, 'o', 'DisplayName', 'Dati');
hold on;
plot(x, y_fit, '-', 'LineWidth', 2, 'DisplayName', sprintf('Fit grado %d', grado));
xlabel(xlab); ylabel(ylab); title(sprintf('Fit polinomiale (grado %d)', grado));
legend show; grid on; xticks(0:200:1000);
end

Funzione 2 — Fit modello y = b*x^a

function [p, y_fit] = fitPotenza(x, y, xlab, ylab)
lx = log(x);
ly = log(y); % Linearizzazione log-log
p = polyfit(lx, ly, 1); % Fit lineare log-log
y_fit = exp(polyval(p, lx)); % Ritorno alla scala originale
stiamo facendo Yfit = a⋅ln(x) + ln(b)

% Grafico dati + fit
plot(x, y, 'o', 'DisplayName', 'Dati'); hold on;
plot(x, y_fit, '-', 'LineWidth', 2, 'DisplayName', 'Fit potenza');
xlabel(xlab); ylabel(ylab); title('Fit modello di potenza');
legend show; grid on; xticks(0:200:1000);
end

• polyfit(lx, ly, 1) calcola i coefficienti della retta: slope intercept= , = ln ()
• In p, p(1) = a, p(2) = ln(b)

Esercitazione 2

clc; % Pulisce la Command Window
clear; % Cancella tutte le variabili
close all;

Parametri noti

Treq = 1e5; % Spinta richiesta (N)
Preq = 2e6; % Potenza richiesta (W)
rho = 1025; % Densità dell'acqua di mare (kg/m^3)
KT = 0.2; % Coefficiente di spinta (adimensionale)
KQ = 0.05; % Coefficiente di coppia (adimensionale)

Definizione del sistema di equazioni non lineari

% La variabile x è un vettore: x(1) = n (velocità di rotazione in rps), x(2) = D(diametro elica in m)
fun = @(x) [KT * rho * x(1)^2 * x(2)^4 - Treq; % Equazione della spinta
2 * pi * KQ * rho * x(1)^3 * x(2)^5 - Preq % Equazione della potenza];

Stima iniziale della soluzione

x0 = [5; 2]; % Valori iniziali: n = 5 rps, D = 2 metri

Opzioni per il solutore numerico fsolve

options = optimoptions('fsolve', 'Display', 'iter'); % Visualizza l'output ad ogni iterazione

Risoluzione del sistema

[x_sol, fval, exitflag] = fsolve(fun, x0, options); % Risolve il sistema non lineare

Visualizzazione dei risultati

if exitflag > 0
fprintf('Soluzione trovata:\n');
fprintf('n (velocità rotazione) = %.4f rps\n', x_sol(1));
fprintf('D (diametro elica) = %.4f m\n', x_sol(2));
else warning('Nessuna soluzione trovata. Controllare le condizioni iniziali o le equazioni.');
end

Esercitazione 3

clc
clear
close all

Parametri del sistema

Zita = 0.1; % Coefficiente di smorzamento
Omega_0 = 0.8; % Frequenza naturale [rad/s]
F = 0.2; % Ampiezza della forza eccitante
Omega = 1; % Frequenza dell’eccitazione [rad/s]
y1 = 1; % Condizione iniziale: angolo di rollio [rad]
y2 = 0; % Condizione iniziale: velocità angolare [rad/s]

Risoluzione con ODE45

[t, y] = ode45(@(t,y) Roll_Function(t,y,Zita,Omega_0,F,Omega), [0 50], [y1; y2]);

Plot

figure
plot(t, y(:,1), 'o-b', 'LineWidth', 1.2, 'MarkerSize', 6, 'DisplayName', 'y_1 = Angolo')
hold on
plot(t, y(:,2), 'o-r', 'LineWidth', 1.2, 'MarkerSize', 6, 'DisplayName', 'y_2 = Velocità')
grid on
xlabel('Tempo t [s]', 'FontSize', 12)
ylabel('Angolo e Velocità di Rollio', 'FontSize', 12)
legend('Location','best', 'FontSize', 11)
title('Risposta Dinamica a Rollio Tramite ODE45', 'FontSize', 14, 'FontWeight', 'bold')
axis([0 50 -1.5 1.5])

Funzione del sistema

function dydt = Roll_Function(t, y, Zita, Omega_0, F, Omega)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -2 * Zita * Omega_0 * y(2) - (Omega_0^2) * y(1) + F * cos(Omega * t);
end

Esercitazione 4 - Gas