Generate a pore network model with grain diameters

37 次查看（过去 30 天）
Mahima Horta
Mahima Horta 2024-11-16，7:32
回答： Arnav 2024-11-19，12:51
Ran in:
Hello everyone,
I would like to generate a pore network model with grain diameters having log normal distribution (Grain diameter range- 20 um to 500 um; Mean= 200 um, Std dev, 0.1*Mean). The centre to centre distance for the grains need to be maintained as 550 um (to ensure there is no overlap). The domain size is 10 mm X 10 mm. I am using the following code:
clc;
clear;
% Parameters
domain_size = 10000; % Domain size in micrometers (10 mm = 10000 µm)
c_c_distance = 550; % Center-to-center distance in micrometers (µm)
N = floor(domain_size / c_c_distance); % Grid size (18 x 18 = 324 pores)
num_pores = N * N; % Total number of pores
% Define range for grain diameters
min_diameter = 20; % Minimum diameter in micrometers (µm)
max_diameter = 500; % Maximum diameter in micrometers (µm)
% Log-normal distribution parameters
mu = log(sqrt(min_diameter * max_diameter)); % Mean in log-space
sigma = sqrt(log(max_diameter / min_diameter)) / 3.29; % Standard deviation in log-space
% Generate log-normal distributed diameters
diameters = lognrnd(mu, sigma, [num_pores, 1]);
% Enforce the diameter range [20 µm, 500 µm]
diameters(diameters < min_diameter) = min_diameter;
diameters(diameters > max_diameter) = max_diameter;
% Limit the diameters to avoid overlap (should be less than c_c_distance)
max_valid_diameter = c_c_distance - 10;
diameters(diameters >= max_valid_diameter) = max_valid_diameter;
% Generate grid coordinates
x = (0:N-1) * c_c_distance; % X-coordinates
y = (0:N-1) * c_c_distance; % Y-coordinates
[X, Y] = meshgrid(x, y);
% Flatten the grid and assign diameters
X = X(:);
Y = Y(:);
% Plotting the heterogeneous pore network
figure;
hold on;
axis equal;
xlim([0, domain_size]);
ylim([0, domain_size]);
title('Heterogeneous Pore Network Model (10 mm × 10 mm) with Log-normal Grain Diameters');
xlabel('X Position (µm)');
ylabel('Y Position (µm)');
% Plot pores using circles
for i = 1:num_pores
% Plot each pore as a circle
rectangle('Position', [X(i) - diameters(i)/2, Y(i) - diameters(i)/2, diameters(i), diameters(i)], ...
'Curvature', [1, 1], 'EdgeColor', 'b', 'FaceColor', [0.7, 0.8, 1], 'LineWidth', 0.5);
end
grid on;
hold off;
But seems like there are many small grain diameters than the larger ones. Also, it is seen that the grains are overlapping the axis when I change the size distribution. Can someone help me modify it or give ideas for a better one.
Any help is appreciated
Thanks in advance

请先登录，再回答此问题。

回答（1 个）

Arnav
Arnav 2024-11-19，12:51
Hi @Mahima Horta,
It seems like the formula that you are using for finding μ and σ is inaccurate. The current formula seems to be finding μ as the logarithm of the geometric mean of the minimum and maximum diameters, but you want the mean to be 200 µm and the standard deviation to be 0.1 times the mean. Note that the parameters μ and σ are not the mean and standard deviation in logarithmic scale. They are simply parameters of the distribution that need to be calculated from the actual mean and standard deviation in linear scale.
The correct formula for μ and σ in terms of (mean) and (standard deviation) for a Log-normal distribution is shown below:
and
You may refer to the following Wikipedia link for the same: https://en.wikipedia.org/wiki/Log-normal_distribution
Therefore, the code for calculating mu and sigma can be changed to the following:
mean_diameter = 200;
std_dev = 0.1 * mean_diameter;
mu = log(mean_diameter^2 / sqrt(std_dev^2 + mean_diameter^2));
sigma = sqrt(log(1 + (std_dev^2 / mean_diameter^2)));
Regarding the issue of overlapping pores with the axis, you can fix this by adjusting the centres by a distance equal to their radius. Here's how you can do that:
X(i) = min(max(X(i) - diameters(i)/2, 0), domain_size - diameters(i));
Y(i) = min(max(Y(i) - diameters(i)/2, 0), domain_size - diameters(i));
The resultant plot is shown below:
Hope it helps!

类别

在 Help Center 和 File Exchange 中查找有关 Annotations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by