Why do I get NaN error for rho(i,n+1) and v(i,n+1) ?

%% The current code assumes that the nozzles are only present at the first control volume (inlet) and the last control volume (outlet).

%% Section 1:  PARAMETERS 

L = 400; % Length of the hose (m)
N = 100; % Number of control volumes

% Operation and Design Properties
Qfad =4.31; % Air flow rate in m^3/s ; Q= 259 m^3/min = 4.31
D_sup = 104e-3; % Diameter of supply hose & riser in meters
D_nozzle= 2e-3; % Diameter of nozzle 
s=0.3; % Nozzle spacing 
C_d= 0.55; % Nozzle discharge coefficient 
P_FAD = 101325; % Pressure 
T_FAD = 293.15; % Temperature in Kelvin
h = 30; % Water depth in meters
Rho_w = 1025; % Density of water in kg/m^3
g = 9.81; % Gravity acceleration in m/s^2
P_atm = 1.013e5; % Atmospheric pressure in Pa
T_inlet = 291.15; % Air temperature in Celsius
R = 287; % Specific gas constant in J/(kg·K)
mu = 18e-6; % Dynamic viscosity of air in Pa·s
T_FAD_ref = 20; % Reference temperature at FAD condition in Celsius
P_FAD_ref = 1.013e5; % Reference pressure at FAD condition in Pa
P_Fed= 1*10^6; % Bar to Pascals 
theta = pi; %Angle of the hose with respect to the horizontal (radians)
A_hose = pi * (D_sup/2)^2;% Area of the hose (m^2)
A_nozzle= pi * (D_nozzle/2)^2;% Area of the nozzle (m^2)
CFL_max = 0.5; % Define maximum allowable CFL number
A_total= A_hose+A_nozzle;

%% SECTION 2:SPATIAL DISCRETIZATION OF THE CONTROL VOLUMES
%spatial discretization to account for the presence of nozzles in the system.
num_nozzles = floor(L / s); % Calculate the total number of nozzles
%The floor function rounds down the result to the nearest integer
total_control_volumes = N + num_nozzles; % Total number of control volumes including nozzles
% Updating Control Volume Boundaries:
x = linspace(0, L, total_control_volumes + 1); % Update the array of control volume boundaries
%0 is the starting point of the control volume boundaries.
%L is the total length of the hose.
%total_control_volumes + 1 specifies the number of equally spaced points (boundaries) to generate between 0 and L, considering the total number of control volumes (including nozzles).
%The result of this line is an array x that contains the positions of the control volume boundaries, including those of the hose and the nozzles.

%Recalculating Control Volume Lengths:
delta_x = x(2:end) - x(1:end-1); 
%x(2:end) extracts the positions of the right boundaries of each control volume (including nozzles).
%x(1:end-1) extracts the positions of the left boundaries of each control volume (including nozzles).
%Subtracting these two arrays element-wise gives you the lengths of each control volume (including nozzles).
%The result of this line is an array delta_x that contains the lengths of each control volume,considering both the hose and the nozzles.

%% SECTION 3: INITIALIZING ARRAYS 
% Number of time steps
Nt = 10000;
% Time step size
dt = 0.01; % Example time step size (you can adjust it based on stability requirements)  
% Initialize arrays to store the variables at each grid point and time step
%By initializing these arrays with zeros, 
% we create placeholders for the future values of the variables during the simulation
rho = zeros(total_control_volumes, Nt);
v = zeros(total_control_volumes, Nt); 
T = zeros(total_control_volumes, Nt);
P = zeros(total_control_volumes, Nt);
m_in=zeros(total_control_volumes, Nt);
m_out=zeros(total_control_volumes, Nt);
v_in =zeros(total_control_volumes, Nt);
v_out=zeros(total_control_volumes, Nt);
m_n= zeros(total_control_volumes, Nt); 
% Time array
time = linspace(0, Nt*dt, Nt+1); % linspace(start, end, N);Nt*dt=total time duration of the simulation.
% This time array allows us to track the simulation's progress over each time step 

%% SECTION 4: INITIAL CONDITION
%These initial conditions are used as starting values for the simulation.
% initial values remain constant throughout the simulation,no need to recompute them in each time step.
rho_initial = 1.225 ;  % Initial density (kg/m^3) (rho= P/R*T)=(1.013e5/287* 293.15)=1.225 
P_initial = 1.013e5;   % Initial pressure (Pa)
T_inital= 293.15;     % Initial Temperature(K)        
v_in_initial=(Qfad / A_total); %velocity of sysytem( Hose + Nozzle)
v_out_initial= (Qfad / (C_d * pi * (D_nozzle/2)^2));
m_initial=rho_initial * A_total*v_in_initial; 
m_out_initial= rho_initial * (A_nozzle / A_total) * v_out_initial;
m_n_initial = rho_initial* (A_nozzle / A_total) * (v_out_initial-v_in_initial);

% Set initial conditions for each variable at each grid point
% array(row_index, column_index)
rho(:, 1) = rho_initial * ones(total_control_volumes, 1);
% So, rho(:, 1) means all rows (all control volumes) of the first column (first time step) in the rho matrix.
P(:, 1) = P_initial * ones(total_control_volumes, 1);
T(:, 1) =T_inital * ones(total_control_volumes, 1);
v(:,1)=v_in_initial * ones(total_control_volumes, 1);
v_in(:,1)=v_in_initial* ones(total_control_volumes, 1);
v_out(:,1)=v_out_initial* ones(total_control_volumes, 1);
m_in(:,1)=m_initial* ones(total_control_volumes, 1);
m_out(:,1)=v_out_initial* ones(total_control_volumes, 1);
m_n(:,1)=m_n_initial* ones(total_control_volumes, 1);
% Print the value of v_in(1, 1) before it's updated in the loop

%% SECTION 5: TIME INTEGRATION LOOP
%The main reasons for solving the time and space dependent parts together are
% related to the interdependence of the variables and the overall accuracy and stability of the solution
for n = 1:Nt %Loop for Time Steps [iterates from the first time step to last time step]
    % This outer loop iterates over time steps. 
    % Within each time step, you update the variables (density, velocity, and pressure) based on the discretized continuity and momentum equations. 
    % We have included the necessary calculations and updates to account for variations across the control volumes along the length of the hose.  
    % Inlet Boundary Conditions
    P(1, n+1) = 1.33e6; % P_inlet = P_feed + P_hydrostatic9 (Rho_w * g * h)   
    v_in(1, n+1) = (Qfad / ( pi * (D_sup/2)^2)); % Inlet velocity through compressor 
    v_out(1, n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));% Outlet velocity through nozzle  
    rho(1, n+1) = rho_initial * (T_inlet / T_FAD) * (P_FAD / P(1, n+1)); % Inlet density
    %fprintf(' rho(1, n+1): %.3f m/s\n',  rho(1, n+1));
    v(1,n+1) = (Qfad/A_total); %velocity of sysytem( Hose + Nozzle)
    m_in(1, n+1) = rho(1, n+1) * A_hose * v_in(1, n+1); % mass flow rate entering hose
    m_out(1, n+1) = rho(1, n+1) *(A_nozzle /A_total) * v_out(1, n+1); % mass flow rate exiting through nozzle  
    %In the time integration loop, you first calculate the boundary conditions for the first control volume (inlet) and the first time step (n+1). 
    % Then, you use those values as initial conditions for the rest of the control volumes in the loop (i = 2:N-1) and iterate over all time steps (n = 1:Nt-1).
    
    %% SECTION 6:Update density using the continuity equation (Backward Euler)
    for i = 2:N-1
        %Since you've already provided explicit conditions for the first and last control volumes,inlet and outlet BC 
        % you start the loop at i = 2 and go up to i = N-1 to update the remaining control volumes in between. 
        % This allows you to capture the variations within the hose due to the complex interactions between neighboring control volumes.
        
        % Calculate the mass flow rates entering and exiting the control volume
        m_in(i,n+1) = rho(i, n+1) * A_hose * v_in(i, n+1);
        fprintf(' v_in(i, n+1: %.2f m/s\n',  v_in(i, n+1));
        %Here's how the steps would work:
        %At the beginning of the time integration loop:
        %Calculate v_in(1, n+1) based on the given inlet boundary condition.
        %Inside the loop for i = 2:N-1:
        %Calculate v_in(2, n+1) using the value of v_in(1, n+1) that was calculated earlier in the time integration loop.
        %Calculate m_in(2) using the calculated value of v_in(2, n+1) and the updated value of rho(2, n+1) (which is calculated based on the continuity equation).
        m_out(i,n+1) = rho(i, n+1) * (A_nozzle / A_total) * v_out(i, n+1);
        %m_in(i) and m_out(i) are calculated separately to represent the mass flow rates entering and exiting the control volume at grid point i, respectively. 
        % These terms are essential for maintaining mass conservation in the simulation and ensuring the accuracy of the continuity equation calculations.
        
        %% Contunity Equation loop: The first run calculates the predicted density, 
        % and the second run updates the density using the calculated mass flow rates, velocities, and other parameters. 
        % This iterative process ensures that the density values are updated based on the interactions between neighboring control volumes and the physical properties of the system.
        
        rho(i, n+1) = - ((rho(i+1, n+1) * A_total * (v_in(i+1, n+1))^2) - (rho(i-1, n) * A_total * (v_in(i-1, n))^2)) ...
            / (A_total * v_in(i, n+1)) * (dt / (2 * delta_x(i))) ...
            + dt * ((m_in(i, n+1) - m_out(i, n+1)) / (A_total * v_in(i, n+1))) ...
            + (rho(i, n) * v_in(i, n) * A_total ) / (A_total * v_in(i, n+1));
        
    end
    fprintf(' rho(i, n+1: %.2f m/s\n',  rho(i, n+1));
    
    %% SECTION 7: CFL CALCULATION
    % Calculate characteristic velocity at this time step (use the maximum velocity in the flow field)
    v_max = max(max(v(:, n+1))); % Assuming v is the velocity array at all grid points at time step n+1
    % Calculate CFL number for this time step
    CFL_number = v_max * dt / min(delta_x);
    % Check if CFL number exceeds the maximum allowable value
    if CFL_number > CFL_max
        % Reduce the time step to satisfy the CFL condition
        dt = CFL_max * min(delta_x) / v_max;
        fprintf('Time step reduced at time step %d to satisfy CFL condition.\n', n+1);
    end
    
    %% SECTION 8: Update velocity using the momentum equation (Backward Euler)
    for i = 2:N-1 % Loop for Spatial Grid Point
        m_n(i) = rho(i, n+1) * (A_nozzle / A_total) * (v_out(i, n+1)-v_in(i,n+1)); 
        
        % Calculate Reynolds number and friction factor
        Re = rho(i, n+1) * v(i, n+1) * D_sup / mu;
        %v(i, n+1) is the predicted velocity for the next time step (n+1).
        %% how v(i,n+1) is calculated ? 
        %In summary, for control volume 2 (i = 2), the predicted velocity v(2, n+1) is first calculated using the velocity update equation.
        % This predicted velocity is then used in the calculation of the Reynolds number Re. The Reynolds number Re and other terms are used in the momentum equation to calculate the updated velocity v(2, n+1),
        % which represents the final velocity for control volume 2 at the next time step
        f = (0.79 * log(Re) - 1.64)^(-2);
        % Calculate wall friction w
        w = (f * rho(i, n+1) * v(i, n+1)^2) / (8 * (D_sup/2));
        
        %% Momentum Equation : The momentum equation is indeed run twice within the same loop iteration. 
        % The first run calculates the predicted velocity, and the second run updates the velocity using the calculated Reynolds number and other parameters.
        % This iterative process ensures that the velocity values are updated based on the complex interactions between neighboring control volumes and the various physical properties of the system.
        
        % Update the velocity using the discretized momentum equation (Backward Euler)
        v(i, n+1) = sqrt(((rho(i, n) * A_total * v(i, n)^2) / dt - ...
            ((rho(i+1, n+1) * A_total * v(i+1, n+1)^3 + P(i+1,n+1) * A_total) - ...
            (rho(i-1, n) * A_total * v(i-1, n)^3 + P(i-1,n) * A_total) / (2 * delta_x(i)) + ...
            (w / A_total) - rho(i, n) * g * sin(theta) + ...
            (A_nozzle / A_total) * m_n(i) * (v_out(i, n) - v(i, n)) - ...
            ((P(i+1, n+1) * A_nozzle) - (P(i-1, n) * A_nozzle)) /(2 * delta_x(i)) ...
            / (rho(i, n+1) * A_total)/dt)));
    end
    fprintf(' v(i, n+1: %.2f m/s\n',  v(i, n+1));
    
    %% SECTION 9: PRESSURE CALCULATION
    
    % Update pressure using P = rho * R * T
    for i = 2:N-1
        P(i, n+1) = rho(i, n+1) * R *T(i, n+1); % T is not constant
        % the pressure at each spatial control volume (i) and each time step (n+1)
        % is calculated based on the corresponding density (rho(i, n+1)) and
        % temperature (T(i, n+1)) at that grid point and time step.
        % This will account for the changing temperature in the simulation.
    end
    
    %% SECTION 10: Outlet Boundary Conditions at nozzle  
    P(N, n+1) = 402957.5; % P_outlet = P_atm + Rho_w * g * H
    rho(N, n+1) = rho(N-1, n+1); % Density at the outlet is set to the density at the neighboring grid point
    v(N, n+1) = v(N-1, n+1); % Set outlet velocity based on previous control volume
    m_out(N, n+1) = m_out(N-1, n+1); % Set outlet mass flow rate based on previous control volume
    v_out(N, n+1) = v_out(N-1, n+1); % Set outlet velocity through nozzle based on previous control volume  
end

Answer 1

[Voss]

You can put the following immediately after the line where rho(i,n+1) is calculated:

if isnan(rho(i,n+1))
    keyboard
end

Then when you run the program, if rho(i,n+1) is NaN, the program will pause and give you a chance to inspect the value of the relevant variables on the command line or in the editor.

When you do that, you'll find that m_in(i, n+1) is zero, m_out(i, n+1) is zero, A_total is non-zero and finite, and v_in(i, n+1) is zero, so the expression

((m_in(i, n+1) - m_out(i, n+1)) / (A_total * v_in(i, n+1)))

is 0/0 which evaluatates to NaN. (There are other divide-by-zero situations going on in that line, but the others involve non-zero numerators so they evaluate to Inf.) That's where the NaN is coming from in rho(i,n+1).

I imagine you'll find similar things going on in the calculation of v(i,n+1).

[Walter Roberson]

v_in =zeros(total_control_volumes, Nt);

Initialized to all zero

v_in(:,1)=v_in_initial* ones(total_control_volumes, 1);

does not change any column after column 1

    v_in(1, n+1) = (Qfad / ( pi * (D_sup/2)^2)); % Inlet velocity through compressor 

This is inside the for n loop but not yet the for i loop. It changes column n+1 but only of the first row

   for i = 2:N-1

i starts at 2

     m_in(i,n+1) = rho(i, n+1) * A_hose * v_in(i, n+1);

Nothing has been assigned to v_in(2, n+1) other than the initial zeros, so m_in(i,n+1) is going to be 0

     rho(i, n+1) = - ((rho(i+1, n+1) * A_total * (v_in(i+1, n+1))^2) - (rho(i-1, n) * A_total * (v_in(i-1, n))^2)) ...
            / (A_total * v_in(i, n+1)) * (dt / (2 * delta_x(i))) ...
            + dt * ((m_in(i, n+1) - m_out(i, n+1)) / (A_total * v_in(i, n+1))) ...
            + (rho(i, n) * v_in(i, n) * A_total ) / (A_total * v_in(i, n+1));

On the second line, v_in(i,n+1) appears in the denominator, but it has not been assigned any value other than the initial zero, so (A_total * v_in(i, n+1)) is going to be 0, and you have a division by 0.

[MAYUR MARUTI DHIRDE]

%% SECTION 5: TIME INTEGRATION LOOP
for n = 1:Nt
    rho(1, n+1) = rho_initial * (T_inlet / T_FAD) * (P_FAD / P(1, n+1));
    %% SECTION 6:Update density using the continuity equation (Backward Euler)
    for i = 2:N-1
        for i = 2:N-1
            v_in(i,n+1) = ((Qfad / ( pi * (D_sup/2)^2)));
            v_out(i,n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));
            % Calculate the mass flow rates entering and exiting the control volume
            m_in(i,n+1) = rho(i, n+1) * A_hose * v_in(i, n+1);

In the time integration loop, when calculating m_in(1, n+1) for the first control volume, I used the inlet boundary condition directly and not relying on the continuity equation since it's the inlet and we have those values available.

However, for the subsequent control volumes (2 to N-1), I need to calculate m_in(i, n+1) based on the updated density rho(i, n+1) which will be calculate later using the continuity equation. This creates a dependency loop where the density depends on m_in and m_in depends on the density. How should i resolve this problem ?

[Torsten]

Which discretization scheme in time do you use for your PDEs ?

[MAYUR MARUTI DHIRDE]

Implicit Finite Difference Method: The implicit approach involves discretizing the time-dependent part of the PDE using the backward Euler method, and the central difference scheme is used to discretize the spatial derivatives in the PDE.

[Torsten]

I suspected that you use the Implicit Euler. But then, it's usual that the unknown values appear on the left and on the right-hand sides of the equations. Thus you either have to iterate to get the value for the new time stamp or you have to use a linear or nonlinear solver to solve them.

[MAYUR MARUTI DHIRDE]

I have tried solving fsolve approach :

%% Section 1: PARAMETERS

L = 400; % Length of the hose (m)

N = 100; % Number of control volumes

% Operation and Design Properties

Qfad =4.31; % Air flow rate in m^3/s ; Q= 259 m^3/min = 4.31

D_sup = 104e-3; % Diameter of supply hose & riser in meters

D_nozzle= 2e-3; % Diameter of nozzle

s=0.3; % Nozzle spacing

C_d= 0.55; % Nozzle discharge coefficient

P_FAD = 101325; % Pressure

T_FAD = 293.15; % Temperature in Kelvin

h = 30; % Water depth in meters

Rho_w = 1025; % Density of water in kg/m^3

g = 9.81; % Gravity acceleration in m/s^2

P_atm = 1.013e5; % Atmospheric pressure in Pa

T_inlet = 291.15; % Air temperature in Celsius

R = 287; % Specific gas constant in J/(kg·K)

mu = 18e-6; % Dynamic viscosity of air in Pa·s

T_FAD_ref = 20; % Reference temperature at FAD condition in Celsius

P_FAD_ref = 1.013e5; % Reference pressure at FAD condition in Pa

P_Fed= 1030000; % 10.3 Bar to Pascals

theta = pi; %Angle of the hose with respect to the horizontal (radians)

A_hose = pi * (D_sup/2)^2;% Area of the hose (m^2)

A_nozzle= pi * (D_nozzle/2)^2;% Area of the nozzle (m^2)

CFL_max = 0.5; % Define maximum allowable CFL number

%% SECTION 2:SPATIAL DISCRETIZATION OF THE CONTROL VOLUMES

% spatial discretization to account for the presence of nozzles in the system.

num_nozzles = floor(L / s); % Calculate the total number of nozzles

total_control_volumes = N + num_nozzles; % Total number of control volumes including nozzles

% Updating Control Volume Boundaries:

x = linspace(0, L, total_control_volumes + 1); % Update the array of control volume boundaries

%Recalculating Control Volume Lengths:

delta_x = x(2:end) - x(1:end-1);

%% SECTION 3: INITIALIZING ARRAYS

% Number of time steps

Nt = 10000;

% Time step size

dt = 0.001; % Example time step size (you can adjust it based on stability requirements)

% Time array

time = linspace(0, Nt*dt, Nt+1); % linspace(start, end, N);Nt*dt=total time duration of the simulation.

% Initialize arrays to store the variables at each grid point and time step

rho = zeros(total_control_volumes, Nt+1); %Nt+1= ensure that you have enough space to store density values for all time steps from 0 to Nt.

v = zeros(total_control_volumes, Nt+1);

T = zeros(total_control_volumes, Nt+1);

P = zeros(total_control_volumes, Nt+1);

m_in=zeros(total_control_volumes, Nt+1);

m_out=zeros(total_control_volumes, Nt+1);

v_in =zeros(total_control_volumes, Nt+1);

v_out=zeros(total_control_volumes, Nt+1);

m_n= zeros(total_control_volumes, Nt+1);

%% SECTION 4: INITIAL CONDITION

%These initial conditions are used as starting values for the simulation.

rho_initial = 1.225 ; % Initial density (kg/m^3) (rho= P/R*T)=(1.013e5/287* 293.15)=1.225kg/m^3

P_initial = 1030000; % Initial pressure (Pa)

T_inital= 293.15; % Initial Temperature(K)

v_initial=(Qfad / A_total); %velocity of sysytem( Hose + Nozzle)

v_in_initial=(Qfad / A_hose); %velocity of hose

v_out_initial= (Qfad / (C_d * pi * (D_nozzle/2)^2));%velocity of nozzle

m_initial=rho_initial * A_total*v_in_initial;

m_out_initial= rho_initial * (A_nozzle / A_total) * v_out_initial;

m_n_initial = rho_initial* (A_nozzle / A_total) * (v_out_initial-v_in_initial);

% Set initial conditions for each variable at each grid point

% array(row_index, column_index)

rho(:, 1) = rho_initial * ones(total_control_volumes, 1);

% So, rho(:, 1) means all rows (all control volumes) of the first column (first time step) in the rho matrix.

%initializing the density values for all control volumes at the initial time step (time step n = 1) in your simulation

P(:, 1) = P_initial * ones(total_control_volumes, 1);

T(:, 1) =T_inital * ones(total_control_volumes, 1);

v(:,1)=v_initial * ones(total_control_volumes, 1);

v_in(:,1)= v_in_initial* ones(total_control_volumes, 1);

v_out(:,1)=v_out_initial* ones(total_control_volumes, 1);

m_in(:,1)=m_initial* ones(total_control_volumes, 1);

m_out(:,1)=v_out_initial* ones(total_control_volumes, 1);

m_n(:,1)=m_n_initial* ones(total_control_volumes, 1);

%% SECTION 5: TIME INTEGRATION LOOP

for n = 1:Nt %Loop for Time Steps [iterates from the first time step to last time step]

% Inlet Boundary Conditions

P(1, n+1) = P_Fed + (Rho_w*g*h); % P_inlet = P_feed + P_hydrostatic9 (Rho_w * g * h)=

v_in(1, n+1) = (Qfad / ( pi * (D_sup/2)^2)); % Inlet velocity through compressor

v_out(1, n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));% Outlet velocity through nozzle

rho(1, n+1) = rho_initial * (T_inlet / T_FAD) * (P_Fed / P(1, n+1)); % Inlet density

v(1,n+1) = (Qfad/A_total); %velocity of sysytem( Hose + Nozzle)

m_in(1, n+1) = rho(1, n+1) * A_hose * v_in(1, n+1); % mass flow rate entering hose

m_out(1, n+1) = rho(1, n+1) *(A_nozzle /A_total) * v_out(1, n+1); % mass flow rate exiting through nozzle

% Loop over spatial grid points

for i = 2:total_control_volumes-1

v_in(i,n+1) = ((Qfad / ( pi * (D_sup/2)^2)));

v_out(i,n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));

% Update velocity based on average of v_in and v_out

v(i, n+1) = 0.5 * (v_in(i, n+1) + v_out(i, n+1));

% Calculate the mass flow rates entering and exiting the control volume

m_in(i,n+1) = rho(i, n+1) * A_hose * v_in(i, n+1);

m_out(i,n+1) = rho(i, n+1) * (A_nozzle / A_total) * v_out(i, n+1);

% Calculate initial guesses for density and velocity

rho_guess = rho(:, n); % Use the density values from the current time step

%used as the starting point for the iterative solver (fsolve) to find a density value that satisfies the density equation.%rho is a matrix with dimensions (total_control_volumes, Nt+1), and the column n corresponds to the density values at the current time step.

%% Procedure

%The initial guess rho_guess is calculated based on the density values from the current time step.

%The fsolve function is used to find a new density value rho_new that satisfies the density equation, starting from the initial guess.

%The calculated rho_new is then used to update the density array for the next time step, ensuring that the density equation is satisfied at the current control volume and time step.

% Call fsolve to solve the density equation

options = optimoptions('fsolve', 'Display', 'off'); %used to configure the behavior of the optimization solver, specifically the fsolve function in this case.

% optimoptions function allows you to set various options that control how the optimization algorithm behaves during its iterations.

rho_new = fsolve(@(rho_guess) densityEquation(rho_guess,A_total, v, n, i, delta_x, dt, m_in, m_out), rho_guess, options);

%Here, fsolve is called to solve the equation defined by the anonymous function @(rho_guess) densityEquation(...). This equation aims to find the values of rho_guess that satisfy the specified equation, where the equation is defined in the densityEquation function.

% The function densityEquation calculates the residuals of the density equation based on the input arguments.

% Update the density at control volume i and time step n+1

rho(i, n+1) = rho_new(i);

%Finally, the calculated values of rho_new are used to update the density values at control volume i and time step n+1.

% This step ensures that you are using the corrected density values for the next time step in your simulation.

end

% SECTION 8: Update velocity using the momentum equation (Backward Euler)

for i = 2:total_control_volumes-1

% Calculate the momentum transferred between control volumes

m_n(i,n+1) = rho(i, n+1) * (A_nozzle / A_total) * (v_out(i, n+1) - v_in(i,n+1));

% Calculate Reynolds number and friction factor

Re = rho(i, n+1) * v(i, n+1) * D_sup / mu;

f = (0.79 * log(Re) - 1.64)^(-2);

% Calculate wall friction w

w = (f * rho(i, n+1) * v(i, n+1)^2) / (8 * (D_sup/2));

% Calculate initial guess for velocity

v_guess = v(:, n); % Use the velocity values from the current time step (n), not (n+1)

%fprintf('v_guess : %.2f kg/m^3\n', v_guess );

% Call fsolve to solve the momentum equation

options = optimoptions('fsolve', 'Display', 'off');

v_new = fsolve(@(v_guess) momentumEquation(v_guess, A_total, n, i, delta_x, dt, m_n, w, g, theta), v_guess, options);

% Update the velocity at control volume i and time step n+1

v(i, n+1) = v_new;

fprintf('v(i, n+1) : %.1f\n', v(i, n+1));

end

% Update pressure using P = rho * R * T

for i = 2:total_control_volumes-1

P(i, n+1) = rho(i, n+1) * R * T(i, n+1); % T is not constant

% the pressure at each spatial control volume (i) and each time step (n+1)

% is calculated based on the corresponding density (rho(i, n+1)) and

% temperature (T(i, n+1)) at that grid point and time step.

% This will account for the changing temperature in the simulation.

end

P(total_control_volumes, n+1) = 402957.5; % P_outlet = P_atm + Rho_w * g * H

rho(total_control_volumes, n+1) = rho(total_control_volumes-1, n+1); % Density at the outlet is set to the density at the neighboring grid point

v(total_control_volumes, n+1) = v(total_control_volumes-1, n+1); % Set outlet velocity based on previous control volume

m_out(total_control_volumes, n+1) = m_out(total_control_volumes-1, n+1); % Set outlet mass flow rate based on previous control volume

v_out(total_control_volumes, n+1) = v_out(total_control_volumes-1, n+1); % Set outlet velocity through nozzle based on previous control volume

end

% Print the maximum and minimum values of the velocity,Pressure and Density array

%fprintf('Max Density: %.2f kg/m^3\n', max(rho(:)));

%fprintf('Min Desnity: %.2f kg/m^3\n', min(rho(:)));

function residuals = densityEquation(rho_guess,A_total,v, n, i, delta_x, dt, m_in, m_out)

% Calculate the residuals of the density equation

%densityEquation function provides the equation you want to solve for density, and fsolve iteratively adjusts the values of rho_guess to find a solution that satisfies this equation.

% Calculate the terms in your equation based on rho_guess, v, n, and i

term1 = - (((rho_guess(i+1) * A_total * (v(i+1, n+1))^2) - (rho_guess(i-1) * A_total * (v(i-1, n+1))^2)) * dt) ...

/ (A_total * v(i, n+1)) * (2 * delta_x(i));

term2 = (dt * (m_in(i, n+1) - m_out(i, n+1)) / (A_total * v(i, n+1)));

term3 = (rho_guess(i) * v(i, n) * A_total) / (A_total * v(i, n+1));

LHS = term1 + term2 + term3;

RHS = LHS; % The right-hand side is the same as the left-hand side in this case

% Calculate the residuals (LHS - RHS)

residuals = LHS - RHS;

% calculates the difference between the left-hand side and the right-hand side of the density equation. This difference represents the residual of the equation.

% If the residual is close to zero, it indicates that the equation is satisfied and the density guess is reasonable.

% If the residual is not close to zero, it indicates that the guessed density needs to be adjusted to satisfy the equation.

end

function residuals = momentumEquation(v_guess,A_total, n, i, delta_x, dt, m_n, w, g, theta)

% Calculate the residuals of the momentum equation

% Calculate the terms in your equation based on v_guess, n, i, and other parameters

term1 = rho(i, n+1) * A_total * v_guess^2 / dt; % You may need to include "rho" here

term2 = (((rho(i+1, n+1) * A_total * v_guess(i+1, n+1)^2) + (P(i+1, n+1) * A_total)) - ...

((rho(i-1, n+1) * A_total * v_guess(i-1, n+1)^2) + (P(i-1, n+1) * A_total)) / (2 * delta_x(i)));

term3 = w / A_total;

term4 = rho(i, n+1) * g * sin(theta);

term5 = (A_nozzle / A_total) * m_n(i, n+1) * (v_out(i, n) - v_guess(i, n));

term6 = ((P(i+1, n+1) * A_nozzle) - (P(i-1, n+1) * A_nozzle)) / (2 * delta_x(i));

term7 = rho(i, n+1) * A_total / dt;

% Calculate final term

final_term = (term1 - term2 + term3 - term4 + term5 - term6) / term7;

% Update the velocity using the calculated final term

LHS = v_guess - sqrt(final_term);

% Define the LHS and RHS of the momentum equation

RHS = 0; % This is for backward Euler

% Calculate the residuals (LHS - RHS)

residuals = LHS - RHS;

end

%%%%ERROR

The above code stops running at this line rho_new = fsolve(@(rho_guess) densityEquation(rho_guess,A_total, v, n, i, delta_x, dt, m_in, m_out), rho_guess, options);

And open this fsolve.m file

[Torsten]

% Calculate the mass flow rates entering and exiting the control volume
m_in(i,n+1) = rho(i, n+1) * A_hose * v_in(i, n+1);
m_out(i,n+1) = rho(i, n+1) * (A_nozzle / A_total) * v_out(i, n+1);

I just picked these two lines of your code. rho(i,n+1) is not known !

Start with Euler implicit in time for an ordinary differential equation, e.g.

y'(t) = y^2, y(1) = 1

to learn how to proceed to couple your discretization with the solution of the nonlinear systems in each time step.

[MAYUR MARUTI DHIRDE]

%% Section 1:  PARAMETERS 

L = 400; % Length of the hose (m)
N = 100; % Number of control volumes

% Operation and Design Properties
Qfad =0.33; % Air flow rate in m^3/s ; Q= 259 m^3/min = 4.31
D_sup = 104e-3; % Diameter of supply hose & riser in meters
D_nozzle= 2e-3; % Diameter of nozzle 
s=0.5; % Nozzle spacing 
C_d= 0.55; % Nozzle discharge coefficient 
P_FAD = 101325; % Pressure 
T_FAD = 293.15; % Temperature in Kelvin
h = 30; % Water depth in meters
Rho_w = 1025; % Density of water in kg/m^3
g = 9.81; % Gravity acceleration in m/s^2
P_atm = 1.013e5; % Atmospheric pressure in Pa
T_inlet = 291.15; % Air temperature in Celsius
R = 287; % Specific gas constant in J/(kg·K)
mu = 18e-6; % Dynamic viscosity of air in Pa·s
T_FAD_ref = 20; % Reference temperature at FAD condition in Celsius
P_FAD_ref = 1.013e5; % Reference pressure at FAD condition in Pa
P_Fed= 1030000; % 10.3 Bar to Pascals 
theta = pi; %Angle of the hose with respect to the horizontal (radians)
A_hose = pi * (D_sup/2)^2;% Area of the hose (m^2)
A_nozzle= pi * (D_nozzle/2)^2;% Area of the nozzle (m^2)
CFL_max = 0.5; % Define maximum allowable CFL number
A_total= A_hose+A_nozzle;


%% SECTION 2:SPATIAL DISCRETIZATION OF THE CONTROL VOLUMES
% spatial discretization to account for the presence of nozzles in the system.
num_nozzles = floor(L / s); % Calculate the total number of nozzles
total_control_volumes = N + num_nozzles; % Total number of control volumes including nozzles
% Updating Control Volume Boundaries:
x = linspace(0, L, total_control_volumes + 1); % Update the array of control volume boundaries
%Recalculating Control Volume Lengths:
delta_x = x(2:end) - x(1:end-1); 
%% SECTION 3: INITIALIZING ARRAYS 
% Number of time steps
Nt = 100;
% Time step size
dt = 0.1; % Example time step size (you can adjust it based on stability requirements)  
% Time array
time = linspace(0, Nt*dt, Nt+1); % linspace(start, end, N);Nt*dt=total time duration of the simulation.

% Initialize arrays to store the variables at each grid point and time step
rho = zeros(total_control_volumes, Nt+1); %Nt+1= ensure that you have enough space to store density values for all time steps from 0 to Nt.
v = zeros(total_control_volumes, Nt+1); 
T = zeros(total_control_volumes, Nt+1);
P = zeros(total_control_volumes, Nt+1);
m_in=zeros(total_control_volumes, Nt+1);
m_out=zeros(total_control_volumes, Nt+1);
v_in =zeros(total_control_volumes, Nt+1);
v_out=zeros(total_control_volumes, Nt+1);
m_n= zeros(total_control_volumes, Nt+1); 

%% SECTION 4: INITIAL CONDITION
%These initial conditions are used as starting values for the simulation.
rho_initial = 1.225 ;  % Initial density (kg/m^3) (rho= P/R*T)=(1.013e5/287* 293.15)=1.225kg/m^3

P_initial = 1030000;   % Initial pressure (Pa)
T_inital= 293.15;     % Initial Temperature(K)        
v_initial=(Qfad / A_total); %velocity of sysytem( Hose + Nozzle)
v_in_initial=(Qfad / A_hose); %velocity of hose 
v_out_initial= (Qfad / (C_d * pi * (D_nozzle/2)^2));%velocity of nozzle 
m_initial=rho_initial * A_total*v_in_initial; 
m_out_initial= rho_initial * (A_nozzle / A_total) * v_out_initial;
m_n_initial = rho_initial* (A_nozzle / A_total) * (v_out_initial-v_in_initial);

% Set initial conditions for each variable at each grid point
% array(row_index, column_index)
rho(:, 1) = rho_initial * ones(total_control_volumes, 1);
% So, rho(:, 1) means all rows (all control volumes) of the first column (first time step) in the rho matrix.
%initializing the density values for all control volumes at the initial time step (time step n = 1) in your simulation
P(:, 1) = P_initial * ones(total_control_volumes, 1);
T(:, 1) =T_inital * ones(total_control_volumes, 1);
v(:,1)=v_initial * ones(total_control_volumes, 1);
v_in(:,1)= v_in_initial* ones(total_control_volumes, 1);
v_out(:,1)=v_out_initial* ones(total_control_volumes, 1);
m_in(:,1)=m_initial* ones(total_control_volumes, 1);
m_out(:,1)=v_out_initial* ones(total_control_volumes, 1);
m_n(:,1)=m_n_initial* ones(total_control_volumes, 1);

%% SECTION 5: TIME INTEGRATION LOOP
for n = 1:Nt %Loop for Time Steps [iterates from the first time step to last time step]

     % CFL calculation and adjustment
    v_max = max(max(v(:, n+1)));
    CFL_number = v_max * dt / min(delta_x);
    
    if CFL_number > CFL_max
        dt = CFL_max * min(delta_x) / v_max;
        fprintf('Time step reduced at time step %d to satisfy CFL condition.\n', n+1);
    end
    
    
    % Inlet Boundary Conditions
    P(1, n+1) = P_Fed + (Rho_w*g*h); % P_inlet = P_feed + P_hydrostatic9 (Rho_w * g * h)= 
    v_in(1, n+1) = (Qfad / ( pi * (D_sup/2)^2)); % Inlet velocity through compressor 
    v_out(1, n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));% Outlet velocity through nozzle  
    rho(1, n+1) = rho_initial * (T_inlet / T_FAD) * (P_Fed / P(1, n+1)); % Inlet density
    v(1,n+1) = (Qfad/A_total); %velocity of sysytem( Hose + Nozzle)
    m_in(1, n+1) = rho(1, n+1) * A_hose * v_in(1, n+1); % mass flow rate entering hose
    m_out(1, n+1) = rho(1, n+1) *(A_nozzle /A_total) * v_out(1, n+1); % mass flow rate exiting through nozzle  
   
       % Loop over spatial grid points
    % Loop over spatial grid points
     for i = 2:N-1

        v_in(i,n+1) = ((Qfad / ( pi * (D_sup/2)^2)));
        v_out(i,n+1) = (Qfad / (C_d * pi * (D_nozzle/2)^2));
        % Update velocity based on average of v_in and v_out
        v(i, n+1) = 0.5 * (v_in(i-1, n+1) + v_out(i-1, n+1));

        % Calculate initial guesses for density and velocity
        rho_guess = rho(:, n); % Use the density values from the current time step
        %used as the starting point for the iterative solver (fsolve) to find a density value that satisfies the density equation.
        %rho is a matrix with dimensions (total_control_volumes, Nt+1), and the column n corresponds to the density values at the current time step.
        % Calculate the mass flow rates entering and exiting the control volume
        m_in(i,n+1) = rho_guess(i) * A_hose * v_in(i, n+1);
        m_out(i,n+1)= rho_guess(i) * (A_nozzle / A_total) * v_out(i, n+1);
        %fprintf('m_in(i,n+1): %.1f\n',m_in(i,n+1));
      

       %% Procedure
       %The initial guess rho_guess is calculated based on the density values from the current time step.
       %The fsolve function is used to find a new density value rho_new that satisfies the density equation, starting from the initial guess.
       %The calculated rho_new is then used to update the density array for the next time step, ensuring that the density equation is satisfied at the current control volume and time step.

        % fprintf('rho_guess : %.2f kg/m^3\n', rho_guess);
        % Call fsolve to solve the density equation
        options = optimoptions('fsolve', 'Display', 'off'); %used to configure the behavior of the optimization solver, specifically the fsolve function in this case.
        % optimoptions function allows you to set various options that control how the optimization algorithm behaves during its iterations.
        rho_new = fsolve(@(rho_guess) densityEquation(rho_guess,A_total, v, n, i, delta_x, dt, m_in, m_out), rho_guess, options);
        %Here, fsolve is called to solve the equation defined by the anonymous function @(rho_guess) densityEquation(...). This equation aims to find the values of rho_guess that satisfy the specified equation, where the equation is defined in the densityEquation function. 
        % The function densityEquation calculates the residuals of the density equation based on the input arguments.
        
        % Update the density at control volume i and time step n+1
        rho(i, n+1) = rho_new(i);
       % fprintf('rho(%d, %d) : %.1f\n', i, n+1, rho(i, n+1))
        %Finally, the calculated values of rho_new are used to update the density values at control volume i and time step n+1. 
        % This step ensures that you are using the corrected density values for the next time step in your simulation.
     end
    % SECTION 8: Update velocity using the momentum equation (Backward Euler)
for i = 2:N-1
    % Calculate the momentum transferred between control volumes
   m_n(i,n+1) = rho(i, n+1) * (A_nozzle / A_total) * (v_out(i, n+1) - v_in(i,n+1));
    
    % Calculate initial guess for velocity
    v_guess = v(:, n); % Use the velocity values from the current time step (n), not (n+1)
    %fprintf('v_guess : %.2f kg/m^3\n', v_guess );

    % Calculate Reynolds number and friction factor
    Re = rho(i, n+1) * v_guess(i) * D_sup / mu;
    f = (0.79 * log(Re) - 1.64)^(-2);
    % Calculate wall friction w
    w = (f * rho(i, n+1) * v_guess(i)^2) / (8 * (D_sup/2));

   % Call fsolve to solve the momentum equation
   options = optimoptions('fsolve', 'Display', 'off');
  v_new = fsolve(@(v_guess) momentumEquation(v_guess, rho, P, A_nozzle, v_out, A_total, n, i, delta_x, dt, m_n, w, g, theta), v_guess, options);

   % Update the velocity at control volume i and time step n+1
   v(i, n+1) = v_new(i);
    %fprintf('v(%d, %d) : %.1f\n', i, n+1, v(i, n+1))
    
%end
 %Update pressure using P = rho * R * T
    for i = 2:N-1 
        P(i, n+1) = rho(i, n+1) * R * T(i, n+1); % T is not constant
        % the pressure at each spatial control volume (i) and each time step (n+1)
        % is calculated based on the corresponding density (rho(i, n+1)) and
        % temperature (T(i, n+1)) at that grid point and time step.
        % This will account for the changing temperature in the simulation.
    end
end
P(N, n+1) =P_atm + Rho_w * g * h; % P_outlet = P_atm + Rho_w * g * H
rho(N, n+1) = rho(N-1, n+1); % Density at the outlet is set to the density at the neighboring grid point
v(N, n+1) = v(N-1, n+1); % Set outlet velocity based on previous control volume
m_out(N, n+1) = m_out(N-1, n+1); % Set outlet mass flow rate based on previous control volume
v_out(N, n+1) = v_out(N-1, n+1); % Set outlet velocity through nozzle based on previous control volume

end 
 % Print the maximum and minimum values of the velocity,Pressure and Density  array

fprintf('Max Density: %.2f kg/m^3\n', max(rho(:)));
fprintf('Min Desnity: %.2f kg/m^3\n', min(rho(:)));

fprintf('Max Velocity: %.2f m/s\n', max(v(:))); %%f with .2 precision specifier would round the number to two decimal places. For example, 3.14159265 would be displayed as 3.14.
fprintf('Min Velocity: %.2f m/s\n', min(v(:)));

fprintf('Max Pressure: %.2f Pa\n', max(P(:)));
fprintf('Min Pressure: %.2f Pa\n', min(P(:)));

%fprintf('Size of x_center: %dx%d\n', size(x));
%fprintf('Size of rho: %dx%d\n', size(rho));

% Plotting the results
figure;
plot(x(1:total_control_volumes), rho);
xlabel('Length(m)');
ylabel('Density (kg/m^3)');
title('Density vs. Length');

% Figure 2: Velocity vs. Length
figure;
plot(x(1:total_control_volumes), v);
xlabel('Length (m)');
ylabel('Velocity (m/s)');
title('Velocity vs. Length');


% Figure 3: Pressure vs. Length
figure;
plot(x(1:total_control_volumes), P);
xlabel('Length (m) ');
ylabel('Pressure (Pascals)');
title('Pressure vs. Length');

function residuals = densityEquation(rho_guess, A_total, v, n, i, delta_x, dt, m_in, m_out)
    term1 = - (((rho_guess(i+1) * A_total * (v(i+1, n+1))^2) - (rho_guess(i-1) * A_total * (v(i-1, n+1))^2)) * dt) ...
                / (A_total * v(i, n+1)) * (2 * delta_x(i));
    term2 = (dt * (m_in(i, n+1) - m_out(i, n+1)) / (A_total * v(i, n+1)));
    term3 = (rho_guess(i) * v(i, n) * A_total) / (A_total * v(i, n+1));

    LHS = term1 + term2 + term3;
    RHS = LHS; % The right-hand side is the same as the left-hand side in this case

    % Calculate the residuals (LHS - RHS)
    residuals = LHS - RHS;
end


function residuals = momentumEquation(v_guess, rho, P, A_nozzle, v_out, A_total, n, i, delta_x, dt, m_n, w, g, theta)
    % Calculate the residuals of the momentum equation
   
    % Calculate the terms in your equation based on v_guess, rho, P, A_nozzle, v_out, A_total, n, i, and other parameters
   
    Term1 = rho(i, n) * A_total * v_guess(i)^2 / dt;
    Term2 = (((rho(i+1, n+1) * A_total * v_guess(i+1)^3) + (P(i+1, n+1) * A_total)) - ...
             ((rho(i-1, n+1) * A_total * v_guess(i-1)^3) + (P(i-1, n+1) * A_total)) / (2 * delta_x(i)));
    Term3 = w / A_total;
    Term4 = rho(i, n+1) * g * sin(theta);
    Term5 = (A_nozzle / A_total) * m_n(i, n+1) * (v_out(i, n) - v_guess(i));
    Term6 = ((P(i+1, n+1) * A_nozzle) - (P(i-1, n+1) * A_nozzle)) / (2 * delta_x(i));
    Term7 = rho(i, n+1) * A_total / dt;

    % Calculate final term
    final_term = (Term1 - Term2 + Term3 - Term4 + Term5 - Term6) / Term7;

    % Update the velocity using the calculated final term
    LHS = v_guess - sqrt(final_term);
    
    % Define the LHS and RHS of the momentum equation
    RHS = 0; % This is for backward Euler

    % Calculate the residuals (LHS - RHS)
    residuals = LHS - RHS;
end

why this code take too much time to run ? I tried reducing value of N, Nt and dt but still facing same issue. Please help!

[Torsten]

If you solve momentum and continuity equation separately as you do, you have to iterate to get a converged solution - even if you use a nonlinear solver like fsolve in between. I miss this iteration part in your code.

My advice:

Make a function f(t,y) where y is the vector of all your solution variables y = (density(1:nx), velocity(1:nx),...). Then call fsolve in a loop for the function F(y) = (y-y_(n))/dt - f(t_(n+1),y) = 0 where y_n is the solution of the last time step. This call should give you the solution y = y_(n+1).