How to solve the system of time dependent coupled PDE's?

5 views (last 30 days)
The system this paper (DOI: 10.1017/S0022112003003835) Thank you.
Fig10a(0.2,2050)
Iteration: 81000 | G(1/4, t): NaN Iteration: 82000 | G(1/4, t): NaN
function Fig10a(delta, Reynolds)
% Main function to solve for F and G and plot G(1/4, t)
% Initialize variables
R = Reynolds;
dt = 1/10;
nmax = 82001;
dy = 1/100;
yTarget = 1/4;
gValues = [];
Delta = delta;
H = @(t) 1 + Delta * cos(2 * t);
dH = @(t) -2 * Delta * sin(2 * t);
% Create the grid and differentiation matrices
ygrid = 0:dy:1;
ny = length(ygrid);
% Finite difference differentiation matrices (for dy1 and dy2)
dy2 = fdcoeffFDM2(ny, dy); % Second-order finite difference matrix
dy1 = fdcoeffFDM1(ny, dy); % First-order finite difference matrix
% Initialize variables for F and G
fvar0 = zeros(ny, nmax);
gvar0 = zeros(ny, nmax);
% Time-stepping loop
for i = 1:nmax-1
% Update variables for F and G at the current time step
fvar = fvar0(:, i);
gvar = gvar0(:, i);
% Define the equations
eqf = dy2 * fvar + (dy1 * fvar)./ygrid' - fvar./(ygrid'.^2) + ...
gvar * H((i + 1) * dt)^2;
eqg = (gvar - gvar0(:, i)) / dt - ...
dH((i + 1) * dt) / H((i + 1) * dt) * (dy1 * gvar)./ygrid' - ...
(dy1 * gvar) .* fvar / H((i + 1) * dt) + ...
1/H((i + 1) * dt) * (dy1 * fvar) .* gvar + ...
2./(ygrid' * H((i + 1) * dt)) .* fvar .* gvar - ...
(1/(R * H((i + 1) * dt)^2)) * (dy2 * gvar + dy1 * gvar./ygrid' - gvar./(ygrid'.^2));
% Apply boundary conditions
eqf(1) = fvar(1); % F(0, t) = 0 at y = 0
eqf(end) = fvar(end) + dH((i + 1) * dt); % F(1, t) = -H'(t) at y = 1
eqg(1) = gvar(1); % G(0, t) = 0 at t = 0
eqg(end) = dy1(end, :) * fvar - dH((i + 1) * dt); % F'(1, t) = H'(t) at y = 1
% Combine eqf and eqg into a single system
eqns = [eqf; eqg];
% Solve the system using backslash operator
sol = eqns; % Since we are already evaluating eqns as the result
% Check that the solution vector matches the expected size
if length(sol) ~= 2 * ny
error('The number of equations does not match the number of unknowns');
end
% Update variables for the next step
fvar0(:, i+1) = sol(1:ny);
gvar0(:, i+1) = sol(ny+1:end);
% Interpolate G(y, t) and store G(1/4, t)
if i >= 81000 && i <= 82001
gInterp = interp1(ygrid, gvar0(:, i+1), yTarget);
gValues = [gValues; i, gInterp];
end
% Display debugging information every 1000 iterations
if mod(i, 1000) == 0 & i>=81000
disp(['Iteration: ', num2str(i), ' | G(1/4, t): ', num2str(gInterp)]);
end
end
% Plot the values of G(1/4, t)
figure;
if isempty(gValues)
disp('No values of G(1/4, t) were recorded.');
else
plot(gValues(:,1), gValues(:,2), 'r-', 'LineWidth', 2);
grid on;
xlabel('t');
ylabel('G(1/4, t)');
title('G(1/4, t) from t = 8100 to t = 8200');
end
end
% Auxiliary function for second-order finite difference matrix
function D2 = fdcoeffFDM2(ny, dy)
e = ones(ny, 1);
D2 = spdiags([e -2*e e], -1:1, ny, ny) / (dy^2);
D2(1, :) = 0; D2(end, :) = 0; % Apply boundary conditions
end
% Auxiliary function for first-order finite difference matrix
function D1 = fdcoeffFDM1(ny, dy)
e = ones(ny, 1);
D1 = spdiags([-e e], [-1 1], ny, ny) / (2*dy);
D1(1, :) = 0; D1(end, :) = 0; % Apply boundary conditions
end
  4 Comments
sajjad
sajjad on 22 Sep 2024
Moved: Torsten on 22 Sep 2024
First i have to draw fig 10a and then Fig 8. Then I will use the outputs for my further research work.
sajjad
sajjad on 22 Sep 2024
In equation 3.9 he is describing the boundary conditions for G as well.

Sign in to comment.

Answers (1)

Torsten
Torsten on 22 Sep 2024
Edited: Torsten on 23 Sep 2024
ode23t seems to work for high Reynolds numbers. In case ode23t fails (e.g. for the low Reynolds number regime), I recommend the solver "radau" for MATLAB which seems to works for the complete Reynolds number regime, but is much slower than ode23t.
It can be downloaded from:
nx = 101;
xstart = 0.0;
xend = 1.0;
x = linspace(xstart,xend,nx).';
dx = (xend-xstart)/(nx-1);
tstart = 0;
tend = 8200;
tspan = [tstart,tend];
G0 = zeros(nx,1);
F0 = zeros(nx,1);
y0 = [G0;F0];
delta = 0.2;
Reynolds = 2050;
H = @(t) 1 + delta*cos(2*t);
Hdot = @(t) -delta*2*sin(2*t);
M = [eye(nx),zeros(nx);zeros(nx,2*nx)];
M(1,1) = 0;
M(nx,nx) = 0;
options = odeset('Mass',M);
[T,Y] = ode23t(@(t,y)fun(t,y,x,nx,dx,delta,Reynolds,H,Hdot),tspan,y0,options);
G = Y(:,1:nx);
F = Y(:,nx+1:2*nx);
idx = T >= 8100;
figure(1)
plot(T(idx),G(idx,26))
figure(2)
plot(T(idx),F(idx,26))
function dydt = fun(t,y,x,nx,dx,delta,Reynolds,H,Hdot)
%persistent iter
%if isempty(iter)
% iter = 0;
%end
%iter = iter + 1;
%if mod(iter,1000)==0
% t
% iter = 0;
%end
G = y(1:nx);
F = y(nx+1:2*nx);
% F
% Compute spatial derivatives
dFdx = zeros(nx,1);
d2Fdx2 = zeros(nx,1);
dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);
dFdx(nx) = Hdot(t);
d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;
%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;
Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);
d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;
% Compute temporal derivatives
dFdt = zeros(nx,1);
dFdt(1) = F(1);
dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...
F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);
dFdt(nx) = F(nx) + Hdot(t);
% G
% Compute spatial derivatives
dGdx = zeros(nx,1);
d2Gdx2 = zeros(nx,1);
dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);
d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;
% Compute temporal derivatives
dGdt = zeros(nx,1);
dGdt(1) = G(1);
dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...
1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...
1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...
2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...
(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...
(H(t)^2*Reynolds);
dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);
%Taken from the article, should be the same as set
%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));
dydt = [dGdt;dFdt];
end
  23 Comments
Torsten
Torsten on 8 Oct 2024
Can we modify our code or system for Non-Newtonian fluid case or some other type of flow?
Maybe - if you know the new equations :-)

Sign in to comment.

Products


Release

R2013a

Community Treasure Hunt

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

Start Hunting!