I wanna make 2D compressible fluid simulation when off gas ventilation at battery thermal runaway

clear; clc; close all;

% Simulation type:

% 1 - Particle Advection

% 2 - Velocity Heat Map

% 3 - Curl Heat Map

TYPE = 1;

% Simulation parameters

s = 100; % Grid size

ar = 5; % Aspect ratio

n = 120000; % Maximum number of particles

timestep = 0.01; % Time step

% Physical constants

initial_density = 1.225; % kg/m^3 (Air at sea level)

initial_temperature = 300; % K

initial_pressure = 101325; % Pa

R = 287; % Specific gas constant for air (J/kg/K)

viscosity = 1.8e-5; % Dynamic viscosity (Pa.s)

heat_diffusion = 0.001; % Heat diffusion coefficient

% Vent properties for thermal runaway

vent_temperature = 500; % High temperature vented gas (K)

vent_density = 0.5; % Low density vented gas (kg/m^3)

vent_velocity = 1; % High velocity out of the vent (m/s)

%vent_pressure_threshold = 2 * initial_pressure; % Vent opens when pressure exceeds this value

% Initialize colormap

jetMap = colormap(jet);

negJetMap = flipud(jetMap);

% Define vent/cube vertices

a = 40;

b = 60;

c = 50;

d = 55;

% Create grid

[X, Y] = meshgrid(1:s*ar, 1:s);

% Initialize fields for compressible flow

rho = ones(s, s*ar) * initial_density; % Density

T = ones(s, s*ar) * initial_temperature; % Temperature

p = ones(s, s*ar) * initial_pressure; % Pressure

vx = zeros(s, s*ar); % Velocity in x

vy = zeros(s, s*ar); % Velocity in y

% Initial positions of particles

[px, py] = meshgrid(0:d-c, a:b);

% Store initial positions for inflow

pxo = px;

pyo = py;

% Set up figure for visualization

f = figure(1);

set(f, 'WindowState', 'maximized');

% Simulation parameters

total_time = 10; % Total simulation time (adjust this as necessary)

num_steps = total_time / timestep; % Total number of time steps

% Set up figure for visualization

f = figure(1);

set(f, 'WindowState', 'maximized');

% Main simulation loop (for a fixed number of time steps)

for step = 1:num_steps

start = tic;

vent_pressure_threshold = sum(abs(p(:))); % Vent opens when pressure exceeds this value

if max(sum(p(1:d-c, a:b))) >= vent_pressure_threshold

% Update the venting region with high-temperature, high-velocity gas

T(1:d-c, a:b) = vent_temperature; % High temperature for vented gas

rho(1:d-c, a:b) = vent_density; % Low density for vented gas

vx(1:d-c, a:b) = vent_velocity; % High velocity out of the vent in x-direction

vy(1:d-c, a:b) = 0; % No vertical velocity for now

end

% Apply boundary conditions (no-slip walls and open boundaries)

vx(1, :) = 0; vy(1, :) = 0; % Top boundary (no-slip)

vx(end, :) = 0; vy(end, :) = 0; % Bottom boundary (no-slip)

vx(:, 1) = 0; vy(:, 1) = 0; % Left boundary (no-slip)

vx(:, end) = vx(:, end-1); % Open boundary (allow gas to exit freely)

vy(:, end) = vy(:, end-1); % Open boundary (allow gas to exit freely)

% Update pressure using the equation of state (Ideal Gas Law)

p = rho .* R .* T;

% Continuity equation (mass conservation)

drho_dt = -divergence(rho .* vx, rho .* vy);

% Momentum equations for compressible Navier-Stokes (simplified)

dvx_dt = -(vx .* gradient(vx) + vy .* gradient(vx)) ...

- (1 ./ rho) .* gradient(p) ...

+ viscosity * del2(vx); % Viscous term

dvy_dt = -(vx .* gradient(vy) + vy .* gradient(vy)) ...

- (1 ./ rho) .* gradient(p) ...

+ viscosity * del2(vy);

% Energy equation (heat transfer due to venting gas) (simplified)

dT_dt = -vx .* gradient(T) - vy .* gradient(T) ...

+ heat_diffusion * del2(T);

% Update fields for the next timestep(matrix..?)

rho = rho + drho_dt * timestep;

vx = vx + dvx_dt * timestep;

vy = vy + dvy_dt * timestep;

T = T + dT_dt * timestep;

% Advect velocity field using Runge-Kutta 4th order method =>입자의 속도=vx,

% vy/ 속도장=pvx, pvy

[pvx, pvy] = RK4(X, Y, vx, vy, -1);

vx = interp2(vx, pvx, pvy, 'cubic', 0);

vy = interp2(vy, pvx, pvy, 'cubic', 0);

% Visualization and display updates can remain the same

if TYPE == 1

% Advect particles using Runge-Kutta 4th order method

[px, py] = RK4(px, py, vx, vy, timestep);

% Ensure px, py remain within the grid

px = max(min(px, s*ar), 1); % Constrain px within grid range

py = max(min(py, s), 1); % Constrain py within grid range

% 벤트 영역에 있는 입자들을 벤트 외부로 이동시키기 위한 조건 추가

% 벤트 내에 있는 입자들의 위치를 벤트 외부로 업데이트

inside_vent = (px >= 1 & px <= d-c & py >= a & py <= b);

px(inside_vent) = px(inside_vent) + vent_velocity * timestep;

% 벤트 영역을 벗어난 입자들이 계속 밖으로 나갈 수 있도록 처리

[pxo, pyo] = meshgrid(0:d-c, a:b);

% Add inflow particles only if vent opens (pressure threshold exceeded)

if max(max(p(a:b, c:d))) >= vent_pressure_threshold

px = [px; pxo];

py = [py; pyo];

end

% Plot particles

scatter(px(:), py(:), 1, 'filled');

axis equal;

axis([0 s*ar 0 s]);

xticks([]);

yticks([]);

hold on;

% Highlight vent/cube region

rectangle('Position', [0, a, d-c, b-a], 'EdgeColor', 'r', 'LineWidth', 2);

stop = toc(start);

FPS = 1/stop;

% Update title with simulation metrics

title(sprintf("Pressure: %.5f # Particles: %d FPS: %.2f", ...

sum(abs(p(:))) / (s*ar), length(px), FPS));

hold off;

drawnow;

elseif TYPE == 2

% Visualize velocity magnitude

velocity_Mag = sqrt(vx.^2 + vy.^2);

velocity_Mag = imresize(velocity_Mag, 10, 'bicubic');

imagesc(velocity_Mag);

colormap(negJetMap);

xticks([]);

yticks([]);

axis equal;

stop = toc(start);

FPS = 1/stop;

title(sprintf("FPS: %.2f", FPS));

drawnow;

elseif TYPE == 3

% Visualize curl field (vorticity)

CURL = abs(curl(vx, vy));

CURL = imresize(CURL, 10, 'bicubic');

imagesc(CURL);

colormap(negJetMap);

xticks([]);

yticks([]);

axis equal;

drawnow;

end

end

% Function for Runge-Kutta 4th order method for advection

function [x_new, y_new] = RK4(px, py, vx, vy, h)

k1x = interp2(vx, px, py, 'linear', 0);

k1y = interp2(vy, px, py, 'linear', 0);

k2x = interp2(vx, px + h/2 * k1x, py + h/2 * k1y, 'linear', 0);

k2y = interp2(vy, px + h/2 * k1x, py + h/2 * k1y, 'linear', 0);

k3x = interp2(vx, px + h/2 * k2x, py + h/2 * k2y, 'linear', 0);

k3y = interp2(vy, px + h/2 * k2x, py + h/2 * k2y, 'linear', 0);

k4x = interp2(vx, px + h * k3x, py + h * k3y, 'linear', 0);

k4y = interp2(vy, px + h * k3x, py + h * k3y, 'linear', 0);

x_new = px + h/6 * (k1x + 2*k2x + 2*k3x + k4x);

y_new = py + h/6 * (k1y + 2*k2y + 2*k3y + k4y);

end

Hi everyone, I am struggling with making code which can simulate 2D compressible fluid which is off gas when battery thermal runaway situation. I made initial condition for particles inside the rectangle area before vent isn't destructed. After pressure is over the certain pressure threshold, vent become destructed and particles should move outside of rectangle area.

Right now, code upon is just showing particles inside vent(red rectangle), and calculating pressure value is too massive. I wanna modify this code to make particles move outside of vent after certain pressure value.

Method of updating model is based on continuity equation, Navier-Stokes momentum equation, Energy equation & Runge-Kutta 4th method. Main problem is with type 1. Mainly, what I wanna show is particle movement after the vent is destructed. Please help me, I'm waiting for u guys advise.

Answer 1

[Yifeng Tang]

Just a suggestion ...

Right now, you have Neumann type of boundary condition for your normal velocity at the right boundary, namely

vx(:, end) = vx(:, end-1); % Open boundary (allow gas to exit freely)

vy(:, end) = vy(:, end-1); % Open boundary (allow gas to exit freely)

if you instead prescribe a boundary velocity based on the local pressure, i.e. to change it to a Dirichlet type of boundary condition, you may be able to achieve what you want.

Have a if statement based on pressure value at the boundary grid point. In case it exceed a threshold value, apply vx & vy to be a value calculated from the local pressure. You'll need some reasonable empirical formula for this.