Plotting Elastic Modulus Surface for TPMS Structure

17 次查看（过去 30 天）
Kamran
Kamran 2024-3-21
编辑： VBBV 2024-3-22
Ran in:
Hello,
I have been trying to plot the elastic modulus surface of a TPMS-based solid structure. I have computed the stiffness matrix through numerical homogenization while applying the peroidic conditions but unable to think of a way to acquire the elastic modulus surface plots (sphere like shapes) for TPMS structure to characterize isotropic behaviour.
I have attempted plotting the obtained stiffness matrix for polar plots and elastic modulus the figures to which I am attaching below but I am not sure if this is the right way as my elastic modulus plot does not makes sense. Also,
Can anyone please help me with this!
clc
clear all
close all
% Define the stiffness matrix (evaluated using numerical homogenization method)
stiffness_matrix = [
4.03E+00 1.81E+00 1.81E+00 2.19E-03 7.34E-04 9.74E-03;
1.81E+00 4.02E+00 1.81E+00 7.35E-03 5.62E-03 2.56E-03;
1.81E+00 1.81E+00 4.03E+00 -1.97E-03 9.29E-03 4.43E-03;
2.19E-03 7.35E-03 -1.97E-03 1.11E+00 4.06E-03 2.89E-03;
7.34E-04 5.62E-03 9.29E-03 4.06E-03 1.11E+00 2.03E-03;
9.74E-03 2.56E-03 4.43E-03 2.89E-03 2.03E-03 1.11E+00
];
% Angles (in radians) for polar plot
angles = linspace(0, 2*pi, 100);
% Calculate stiffness values for each angle
stiffness = zeros(size(angles));
for i = 1:length(angles)
% Compute stiffness in the direction of the current angle
direction = [cos(angles(i)); sin(angles(i)); 0]; % Direction vector in 3D
stiffness(i) = direction' * stiffness_matrix(1:3,1:3) * direction;
end
% Plot polar plot for stiffness
figure;
polarplot(angles, stiffness);
title('Stiffness Polar Plot');
% Extract the 3x3 upper-left submatrix (representing the elastic modulus tensor)
elastic_modulus_matrix = stiffness_matrix(1:3,1:3);
% Define the range of angles (in radians) for elastic modulus surface
theta = linspace(0, pi, 200);
phi = linspace(0, 2*pi, 200);
[Theta, Phi] = meshgrid(theta, phi);
% Calculate the elastic modulus in different directions
E = zeros(size(Theta));
for i = 1:numel(Theta)
direction = [sin(Theta(i))*cos(Phi(i)); sin(Theta(i))*sin(Phi(i)); cos(Theta(i))];
E(i) = direction' * elastic_modulus_matrix * direction;
end
% Convert spherical coordinates to Cartesian coordinates
X = sin(Theta) .* cos(Phi);
Y = sin(Theta) .* sin(Phi);
Z = cos(Theta);
figure;
surf(X, Y, Z, E, 'EdgeColor', 'none');
colormap('jet');
title('Elastic Modulus Surface');
xlabel('X');
ylabel('Y');
zlabel('Z');
colorbar;
c = colorbar;
c.Label.String = 'Elastic Modulus (GPa)';
c.Ticks = linspace(min(E(:)), max(E(:)), 11);
set(gca, 'FontSize', 12);
set(c.Label, 'FontSize', 12);

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

回答（1 个）

VBBV
VBBV 2024-3-21 
direction = [cos(Theta(i))*sin(Phi(i)); sin(Theta(i))*sin(Phi(i)); cos(Phi(i))]
  3 个评论
显示 1更早的评论
Kamran
Kamran 2024-3-21
编辑：Kamran 2024-3-21
Hello, thank you for you response.
I have included this change into the code but the plot still returns a modulus surface that is not close to isotropic. For it to be close to isotropic the shade throughout the sephere should be very similar with a minor variation. The Zener ratio, which is a measure of isotropy, returns a value of 0.99 (1 being isotropic and 0 being anisotropic). Also, I am using just 3×3 (upper left part) of the stiffness matrix becasue I cannot comprehend a way to to match the direction vector array size of 3×1 to the stiffness matrix array size of 6×6.
What is possibly going wrong in my code?
VBBV
VBBV 2024-3-22
编辑：VBBV 2024-3-22
For isotropic materials, you need to multply the resultant elastic modulus with poissons ratio, as
nu = 0.25; % poisson ratio for isotripic materials
E = E*nu/((1+nu)*(1-2*nu));
which would result in smaller changes as you'd expect

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by