清空筛选条件
清空筛选条件

Why do I get NaN error in A_momentum and (-F)?

6 次查看（过去 30 天）
MAYUR MARUTI DHIRDE
评论： MAYUR MARUTI DHIRDE 2023-9-17
Ran in:
; Momentum Equation.
Problem: In this discretize form v_i, v_i+1 and v_i-1are unknow. Which make equation non linear, To solve this non linear equation, i used iterative Newton Raphason Method, which uses initial guess value to solve. The matrix equation were formed to obtain the unknown values. But i am getting error at v(i) , v(i+1) , v(i-1) as NaN , because they are not getting updated using iterative solver.
% Coefficients based on the discretized momentum equation
A_momentum(i, i) = rho_solution(i) * v(i) / (2 * dx);
A_momentum(i, i + 1) = -rho_solution(i) * v(i + 1) / (2 * dx);
A_momentum(i, i - 1) = -rho_solution(i) * v(i - 1) / (2 * dx);
Below is the complete code :
% Parameters
L = 400; % Length of the domain
Nx = 50; % Number of spatial grid points
h = 30;
P_atm = 1.013e5; % Atmospheric pressure in Pa;
R = 287; % Specific gas constant
T_FAD = 293.15; % Reference temperature (e.g., room temperature)
Qfad = 4.31; % Flow rate
D_sup = 104e-3; % Diameter of the supply hose
T_inlet = 291.15; % Inlet temperature
P_FAD = 101325;
dx = L / (Nx-1); % Spatial step size
P_Fed = 1.33e6;
% Define spatial grid
x = linspace(0, L, Nx); % Spatial points
% Initialize variables
v = zeros(Nx, 1); % Initialize velocity field
P = zeros(Nx, 1); % Initialize pressure field
rho = zeros(Nx, 1); % Initialize density field
% Initialize the coefficient matrix and RHS vector for continuity equation
A_cont = zeros(Nx, Nx);
b_cont = zeros(Nx, 1);
A_momentum = zeros(Nx, Nx);
F = zeros(Nx, 1);
% Initialize initial conditions
rho_initial = P_atm / (R * T_FAD);
P_initial = P_atm;
% Set initial conditions
rho(:) = rho_initial;
P(:) = P_initial;
% Set initial guess for velocity
v_guess = ones(Nx, 1); % Initial guess for velocity field
% Convergence tolerance
tolerance = 1e-6;
% Maximum number of iterations
max_iterations = 100;
% Initialize b_cont based on initial conditions and boundary condition
b_cont(1) = rho_initial * (T_inlet / T_FAD) * (P_FAD / P_initial); % Use a suitable value based on your problem's boundary condition
b_cont(2:Nx-1) = rho_initial; % Assuming no source terms in the interior
b_cont(Nx) = rho(Nx-1); % Use a suitable value based on your problem's boundary condition
% Boundary condition of first grid point
P(1) = 1.33e6;
rho(1) = rho_initial * (T_inlet / T_FAD) * (P_FAD / P_initial);
v(1) = Qfad / (pi * (D_sup/2)^2);
% Avoid singularities (small regularization )
A_cont = A_cont + 1e-6 * eye(Nx);
A_momentum = A_momentum + 1e-6 * eye(Nx);
for i = 2:Nx-1
% Fill the diagonal and off-diagonal elements of A_cont
A_cont(i, i) = 0; % Diagonal elements (zero for central difference)
A_cont(i, i+1) = 1 / (2 * dx); % Right neighbor
A_cont(i, i-1) = -1 / (2 * dx); % Left neighbor
end
rho_solution = A_cont \ b_cont;
% Newton-Raphson iteration for the momentum equation
for iter = 1:max_iterations
% Fill the coefficient matrix A_momentum and right-hand side vector b_momentum
for i = 2:Nx - 1
% Coefficients based on the discretized momentum equation
A_momentum(i, i) = rho_solution(i) * v(i) / (2 * dx);
A_momentum(i, i + 1) = -rho_solution(i) * v(i + 1) / (2 * dx);
A_momentum(i, i - 1) = -rho_solution(i) * v(i - 1) / (2 * dx);
% Calculate pressure using the ideal gas law
P(i) = rho_solution(i) * R * T_FAD;
% Compute the nonlinear function F(v)
F(i) = rho_solution(i) * v(i) * (v(i + 1) - v(i - 1)) / (2 * dx) + ...
(P(i + 1) - P(i - 1)) / (2 * dx);
end
%%%EDIT%%%
%%%This line was below the fprintf commands
% Solve the linearized system A_momentum * delta_v = -F for delta_v
delta_v = A_momentum \ (-F);
% Print values for debugging
fprintf('Iteration %d:\n', iter);
fprintf('delta_v: %.6f\n', delta_v);
fprintf('v: %.6f\n', v);
fprintf('F: %.6f\n', F);
% Update the velocity field v
v = v + delta_v;
% Check for convergence
if max(abs(delta_v)) < tolerance
fprintf('Converged after %d iterations.\n', iter);
break;
end
end
Warning: Matrix is singular to working precision.
Iteration 1:
delta_v: -0.000000 delta_v: Inf delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN
v: 507.365240 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000 v: 0.000000
F: 0.000000 F: -75257.875000 F: -6202187595.375010 F: -6163910180.442371 F: -6202288895.375009 F: -6164011480.442371 F: -6202390195.375009 F: -6164112780.442370 F: -6202491495.375009 F: -6164214080.442370 F: -6202592795.375008 F: -6164315380.442370 F: -6202694095.375009 F: -6164416680.442369 F: -6202795395.375009 F: -6164517980.442369 F: -6202896695.375008 F: -6164619280.442369 F: -6202997995.375006 F: -6164720580.442369 F: -6203099295.375007 F: -6164821880.442369 F: -6203200595.375006 F: -6164923180.442369 F: -6203301895.375005 F: -6165024480.442368 F: -6203403195.375006 F: -6165125780.442367 F: -6203504495.375004 F: -6165227080.442366 F: -6203605795.375004 F: -6165328380.442366 F: -6203707095.375003 F: -6165429680.442366 F: -6203808395.375003 F: -6165530980.442365 F: -6203909695.375003 F: -6165632280.442363 F: -6204010995.375002 F: -6165733580.442364 F: -6204112295.375001 F: -6165834880.442363 F: -6204213595.375001 F: -6165936180.442362 F: -6204314895.375001 F: -6166037480.442363 F: -6204416195.375001 F: -6166138780.442363 F: -6204517495.375001 F: 0.000000
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
Iteration 2:
delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: NaN delta_v: -0.000000
v: 507.365240 v: Inf v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN v: NaN
F: 0.000000 F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: NaN F: 0.000000
Converged after 2 iterations.
% Boundary condition for last grid point
P(Nx) = 402957.5; % P_atm + Rho_w * g * h
rho(Nx) = rho(Nx - 1);
v(Nx) = v(Nx - 1);
% Plot density profile
%figure;
%plot(x, rho_solution, 'b-', 'LineWidth', 2);
%xlabel('Spatial Position (x)');
%ylabel('Density (rho)');
%title('Density Profile');
%grid on;
% Plot velocity profile
%figure;
%plot(x, v, 'r-', 'LineWidth', 2);
%xlabel('Spatial Position (x)');
%ylabel('Velocity (v)');
%title('Velocity Profile');
%grid on;

请先登录，再回答此问题。

回答（1 个）

Yash
Yash 2023-9-13
Hi Mayur,
I understand that you are interested in using the Newton Raphson Method to solve the equations. This method involves making an initial guess and then iteratively converging to a solution by minimizing the error.
Upon reviewing your code, I have identified that it produces an inconsistent set of equations, resulting in "NaN" values in the output. Here are the points that led me to this conclusion:
  1. The variable "v" is initialized as a 50*1 vector of zeros, and only the first element "v(1)" is updated with a specific value. The remaining elements of "v" remain as zeros.
  2. Consider the first iteration, only the element "A_momentum(i, i-1)" for "i=2" is non-zero, as it involves multiplication with "v(i-1)". All other elements of "A_momentum" are zeros.
  3. The vector "F" has non-zero elements for most of its entries, and it has a size of 50x1.
  4. Considering that "A" is a 50x50 matrix and "B" is a 50x1 column vector with most elements non-zero, the equation "A*x = B" will have no solution if A has only one non-zero element.
The reason for this inconsistency, as I can understand, is that you have defined an initial guess "v_guess" using the below line, but you haven't assigned this value to "v".
v_guess = ones(Nx, 1); % Initial guess for velocity field
By assigning "v_guess" to "v", you will now notice that there are multiple non-zero values in "A_momentum", and the "NaN" values are absent in the first few iterations.
However, note that even after properly initializing "v", you still need to review your equations because the solution is still not converging as the value of "delta_v" keeps on increasing with each iteration.
I hope this helps you in resolving the issue.

类别

在 Help Center 和 File Exchange 中查找有关 Newton-Raphson Method 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by