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.
