Plot directional Young's modulus and it's inverse in 3D #mtexScript

stiffness3D

%% Plot directional Young's modulus and it's inverse in 3D for HCP
%this script takes values of the stiffness tensor and creates a 3d plot of
%the values of S, or 1/S which equals the directional Young's modulus.  It
%is currently set up for HCP, and two sets of values are provided.

%note that the surface plot method used here has trouble with 'undercut' shapes.
%since it works fine for my purposes I do not intend to alter it further,
%and have plotted by splitting into two surfaces and using hold on.  This
%does not work well for the Zn figures which are dumbbell shaped.

%If you plan to use this with eg cubic materials, significant changes will
%be necessary, as HCP only depends on the angle between the applied force
%and the C axis, while cubic depends on the angle to all three axes.  Also
%the equation in the loop will be different (see references at end), and of
%course the values for the compliance tensor will all be different.


%values for Magnesium
S11=22;
S33=19.7;
S13=-4.96;
S44=60.98;
% %values for Zinc
% S11=8.07;
% S33=27.55;
% S13=-7.02;
% S44=25.25;

% create array, of 3d angles
%populate list of inclinations (angle to c axis of crystal)
theta=-90:5:90;
%convert to radians
theta=deg2rad(theta);
%populate list of azimuth angles, in xy plane-note 0=360
az=0:5:355;
az=deg2rad(az);
%mesh for full sphere
[az,theta]=meshgrid(az,theta);
%linear arrays
theta = reshape(theta,[],1);
az = reshape(az,[],1);
%convert to cartesian both to allow plotting, and to use as input to the S
%calculation which uses unit vectors
[x,y,z] = sph2cart(az,theta,1);

% %can check sphere coverage here if desired
% figure;scatter3(x,y,z)
% axis equal
% xlabel('x')
% ylabel('y')
% zlabel('z')

%Equation uses cosine of the angle between the c axis and the direction.  Theta=0 in x-y
%plane, want a value of 90 here instead.
ang=cos(deg2rad(90)-theta);

%create a blank variable
S=zeros(length(ang),1);

%calculate local value of S
for i=1:length(theta)
%equation from R295, see references
S(i)=(1-ang(i)^2)^2*S11+ang(i)^4*S33+(ang(i)^2)*(1-ang(i)^2)*(2*S13+S44);
end

%convert to directional young's modulus
 e=1./S;

%select which one you want to plot
%note; if issues with surface plot, try swapping x and z
[x,y,z] = sph2cart(az,theta,S);
figure;scatter3(x,y,z,'.')
title('S relative to c axis- Mg')
axis equal

%to plot a surface, we need to split it into two halves
test1=z>=0;
x1=x(test1);
y1=y(test1);
z1=z(test1);
test2=z<=0;
x2=x(test2);
y2=y(test2);
z2=z(test2);

figure;

%plot3(x1,y1,z1,'.-')
tri1 = delaunay(x1,y1);
%figure;plot(x1,y1,'.')
%[r,c] = size(tri1);
%disp(r)
h = trisurf(tri1, x1, y1, z1);
tri2 = delaunay(x2,y2);
hold on
h = trisurf(tri2, x2, y2, z2);
hold off
axis equal
title('S relative to c axis- Mg')
axis off


[x,y,z] = sph2cart(az,theta,e);
figure;scatter3(x,y,z)
title('1/S relative to c axis- Mg')
axis equal


%to plot a surface, we need to split it into two halves
test1=z>=0;
x1=x(test1);
y1=y(test1);
z1=z(test1);
test2=z<=0;
x2=x(test2);
y2=y(test2);
z2=z(test2);

figure;
%plot3(x1,y1,z1,'.-')
tri1 = delaunay(x1,y1);
%figure;plot(x1,y1,'.')
%[r,c] = size(tri1);
%disp(r)
h = trisurf(tri1, x1, y1, z1);
tri2 = delaunay(x2,y2);
hold on
h = trisurf(tri2, x2, y2, z2);
hold off
title('1/S relative to c axis- Mg')
axis equal
 axis off
 
% REFERENCES - INFO ON ON EQUATIONS NECESSARY
%R295- includes equations for cubic and HCP (10/1 and 10/4) on p21 and
%figures on p.193.  Best explanation of the equation inputs.
%Schmid & Boas Plasticity of Crystals Chapman & Hall, 1968

%R284- equations for many crystallographic systems on p.2420 and figures 
%for many materials on p.2421,2422
%note definition of l incorrect- using unit vector gives sin of angle, need cos instead 
%Buschow, K. (Ed.) Elastic Behaviour of Polycrystals The Encyclopedia of Materials: Science and Technology, Elsevier, 2001

%R246- avalable online, but only HCP.  Many material constants.  HCP
%equation (13) different notation but same result.
%Tromans, D. ELASTIC ANISOTROPY OF HCP METAL CRYSTALS AND POLYCRYSTALS International Journal of Recent Research and Applied Studies, 2011, 6, 462-483

Hi Jessica,

what is the difference to this one? Can your code be translated into something similar to this?

cs = loadCIF('Mg-Magnesium.cif')
S = zeros(6,6);
S(1,1)=22;
S(2,2)=22;
S(3,3)=19.7;
S(1,3)=-4.96;
S(1,3)=-4.96;
S(2,3)=-4.96;
S(3,2)=-4.96;
S(1,2)=-7.75;
S(4,4)=60.98;
S(5,5)=60.98;
S(6,6)=2*(S(1,1)-S(1,2));
S =stiffnessTensor(S,cs)

y = YoungsModulus(S)

surf(1./y)