Plotting Elastic Modulus Surface for TPMS Structure

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);

Answer 1

[VBBV]

direction = [cos(Theta(i))*sin(Phi(i)); sin(Theta(i))*sin(Phi(i)); cos(Phi(i))]

[VBBV]

The direction in spherical coordinates is something different in your case.

[Kamran]

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]

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