solution of system of coupled partial differential equations

Answer 1

[Precise Simulation]

One way to approach this is, due to the 2nd derivative in time, to split the time dependent variables into

du1/dt = u1h
du1h/du = rest of the equation

as typical for the wave equation fem simulations, and using Hermite finite element shape functions to account for the 4th order derivative, as done here for the Euler-Bernoulli beam deformation simulation.

Answer 2

[mohamed]

clear; clc;

%% Parameters

L = 1; % Beam length (m)

EI = 2.1*10^5; % Flexural rigidity (Pa.m^4)

% rhoA = 2.5*2450*10^(-4); % Mass per unit length (kg/m)

C = .5; % Stability constant

p=100;

omega=2*pi;

Nx = 10; % Number of spatial points

x = linspace(0, L, Nx); % Spatial discretization

dx = x(2) - x(1); % Spatial step size

k = 0; % Stiffness coefficient

c = 0; % Damping coefficient

m = 78; % rhoA Mass coefficient

% Calculate time step size for stability

dt = C * (dx^2) / sqrt(EI / m); % Time step for stability

T = 1.0; % Total time

Nt = floor(T/dt); % Number of time steps

t = linspace(0,T, Nt ); %t values

% Beam initial conditions

w= zeros(Nt, Nx); % Displacement over time

v = zeros(1,Nx); % Initial velocity

%% Initial deflection (example: sine wave)

w(1, :) = sin(pi*x/L); % Initial condition for beam deflection

w( 2,:) = w(1,:);

% Load function

f = zeros(Nt , Nx );

for n = 1:Nt

f(n, :) =p* cos(omega * t(n));

end

%% Finite Difference Method

for n = 2:Nt-1

for i = 2:Nx-1 % Spatial points excluding boundaries

% Boundary conditions

if(i == 1 || i == Nx) % Outside Boundary

w(n,i) = 0;

elseif(i==2) % Near left boundary

w(n+1,i) = ((-EI / dx^4) * (w(n,i+2) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) - w(n,i)) - ...

2*k*w(n,i) + ...

(c/dt)*w(n-1,i) + ...

(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...

p * cos(omega * t(n))) / ...

((m / dt^2) + (2*c / (2*dt)));

elseif i == Nx-1 % Near right boundary

w(n+1,i) = ((-EI / dx^4) * (-w(n,i) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) + w(n,i-2)) - ...

2*k*w(n,i) + ...

(c/dt)*w(n-1,i) + ...

(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...

p * cos(omega * t(n))) / ...

((m / dt^2) + (2*c / (2*dt)));

else % Interior points

w(n+1,i) = ((-EI / dx^4) * (w(n,i+2) - 4*w(n,i+1) + 6*w(n,i) - 4*w(n,i-1) + w(n,i-2)) - ...

2*k*w(n,i) + ...

(c/dt)*w(n-1,i) + ...

(m / dt^2) * (2*w(n,i) - w(n-1,i)) + ...

p * cos(omega * t(n))) / ...

((m / dt^2) + (2*c / (2*dt)));

end

end

end

%% Plot results

figure;

for n = 1:100:Nt

plot(x, w(n,:));

title(['Time = ', num2str(n*dt), ' sec']);

xlabel('x (m)'); ylabel('Deflection (m)');

grid on;

pause(0.1);

end