newton raphson for nonlinear heat equation scheme
15 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I am trying to apply newton raphson method to linearize a scheme arise from non linear heat equation of the form
and the scheme is :
where are non linear term that was splitted to convex-concave function.
I have written this code but I couldn't get what is my mistake and how to make it work.
I appreciate any one can help me with that
close all
clear all
clc
%%
%PR parameters
R=8.314;
T=352.49;
% Critical constants for normal butane (c_6)
Tc = 507.44; % Critical temperature in Kelvin
Pc = 30.31; % Critical pressure in bar
w = 0.299; % Acentric factor
% MW=0.05812;
m=0.37464+1.54226*w-0.26992*w^2;
a=0.45724*R^2*Tc^2/Pc*(1+m*(1-sqrt(T/Tc)))^2;
b=0.07780*R*Tc/Pc;
%Influenceparameter
m1=10^(-16)/(1.2326+1.3757*w);
m2=10^(-16)/(0.9051+1.5410*w);
c=a*b^(2/3)*(m1*(1-T/Tc)+m2);
%% Parameters
L = 0.000000008; % Length of the domain
T = 0.1; % Total time
Nx = 5; % Number of spatial grid points
Nt = 10; % Number of time steps
% c = 1; % Thermal diffusivity constant
dx = L / (Nx - 1); % Spatial step size
dt = T / Nt; % Time step size
% Initial and boundary conditions
u0 = @(x) 7000-(7000/L)*x; % Initial condition
uL = 7000; % Dirichlet boundary condition at x=0
uR = 0; % Dirichlet boundary condition at x=L
% Nonlinear source term and its convex/concave splitting
f = @(u) R.*T*(log(u./(1-u.*b))+u.*b./(1-u.*b))-u.*a./(1+2.*u.*b-u.^2*b.^2)-a./(2.*sqrt(2).*b).*log((1+(1-sqrt(2)).*u.*b)./(1+(1+sqrt(2)).*u.*b)); % Example: cubic nonlinearity
fc = @(u) R.*T.*(log(u./(1-u.*b))+u.*b./(1-u.*b)); % Convex part
fd = @(u) f(u) - fc(u); % Concave part
fplot(f)
% Derivative of the convex part for Newton's method
dfc = @(u) R.*T*(1/(u.*(b.*u-1)).^2); % Derivative of the convex part
% Exact solution function (if available)
% For demonstration, assuming an exact solution function
% exact_solution = @(x, t) 1./(1+exp(x-t));
% Spatial grid
x = linspace(0, L, Nx);
% u=zeros(Nt,Nx);
% Initial condition
u = u0(x);
figure(1)
plot(x,u0(x))
grid on
title('Initial condition')
xlabel('x')
ylabel('u0')
%%
% Laplacian matrix (central difference)
e = ones(Nx,1);
A = spdiags([e -2*e e], -1:1, Nx, Nx) / dx^2;
% Apply boundary conditions to A
A(1,:) = 7000; A(1,1) = 7000; % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;
%%
% Time integration loop
for n = 1:Nt
% Apply boundary conditions at each time step
u(1) = 0; % Boundary at x=0
u(Nx) = 0; % Boundary at x=L
% Initial guess for Newton's method (use the previous time step)
U = u;
% Newton iteration loop
for k = 1:100
% Calculate the residual F(U^k)
F_Uk = U - dt * c * A * U - dt * fc(U) - (u + dt * fd(u));
% Calculate the Jacobian J(U^k)
J_Uk = speye(Nx) - dt * c * A - dt * spdiags(dfc(U), 0, Nx, Nx);
% Solve for Newton update dU
dU = J_Uk \ (-F_Uk);
% Update U with the Newton step
U = U + dU;
% Check for convergence
if norm(dU, inf) < 1e-6
break;
end
end
% Update solution for the next time step
u = U;
% Optionally, plot the solution at selected times
if mod(n, Nt/10) == 0 % Plot every 10% of the total time steps
figure;
plot(x, u, '-o', 'DisplayName', 'Numerical Solution');
xlabel('x');
ylabel('u');
title(sprintf('Solution at time t = %.2f', n * dt));
grid on;
end
end
0 Kommentare
Antworten (1)
Torsten
am 9 Sep. 2024
Bearbeitet: Torsten
am 9 Sep. 2024
Why don't you use "pdepe" ? Ten lines of code maximum, and you are done.
To get rid of the error in the above code, replace
% Spatial grid
x = linspace(0, L, Nx);
by
% Spatial grid
x = linspace(0, L, Nx).';
But other problems follow.
Why do you set A as below ? Shouldn't this be u or U ?
A(1,:) = 7000; A(1,1) = 7000; % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;
Are you sure about u(1) = 0 ? I thought you want to keep it at 7000.
% Apply boundary conditions at each time step
u(1) = 0; % Boundary at x=0
u(Nx) = 0; % Boundary at x=L
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!