Fourth order PDE describing the fluid motion

9 views (last 30 days)
I am trying to solve fourth order partial differential equation describing a fluid problem, by using MATLAB. The PDE along with boundary conditions is provided here as:
Initial conditon is not provided but they have choose the initial condition as mentioned in the article Page. No: 61 in the paragraph after equation (2.6) and at the end of Page. No: 61. So we choose
And,
is the derivative of H with respect to t.
My working for the given problem is:
% Define parameters
H = @(t) ...; % Function for H as a function of time
H_dot = @(t) ...; % Derivative of H with respect to time
R = ...; % Value of R
% Define the spatial variable and time variable
eta = linspace(eta_min, eta_max, num_points); % Range of eta
t = linspace(t_min, t_max, num_time_steps); % Time steps
% Define V as a matrix where rows correspond to eta and columns to time
V = zeros(num_points, num_time_steps);
% Initial and boundary conditions for V (example)
V(:, 1) = ...; % Initial condition for V
V(1, :) = ...; % Boundary condition at eta = eta_min
V(end, :) = ...; % Boundary condition at eta = eta_max
% Loop through time steps using finite difference method
for n = 1:num_time_steps-1
% Compute derivatives using finite differences
V_eta = diff(V(:, n))./diff(eta);
V_eta_eta = diff(V_eta)./diff(eta(1:end-1));
V_eta_eta_eta = diff(V_eta_eta)./diff(eta(1:end-2));
V_eta_eta_eta_eta = diff(V_eta_eta_eta)./diff(eta(1:end-3));
% Compute the right-hand side of the equation
RHS = (1/H(t(n))) * (2 * H_dot(t(n)) * V_eta_eta + eta .* H_dot(t(n)) .* V_eta_eta_eta) ...
+ (1 / (R * H(t(n))^2)) * V_eta_eta_eta_eta ...
+ (1 / H(t(n))) * (V .* V_eta_eta_eta - V_eta .* V_eta_eta);
% Update V for the next time step using Euler's method (or other methods)
V(2:end-3, n+1) = V_eta_eta + dt * RHS; % Adjust indices for interior points
end
% Plot the results
mesh(t, eta, V);
xlabel('Time (t)');
ylabel('Eta');
zlabel('V');
title('Solution for V(\eta, t)');
The Manuscript with title "The Onset of Chaos in a class of Navier-Stokes Solutions" is also provided here:
The above stated PDE is on page number 61 of the given manuscript against the equation number (2.4) along with linear and nonlinear terms presented in equations (2.5) and (2.6) respectively.
The code is complete in all aspect but still not working.
Need the guideline.
If you need more information or think I have missed something, please let me know.

Answers (1)

Umar
Umar on 13 Oct 2024

Hi @Muhammad Zeeshan ,

After going through your comments, analysis of your code and documentation provided at the link below.

https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/onset-of-chaos-in-a-class-of-navierstokes-solutions/ED5266323D6E6362E275CD809899A2C5

Since you have to purchase journal and not able to provide full details, this is how I would tackle the fourth-order partial differential equation (PDE) you have described, you need to set up a MATLAB script that incorporates the specified linear and nonlinear operators, boundary conditions, and initial conditions derived from the manuscript you provided. Your problem revolves around fluid dynamics, specifically a flow between parallel walls influenced by oscillatory motion. The behavior of the flow is affected by parameters such as amplitude (Δ) and Reynolds number (R), with chaotic behavior emerging under certain conditions. The PDE is significant as it models these dynamics and serves as a foundation for understanding complex flow phenomena. Below is an updated MATLAB code that defines the necessary parameters and implements a numerical method to solve the given PDE. The code uses finite difference methods or similar approaches for spatial discretization, and can be adjusted based on your specific needs for stability and accuracy.

% Parameters
L = 10;          % Length of the domain
T = 20;          % Total time duration
Nx = 100;        % Number of spatial points
Nt = 2000;       % Number of time steps
dx = L/(Nx-1);   % Spatial step size
dt = T/Nt;       % Time step size
% Physical parameters
delta = 0.1;     % Amplitude of oscillation
Re = 100;        % Reynolds number
% Define spatial grid
x = linspace(0, L, Nx);
% Initial condition: H(t)
H = @(t) 1 + delta * cos(2*t);
dHdt = @(t) -2 * delta * sin(2*t);  % Derivative of H with respect to t
% Initialize V field
V = zeros(Nx, Nt);
V(:, 1) = H(0);  % Set initial condition at t=0
% Time-stepping loop
for n = 1:Nt-1
  t = n * dt;
    % Compute V at next time step using explicit method (pseudo-code)
    for i = 2:Nx-1
        % Linear operator LV (eq. 2.5)
        LV = (V(i+1,n) - 2*V(i,n) + V(i-1,n)) / dx^2;
        % Nonlinear operator NV (eq. 2.6)
        NV = (1/H(t)) * (V(i,n) * V(i,n) * V(i,n) - V(i,n)^2);
        % Update V based on the PDE
        V(i,n+1) = V(i,n) + dt * (LV + NV);
    end
    % Apply boundary conditions
    V(1,n+1) = H(t);   % V(0,t) = H(t)
    V(end,n+1) = 0;    % V(L,t) = 0
  end
% Visualization of results
figure;
surf(linspace(0,T,Nt), x, V);
xlabel('Time');
ylabel('Position');
zlabel('Velocity Field');
title('Fluid Flow Between Parallel Walls');
colorbar;

Please see attached.

In the above script, parameters are initialized such as the length of the domain, total time duration, number of spatial points, and time. The initial condition for H(t) is set as specified in your context. A finite difference approach is utilized to approximate derivatives in space and time: The linear operator LV is calculated using central differences. The nonlinear operator NV is computed directly from the defined equation. Boundary conditions are applied after updating V at each time step. Finally, the results are visualized in a surface plot showing how velocity changes over time and space.

Make sure that your chosen time step dt and spatial resolution dx satisfy stability criteria for numerical methods employed. Depending on the specifics of your problem and computational resources, you may consider more sophisticated methods such as implicit schemes or spectral methods for better stability in solving higher-order PDEs. Conduct sensitivity analyses on delta and Re to explore different regimes of flow behavior, including stability and chaos.

This code at the moment provides a foundational framework to address the fluid dynamics problem you are investigating while adhering to the details outlined in the manuscript. Adjustments may be needed based on specific requirements or further insights from your research context.

  11 Comments
Muhammad  Zeeshan Khan
Muhammad Zeeshan Khan on 19 Oct 2024
Hi @Torsten @Umar Thank you for such a nice comment. I am facing still some issues regarding this code. The code has still some syntax errors. is it due to the MATLAB version difference?
See the attachment:
Here is the complete code as provided by you:
clc; clear; % Parameters
nu = 0.01; % Kinematic viscosity (cm^2/s)
n = 1; % Frequency of wall oscillation
b = 1; % Half-width of the channel
R_values = [0, 5, 25]; % Reynolds numbers to evaluate
Delta_values = [0.25, 0.45, 0.65]; % Oscillation amplitudes
t_end = 100; % Total time for simulation
dt = 0.01; % Time step
N = 500; % Number of spatial points
eta = linspace(-1, 1, N); % Spatial domain
% Initialize variables for results storage
V_results = cell(length(R_values), length(Delta_values));
% Time-stepping loop
for R_idx = 1:length(R_values)
R = R_values(R_idx);
for Delta_idx = 1:length(Delta_values)
Delta = Delta_values(Delta_idx);
% Initialize stream function V with boundary conditions
V = zeros(N, round(t_end/dt));
for m = 1:round(t_end/dt)-1
t = m * dt;
H_t = 1 + Delta * cos(2 * n * t); % Wall position
L_matrix = construct_L_matrix(nu, H_t, eta, R, N);
V(:, m+1) = L_matrix \ (V(:, m) + compute_N(V(:, m), eta, R));
if mod(m, 100) == 0
fprintf('Computing R=%d, ∆=%.2f at t=%.2f\n', R, Delta, t);
end
end
% Store results for plotting later
V_results{R_idx, Delta_idx} = V;
end
end
% Visualization
for R_idx = 1:length(R_values)
for Delta_idx = 1:length(Delta_values)
V_plot = V_results{R_idx, Delta_idx};
[X, T] = meshgrid(eta, dt*(0:round(t_end/dt)-1));
subplot(length(R_values), length(Delta_values), (R_idx-1)*length(Delta_values)+Delta_idx);
surf(X, T, V_plot', 'EdgeColor', 'none');
title(sprintf('R=%d, ∆=%.2f', R_values(R_idx), Delta_values(Delta_idx)));
xlabel('η'); ylabel('Time'); zlabel('Stream Function V');
view(30, 30); % Adjust view angle for better visualization
end
end
function L_matrix = construct_L_matrix(nu, H_t, eta, R, N)
% Construct the linear operator matrix based on given parameters
h = eta(2) - eta(1); % Spatial step size
L_matrix = zeros(N);
for i = 2:N-1
L_matrix(i,i-1) = nu/h^2; % Left neighbor term
L_matrix(i,i) = -2*nu/h^2 - H_t;
% Diagonal term including wall motion
L_matrix(i,i+1) = nu/h^2; % Right neighbor term
end
% Boundary condition (No slip)
L_matrix(1,:) = 0; L_matrix(1,1) = 1; % Boundary condition at η=-1
L_matrix(N,:) = 0; L_matrix(N,N) = 1; % Boundary condition at η=1
end
function N_term = compute_N(V_m, eta, R)
% Compute the nonlinear terms in the Navier-Stokes equations
N_term = zeros(size(V_m));
dV_deta2 = gradient(gradient(V_m));
dV_deta3 = gradient(dV_deta2);
N_term(2:end-1) = R * (dV_deta2(2:end-1) .* V_m(2:end-1) - V_m(2:end-1) .* dV_deta3(2:end-1));
return;
end
Torsten
Torsten on 19 Oct 2024
Edited: Torsten on 19 Oct 2024
The code is just an example to show you the structure on how you could approach your problem.
If you still didn't recognize this, I suggest you first take the MATLAB Onramp course to get used to the basics of the new software:
The problem is too hard for that you could hope that we do the work for you.

Sign in to comment.

Products


Release

R2018a

Community Treasure Hunt

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

Start Hunting!