How to solve the system of time dependent coupled PDE's?
5 views (last 30 days)
Show older comments
The system this paper (DOI: 10.1017/S0022112003003835) Thank you.
Fig10a(0.2,2050)
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
Answers (1)
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
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 :-)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!