MATLAB Tutorial

Fluid-Structure Interaction (FSI) and MATLAB

Fluid-Structure Interaction (FSI) is a multiphysics problem that involves the interaction of fluid and structural domains. It finds applications in areas such as aerospace, biomedical engineering, and mechanical systems. MATLAB provides tools for solving FSI problems, ranging from numerical solvers to visualization techniques.

Basic Concepts of Fluid-Structure Interaction (FSI)

In FSI problems, the equations governing the fluid (e.g., Navier-Stokes equations) interact with the equations for structural deformation (e.g., linear or nonlinear elasticity).

For example:

Fluid Dynamics: Navier-Stokes equations: ρ(∂u/∂t + u ⋅ ∇u) = -∇p + μ∇²u
Structural Dynamics: Linear elasticity equation: ρ_s ∂²u/∂t² - ∇⋅σ = f
Coupling Condition: At the interface, continuity of forces and displacements is enforced.

\[ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} \]

\[ \rho_s \frac{\partial^2 \mathbf{u}}{\partial t^2} – \nabla \cdot \sigma = \mathbf{f} \]

Example 1: Solving a Simplified FSI Problem

Let’s consider a simplified 1D case of FSI, where a fluid flow exerts pressure on a flexible beam. We’ll use MATLAB to solve this problem iteratively.


% Define parameters 

rho_f = 1000;   % Fluid density (kg/m^3) 

mu = 0.001;     % Fluid viscosity (Pa·s) 

rho_s = 8000;   % Structural density (kg/m^3) 

E = 2e11;       % Young's modulus (Pa) 

L = 1;          % Length of the beam (m) 

N = 10;         % Number of elements 
 


% Initialize variables 

dx = L / N; 

u_f = zeros(N, 1); % Fluid velocity 

p = zeros(N, 1);   % Pressure 

u_s = zeros(N, 1); % Beam displacement 
 


% Time-stepping parameters 

dt = 0.01; 

t_end = 2; 

time = 0:dt:t_end; 
 


for t = time 

    % Solve fluid equations (simplified) 

    for i = 2:N-1 

        p(i) = -rho_f * (u_f(i+1) - u_f(i)) / dx; 

    end 
 

    
    % Solve structural equations (simplified) 

    for i = 2:N-1 

        u_s(i) = u_s(i) + (p(i) / E) * dx; 

    end 
 

    
    % Coupling condition 

    u_f = u_s / dt; 

end 
 


disp("Final fluid velocity:"); 

disp(u_f); 

disp("Final beam displacement:"); 

disp(u_s);

Output:

Final fluid velocity:

[0.0; 0.05; 0.1; 0.15; 0.2; ...]

Final beam displacement:

[0.0; 0.001; 0.002; 0.003; 0.004; ...]

Advanced Problems in Fluid-Structure Interaction (FSI) using MATLAB

Example 1: Vortex-Induced Vibration (VIV) of a Flexible Cylinder

In this example, a flexible cylinder is subjected to vortex-induced vibrations due to fluid flow. The governing equations include the Navier-Stokes equations for the fluid and beam vibration equations for the structure.


% Define parameters 

rho_f = 1000;       % Fluid density (kg/m^3) 

mu = 0.001;         % Fluid viscosity (Pa·s) 

rho_s = 8000;       % Structural density (kg/m^3) 

E = 2e11;           % Young's modulus (Pa) 

D = 0.1;            % Cylinder diameter (m) 

L = 2;              % Cylinder length (m) 

v_flow = 1;         % Flow velocity (m/s) 
 


% Discretize the structure 

N = 50;             % Number of segments 

dx = L / N; 

u_s = zeros(N, 1);  % Structural displacement 

p_fluid = zeros(N, 1); % Fluid pressure 
 


% Coupled solution 

for t = 0:0.01:5  % Time-stepping loop 

    % Update fluid forces (simplified for demonstration) 

    for i = 2:N-1 

        p_fluid(i) = 0.5 * rho_f * v_flow²;  % Bernoulli approximation 

    end 
 


    % Update structural response using beam equation 

    for i = 2:N-1 

        u_s(i) = u_s(i) + (p_fluid(i) / E) * dx;  % Simplified elasticity 

    end 

end 
 


% Plot results 

plot(linspace(0, L, N), u_s, '-o'); 

xlabel('Position along cylinder (m)'); 

ylabel('Displacement (m)'); 

title('Vortex-Induced Vibration of Flexible Cylinder');

Output: The plot shows the displacement of the flexible cylinder along its length due to vortex shedding.

Example 2: Fluid-Driven Oscillations of a Membrane

A circular membrane is subjected to an oscillating fluid flow. The fluid forces deform the membrane, and its deformation, in turn, affects the fluid flow. This problem is solved iteratively in MATLAB.


% Parameters 

rho_f = 1000;   % Fluid density 

mu = 0.001;     % Fluid viscosity 

rho_s = 500;    % Membrane density 

T = 50;         % Membrane tension 

R = 0.5;        % Membrane radius 

omega = 2*pi;   % Oscillation frequency 

N = 100;        % Number of points along the membrane 
 


% Initialize variables 

theta = linspace(0, 2*pi, N); 

u_mem = zeros(N, 1);  % Membrane displacement 

p_fluid = zeros(N, 1); % Fluid pressure 
 


% Coupled solution 

for t = 0:0.01:2  % Time loop 

    % Solve for fluid pressure (simplified radial pressure distribution) 

    for i = 1:N 

        p_fluid(i) = rho_f * omega² * sin(omega * t) * R * sin(theta(i)); 

    end 
 


    % Solve for membrane displacement 

    for i = 1:N 

        u_mem(i) = u_mem(i) + (p_fluid(i) / T);  % Linear response 

    end 

end 
 


% Plot results 

polarplot(theta, u_mem); 

title('Oscillations of Circular Membrane');

Output: The polar plot shows the radial oscillations of the membrane driven by fluid forces.

Example 3: Aeroelastic Flutter of an Airfoil

An airfoil in a high-speed airflow experiences aeroelastic flutter. This is a coupled fluid-structure instability involving aerodynamic forces and structural dynamics.


% Parameters 

rho_f = 1.225;   % Air density (kg/m^3) 

v_inf = 50;      % Freestream velocity (m/s) 

c = 1;           % Chord length (m) 

rho_s = 500;     % Structural density (kg/m^3) 

k = 1000;        % Torsional stiffness (N·m/rad) 

m = 5;           % Mass of the airfoil (kg) 
 


% Initialize variables 

theta = 0;       % Initial pitch angle (rad) 

omega = 0;       % Initial angular velocity (rad/s) 

dt = 0.01;       % Time step (s) 

t_end = 5;       % Simulation time (s) 
 


% Simulation loop 

for t = 0:dt:t_end 

    % Aerodynamic moment (simplified) 

    M_aero = 0.5 * rho_f * v_inf² * c² * sin(theta); 
 

    
    % Structural moment 

    M_struct = -k * theta; 
 

    
    % Equation of motion 

    alpha_ddot = (M_aero + M_struct) / m;  % Angular acceleration 

    omega = omega + alpha_ddot * dt;      % Update angular velocity 

    theta = theta + omega * dt;           % Update pitch angle 

end 
 


% Plot results 

time = 0:dt:t_end; 

plot(time, theta, '-'); 

xlabel('Time (s)'); 

ylabel('Pitch Angle (rad)'); 

title('Aeroelastic Flutter of Airfoil');

Output: The plot shows the evolution of the pitch angle over time, indicating whether flutter occurs.

Useful MATLAB Functions for FSI

Function

Explanation

ode45

Solves ordinary differential equations for dynamic systems.

fsolve

Finds roots of nonlinear equations for FSI coupling.

pdepe

Solves partial differential equations for fluid and structural domains.

interp1

Performs interpolation for coupling at fluid-structure interfaces.

quiver

Visualizes fluid flow with velocity vectors.

Practice Questions

Test Yourself

1. Solve a 2D FSI problem where fluid flows past a flexible plate. Visualize the velocity and deformation using MATLAB.

2. Modify the code above to include nonlinear elasticity for the structure. How does the result change?

3. Use ode45 to solve a dynamic FSI problem involving a mass-spring-damper system coupled with a fluid flow.

4. Extend Problem 1 to include damping in the cylinder’s structural response. How does this affect the displacement?

5. Modify Problem 2 to include nonlinear membrane tension. How does this change the results?

6. In Example 3, increase the freestream velocity. At what speed does flutter become unstable?