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.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. 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. #2d simulation, #navier-stokes, #compressible fluid, #venting system MATLAB Answers — New Questions

​

I want to select the first element of a vector without having to store it in a variable. For example, if i want to know the smallest prime divider of a number, i want to write this
factor(56)(1)
However this is invalid. Any correct way to do this?I want to select the first element of a vector without having to store it in a variable. For example, if i want to know the smallest prime divider of a number, i want to write this
factor(56)(1)
However this is invalid. Any correct way to do this? I want to select the first element of a vector without having to store it in a variable. For example, if i want to know the smallest prime divider of a number, i want to write this
factor(56)(1)
However this is invalid. Any correct way to do this? indexing, vector MATLAB Answers — New Questions

​

Hello,

I am trying to test my model using simulink test, by placing it within a resettable subsystem, however only the MIL mode passes and the SIL mode fails for the reset functionality with a message "Simulation Error: Subsystem Model_xxx contains referenced models, which is not supported for SIL and PIL simulations.". The model within the resettable subsystem contains 2 Unit Delay blocks with the ‘Initial Condition’ set which I want to reset. I would appreciate any hints or alternate options to test the rest functionality in SIL mode.

Thanks,
Jerome.Hello,

I am trying to test my model using simulink test, by placing it within a resettable subsystem, however only the MIL mode passes and the SIL mode fails for the reset functionality with a message "Simulation Error: Subsystem Model_xxx contains referenced models, which is not supported for SIL and PIL simulations.". The model within the resettable subsystem contains 2 Unit Delay blocks with the ‘Initial Condition’ set which I want to reset. I would appreciate any hints or alternate options to test the rest functionality in SIL mode.

Thanks,
Jerome. Hello,

I am trying to test my model using simulink test, by placing it within a resettable subsystem, however only the MIL mode passes and the SIL mode fails for the reset functionality with a message "Simulation Error: Subsystem Model_xxx contains referenced models, which is not supported for SIL and PIL simulations.". The model within the resettable subsystem contains 2 Unit Delay blocks with the ‘Initial Condition’ set which I want to reset. I would appreciate any hints or alternate options to test the rest functionality in SIL mode.

Thanks,
Jerome. simulink, sil, resettale subsystem, mil MATLAB Answers — New Questions

​

Hello , I have two CSV files attached in ZIP format.one of the colums is from -1 to 1.
Is the a way in matlab to import the CSV ,ater that shift -1 to 1 coulms into 0 to 2?
adding 1 to each cell of the column?
then saving the CSV in updated form?
Thanks.Hello , I have two CSV files attached in ZIP format.one of the colums is from -1 to 1.
Is the a way in matlab to import the CSV ,ater that shift -1 to 1 coulms into 0 to 2?
adding 1 to each cell of the column?
then saving the CSV in updated form?
Thanks. Hello , I have two CSV files attached in ZIP format.one of the colums is from -1 to 1.
Is the a way in matlab to import the CSV ,ater that shift -1 to 1 coulms into 0 to 2?
adding 1 to each cell of the column?
then saving the CSV in updated form?
Thanks. csv, shift MATLAB Answers — New Questions

​

Hi
I have a guassian point spread function and i would like to find the equivalent MTF in frequency domain.
I have the following code with assistance from : https://www.mathworks.com/matlabcentral/answers/274848-point-spread-function-translate-to-mtf

x=[1,2,3,4,5,6,7,8,9,10,11,12,13]
y=[0,1,3,21,46,156,222,140,111,25,14,3,2]
LSF=[x,y];
OTF = fftshift(fft(LSF)); % OTF
MTF = abs(OTF); % absolute value of OTF
MTF = MTF./max(MTF); % normalize
% correct sampling frequency for conversion on frequency bins to frequency
Size = length(MTF);
spacing = 1.12e-6;% pixel size
fsx = abs(1/spacing); % turn into sampling frequency
a = linspace(-Size/2,Size/2,Size); % form scale for conversion based on frequency bins
conversionx = fsx./Size; % conversion factor for frequency bin units to frequency (unit^-1)
Psi = a.*conversionx; % frequency (unit^-1)

figure
plot(Psi,MTF)
xlim([0,Psi(end)])
xlabel(‘cycles/unit’)
ylabel(‘MTF’)

could someone please help me convert the X axis from cycles/unit to lp/mm , i am not sure how this conversion is peformedHi
I have a guassian point spread function and i would like to find the equivalent MTF in frequency domain.
I have the following code with assistance from : https://www.mathworks.com/matlabcentral/answers/274848-point-spread-function-translate-to-mtf

x=[1,2,3,4,5,6,7,8,9,10,11,12,13]
y=[0,1,3,21,46,156,222,140,111,25,14,3,2]
LSF=[x,y];
OTF = fftshift(fft(LSF)); % OTF
MTF = abs(OTF); % absolute value of OTF
MTF = MTF./max(MTF); % normalize
% correct sampling frequency for conversion on frequency bins to frequency
Size = length(MTF);
spacing = 1.12e-6;% pixel size
fsx = abs(1/spacing); % turn into sampling frequency
a = linspace(-Size/2,Size/2,Size); % form scale for conversion based on frequency bins
conversionx = fsx./Size; % conversion factor for frequency bin units to frequency (unit^-1)
Psi = a.*conversionx; % frequency (unit^-1)

figure
plot(Psi,MTF)
xlim([0,Psi(end)])
xlabel(‘cycles/unit’)
ylabel(‘MTF’)

could someone please help me convert the X axis from cycles/unit to lp/mm , i am not sure how this conversion is peformed Hi
I have a guassian point spread function and i would like to find the equivalent MTF in frequency domain.
I have the following code with assistance from : https://www.mathworks.com/matlabcentral/answers/274848-point-spread-function-translate-to-mtf

x=[1,2,3,4,5,6,7,8,9,10,11,12,13]
y=[0,1,3,21,46,156,222,140,111,25,14,3,2]
LSF=[x,y];
OTF = fftshift(fft(LSF)); % OTF
MTF = abs(OTF); % absolute value of OTF
MTF = MTF./max(MTF); % normalize
% correct sampling frequency for conversion on frequency bins to frequency
Size = length(MTF);
spacing = 1.12e-6;% pixel size
fsx = abs(1/spacing); % turn into sampling frequency
a = linspace(-Size/2,Size/2,Size); % form scale for conversion based on frequency bins
conversionx = fsx./Size; % conversion factor for frequency bin units to frequency (unit^-1)
Psi = a.*conversionx; % frequency (unit^-1)

figure
plot(Psi,MTF)
xlim([0,Psi(end)])
xlabel(‘cycles/unit’)
ylabel(‘MTF’)

could someone please help me convert the X axis from cycles/unit to lp/mm , i am not sure how this conversion is peformed image analysis, image processing, resolution MATLAB Answers — New Questions

​

Hey!
I am working on a project currently, where I would like to evaluate the conndensed water mass at the outer surface of a pipe.
To give you an example: A cold fluid below 0 °C / 32 °F is flowing through a pipe. The pipe is placed into a container with moist air. Assuming that the wall temperature is equal to the temperature of the fluid inside the pipe, water included in moist air should eventually condensate at the outer wall surface.
Hence my question, does anyone know how to model such a behaviour? I know how to model condensation inside the pipe, which is not of interest now, but I am not aware how to model the condensation inside the container at the outer pipe surface. Is this even possible?
Thank you for helping me out in advance!Hey!
I am working on a project currently, where I would like to evaluate the conndensed water mass at the outer surface of a pipe.
To give you an example: A cold fluid below 0 °C / 32 °F is flowing through a pipe. The pipe is placed into a container with moist air. Assuming that the wall temperature is equal to the temperature of the fluid inside the pipe, water included in moist air should eventually condensate at the outer wall surface.
Hence my question, does anyone know how to model such a behaviour? I know how to model condensation inside the pipe, which is not of interest now, but I am not aware how to model the condensation inside the container at the outer pipe surface. Is this even possible?
Thank you for helping me out in advance! Hey!
I am working on a project currently, where I would like to evaluate the conndensed water mass at the outer surface of a pipe.
To give you an example: A cold fluid below 0 °C / 32 °F is flowing through a pipe. The pipe is placed into a container with moist air. Assuming that the wall temperature is equal to the temperature of the fluid inside the pipe, water included in moist air should eventually condensate at the outer wall surface.
Hence my question, does anyone know how to model such a behaviour? I know how to model condensation inside the pipe, which is not of interest now, but I am not aware how to model the condensation inside the container at the outer pipe surface. Is this even possible?
Thank you for helping me out in advance! pipe (ma), condensation MATLAB Answers — New Questions

​

Hey everybody,
I have some data, that I want ot parse in an mySQL database. The data consists of a timestamp which has millisecond precision and some other colums. The data is stored in a table.
I use mysql-connector-java-5.1.47mysql-connector-java-5.1.47-bin.jar as described to connect to the database and write with:
sqlwrite(conn,’Messdaten’,testtable(1:5,:))
Inserting into the database works, however the millisecond precision is gone altough I declared the column as datetime(3). The data that is being stored has the e.g. the value 2019-02-16 10:59:46.000 instead of 2019-02-16 10:59:46.425
Can somebody please gibe me a hint what I am doing wrong ?
Thanks in advance!Hey everybody,
I have some data, that I want ot parse in an mySQL database. The data consists of a timestamp which has millisecond precision and some other colums. The data is stored in a table.
I use mysql-connector-java-5.1.47mysql-connector-java-5.1.47-bin.jar as described to connect to the database and write with:
sqlwrite(conn,’Messdaten’,testtable(1:5,:))
Inserting into the database works, however the millisecond precision is gone altough I declared the column as datetime(3). The data that is being stored has the e.g. the value 2019-02-16 10:59:46.000 instead of 2019-02-16 10:59:46.425
Can somebody please gibe me a hint what I am doing wrong ?
Thanks in advance! Hey everybody,
I have some data, that I want ot parse in an mySQL database. The data consists of a timestamp which has millisecond precision and some other colums. The data is stored in a table.
I use mysql-connector-java-5.1.47mysql-connector-java-5.1.47-bin.jar as described to connect to the database and write with:
sqlwrite(conn,’Messdaten’,testtable(1:5,:))
Inserting into the database works, however the millisecond precision is gone altough I declared the column as datetime(3). The data that is being stored has the e.g. the value 2019-02-16 10:59:46.000 instead of 2019-02-16 10:59:46.425
Can somebody please gibe me a hint what I am doing wrong ?
Thanks in advance! sql, sqlwrite, millisecond precision, mysql, mysql-connector-java-5.1.47 MATLAB Answers — New Questions

​

When I am ready to use imufilter, "Tune imufilter to Optimize Orientation Estimate",the current latest version R2024b cannot reproduce the documented example results? Please fix it in time!
TEST in MATLAB online and desktop MATLAB
ld = load(‘imufilterTuneData.mat’);
qTrue = ld.groundTruth.Orientation; % true orientation

fuse = imufilter;
qEstUntuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

% reset(fuse);
cfg = tunerconfig(‘imufilter’);
tune(fuse, ld.sensorData, ld.groundTruth, cfg)

qEstTuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

dUntuned = rad2deg(dist(qEstUntuned, qTrue));
dTuned = rad2deg(dist(qEstTuned, qTrue));
rmsUntuned = sqrt(mean(dUntuned.^2))
rmsTuned = sqrt(mean(dTuned.^2))

N = numel(dUntuned);
t = (0:N-1)./ fuse.SampleRate;
plot(t, dUntuned, ‘r’, t, dTuned, ‘b’);
legend(‘Untuned’, ‘Tuned’);
title(‘imufilter – Tuned vs Untuned Error’)
xlabel(‘Time (s)’);
ylabel(‘Orientation Error (degrees)’);When I am ready to use imufilter, "Tune imufilter to Optimize Orientation Estimate",the current latest version R2024b cannot reproduce the documented example results? Please fix it in time!
TEST in MATLAB online and desktop MATLAB
ld = load(‘imufilterTuneData.mat’);
qTrue = ld.groundTruth.Orientation; % true orientation

fuse = imufilter;
qEstUntuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

% reset(fuse);
cfg = tunerconfig(‘imufilter’);
tune(fuse, ld.sensorData, ld.groundTruth, cfg)

qEstTuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

dUntuned = rad2deg(dist(qEstUntuned, qTrue));
dTuned = rad2deg(dist(qEstTuned, qTrue));
rmsUntuned = sqrt(mean(dUntuned.^2))
rmsTuned = sqrt(mean(dTuned.^2))

N = numel(dUntuned);
t = (0:N-1)./ fuse.SampleRate;
plot(t, dUntuned, ‘r’, t, dTuned, ‘b’);
legend(‘Untuned’, ‘Tuned’);
title(‘imufilter – Tuned vs Untuned Error’)
xlabel(‘Time (s)’);
ylabel(‘Orientation Error (degrees)’); When I am ready to use imufilter, "Tune imufilter to Optimize Orientation Estimate",the current latest version R2024b cannot reproduce the documented example results? Please fix it in time!
TEST in MATLAB online and desktop MATLAB
ld = load(‘imufilterTuneData.mat’);
qTrue = ld.groundTruth.Orientation; % true orientation

fuse = imufilter;
qEstUntuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

% reset(fuse);
cfg = tunerconfig(‘imufilter’);
tune(fuse, ld.sensorData, ld.groundTruth, cfg)

qEstTuned = fuse(ld.sensorData.Accelerometer, …
ld.sensorData.Gyroscope);

dUntuned = rad2deg(dist(qEstUntuned, qTrue));
dTuned = rad2deg(dist(qEstTuned, qTrue));
rmsUntuned = sqrt(mean(dUntuned.^2))
rmsTuned = sqrt(mean(dTuned.^2))

N = numel(dUntuned);
t = (0:N-1)./ fuse.SampleRate;
plot(t, dUntuned, ‘r’, t, dTuned, ‘b’);
legend(‘Untuned’, ‘Tuned’);
title(‘imufilter – Tuned vs Untuned Error’)
xlabel(‘Time (s)’);
ylabel(‘Orientation Error (degrees)’); sensorfusion, imu MATLAB Answers — New Questions

​

I am connected to a database using mym. I want to update row based on conditions the sql for that would be
update dbname.dbtable
set Name = ‘A’
where Age = 20
What is the easiest most elegant way to do that on MATLAB? I want to update table in a MATLAB script where some calculations are being done.I am connected to a database using mym. I want to update row based on conditions the sql for that would be
update dbname.dbtable
set Name = ‘A’
where Age = 20
What is the easiest most elegant way to do that on MATLAB? I want to update table in a MATLAB script where some calculations are being done. I am connected to a database using mym. I want to update row based on conditions the sql for that would be
update dbname.dbtable
set Name = ‘A’
where Age = 20
What is the easiest most elegant way to do that on MATLAB? I want to update table in a MATLAB script where some calculations are being done. mysql MATLAB Answers — New Questions

​

Post Content Post Content aek rev2 project files MATLAB Answers — New Questions

​