initial work on bulkmodulus

TensorAnalysis/@complianceTensor/HashinShtrikmanModulus.m

@@ -0,0 +1,85 @@
+function [khs, ghs, minmax] = HashinShtrikmanModulus(C,ko,go)
+% Hashin Shtrikman moduli 
+%
+% The equations are from Peselnick and Meister 1965 through Watt and
+% Peselnick 1980 which are based on Hashin and Shtrikman 1963 JMB 8/2013
+%
+% Syntax:
+%
+%   [hs,value] = hscalc(ko,go,cij)
+%
+% Input
+%
+%  C  - @stiffnessTensor
+%  ko - bulk modulus of the reference isotropic material
+%  go - shear modulus of the reference isotropic material
+
+%
+% Output
+%  khs - Hashin Shtrikman bulk modulus
+%  ghs - Hashin Shtrikman shear modulus
+%  minmax - is +1 for positive definite and -1 for negative definite R
+%
+
+% the isotropic compliances elements  (eq. 10-11 Watt Peselnick 1980)
+alpha = -3 / (3*ko + 4*go);
+beta = -3 * (ko+2*go) / ( 5*go * (3*ko+4*go) );
+gamma = (alpha - 3 * beta) / 9;
+
+% set up isotropic elastic matrix 
+co = 2 * go * stiffnessTensor.eye;
+
+co{1:3,1:3} = co{1:3,1:3} + ko - 2/3 * go;
+
+% take difference between anisotropic and isotropic matrices and take
+% inverse to get the residual compliance matrix (eq 4 & 8 Watt Peselnick
+% 1980)
+R = C - co;
+H = inv(R);
+
+% Subtract the isotropic compliances from the residual compliance matrix
+% and invert back to moduli space.  The matrix B is the moduli matrix that
+% when multiplied by the strain difference between the reference isotropic
+% response and the isotropic response of the actual material gives the
+% average stres <pij> (eq 15 - 16 Watt Peselnick 1980) note identity matrix
+% is Iinv
+A = H - beta * complianceTensor.eye;
+A{1:3,1:3} = A{1:3,1:3} - gamma;
+B=matrix(inv(A),'Voigt');
+
+% now perform the averaging of the B matrix to get two numbers (B1 and B2)
+% that are related to the two isotropic moduli - perturbations of the
+% reference values (eq 21 and 22 Watt Peselnick 1980)
+sB1 = sum(B(1:3,1:3,:),[1,2]);
+
+sB2 = B(1,1,:) + B(2,2,:) + B(3,3,:) + 2 * (B(4,4,:) + B(5,5,:) + B(6,6,:));
+B1 = (2*sB1 - sB2) / 15;
+B2 = (3*sB2 - sB1) / 30;
+
+% The Hashin-Shtrikman moduli are then calculated as changes from the
+% reference isotropic body. (eq 25 & 27 Watt and Peselnick 1980)
+khs = ko + ( 3 * B1 + 2 * B2 ) / ( 3 + alpha * (3 * B1 + 2 * B2) );
+ghs = go + B2 / ( 1 + 2 * beta * B2 );
+
+% The moduli are valid in two limits - minimizing or maximizing the
+% anisotropic difference elastic energy. Think of it as the maximum
+% positive deviations from the reference state or the maximum negative
+% deviations from the reference state. These are either in the positive
+% definite or negative definite regime of the matrix R. Here is a test of
+% the properties of R.  A +1 is returned if positive definite, a -1 is
+% retunred if negative definite.  0 returned otherwise.
+
+D = eig(R);
+if all(D>0)
+  minmax=1;
+elseif all(Dy0)
+  minmax=-1;
+else
+  minmax = 0;
+end
+
+end
+
+
+
+end

TensorAnalysis/@complianceTensor/VRHBulkModulus.m

@@ -0,0 +1,29 @@
+function [KV,KR,KVRH] = VRHBulkModulus(S)
+% Voigt-Reuss-Hill elastic bulk moduli
+%
+% Syntax
+%
+%   [KV,KR,KVRH] = VRHBulkModulus(S)
+%
+% Input
+%  S - @complianceTensor
+%
+% Output
+%  KV - bulk modulus Voigt average, upper bound
+%  KR - bulk modulus Reuss average, lower bound
+%  KVRH - bulk modulus Voigt Reuss Hill average, 
+
+% compute stifness tensor as 6x6 matrices
+C = matrix(inv(S),'voigt');
+S = matrix(S,'voigt');
+
+% KVoigt = ((C(1,1)+C(2,2)+C(3,3))+(2*(C(1,2)+C(2,3)+C(3,1))))/9
+KV = mean(C,[1,2]);
+
+% KReuss = 1/((S(1,1)+S(2,2)+S(3,3))+2*(S(1,2)+S(2,3)+S(3,1))) 
+KR = 1./sum(S,[1,2]);
+
+% Voigt Reuss Hill average
+KVRH = 0.5.*(KV + KR);
+
+end

TensorAnalysis/@complianceTensor/VRHShearModulus.m

@@ -0,0 +1,33 @@
+function [GV,GR,GVRH] = VRHShearModulus(S)
+% Voigt-Reuss-Hill elastic shear moduli
+%
+% Syntax
+%
+%   [GV,GR,GVRH] = VRHShearModulus(S)
+%
+% Input
+%  S - @complianceTensor
+%
+% Output
+%  GV - shear modulus Voigt average, upper bound
+%  GR - shear modulus Reuss average, lower bound
+%  GVRH - shear modulus Voigt Reuss Hill average, 
+
+% compute stifness tensor as 6x6 matrices
+C = matrix(inv(S),'voigt');
+S = matrix(S,'voigt');
+
+% GV = ((C(1,1)+C(2,2)+C(3,3))-(C(1,2)+C(2,3)+C(3,1))+3*(C(4,4)+C(5,5)+C(6,6)))/15
+GV= ((C(1,1,:) + C(2,2) + C(3,3,:)) ...
+  - (C(1,2,:) + C(2,3,:) + C(3,1,:)) ...
+  + 3 * (C(4,4,:) + C(5,5,:) + C(6,6,:))) ./ 15;
+
+% GR = 15/(4*(S(1,1)+S(2,2)+S(3,3))-4*(S(1,2)+S(2,3)+S(3,1))+3*(S(4,4)+S(5,5)+S(6,6)))
+GR = 15 ./ (4 * (S(1,1,:) + S(2,2,:) + S(3,3,:)) ...
+  - 4 * (S(1,2,:) + S(2,3,:) + S(3,1,:)) ...
+  + 3 * (S(4,4,:) + S(5,5,:) + S(6,6,:)));
+
+% Voigt Reuss Hill average
+GVRH = 0.5.*(GV + GR);
+
+end

TensorAnalysis/@stiffnessTensor/HashinShtrikmanModulus.m

@@ -0,0 +1,79 @@
+function [minmax, khs, ghs] = HashinShtrikmanModulus(C,K0,G0)
+% Hashin Shtrikman moduli 
+%
+% The equations are from Peselnick and Meister 1965 through Watt and
+% Peselnick 1980 which are based on Hashin and Shtrikman 1963 JMB 8/2013
+%
+% Syntax:
+%
+%   [minmax, khs, ghs] = HashinShtrikmanModulus(C,K0,G0)
+%
+% Input
+%  C  - @stiffnessTensor
+%  K0 - bulk modulus of the reference isotropic material
+%  G0 - shear modulus of the reference isotropic material
+%
+% Output
+%  khs - Hashin Shtrikman bulk modulus
+%  ghs - Hashin Shtrikman shear modulus
+%  minmax - +1 positive definite, -1 for negative definite R
+%
+
+% ensure K0 and G0 have same size
+if numel(G0) == 1, G0 = repmat(G0,size(K0)); end
+
+% the isotropic compliances elements  (eq. 10-11 Watt Peselnick 1980)
+alpha = -3 ./ (3*K0 + 4*G0);
+beta = -3 * (K0 + 2*G0) ./ ( 5*G0 .* (3*K0 + 4*G0) );
+gamma = (alpha - 3 * beta) ./ 9;
+
+% set up isotropic elastic matrix 
+C0 = matrix(2 * G0 .* stiffnessTensor.eye,'Voigt');
+
+C0(1:3,1:3,:,:) = C0(1:3,1:3,:,:) + shiftdim(K0 - 2/3 * G0,-2);
+
+% take difference between anisotropic and isotropic matrices and take
+% inverse to get the residual compliance matrix (eq 4 & 8 Watt Peselnick
+% 1980)
+R = C - C0;
+H = inv(R);
+
+% Subtract the isotropic compliances from the residual compliance matrix
+% and invert back to moduli space.  The matrix B is the moduli matrix that
+% when multiplied by the strain difference between the reference isotropic
+% response and the isotropic response of the actual material gives the
+% average stres <pij> (eq 15 - 16 Watt Peselnick 1980) note identity matrix
+% is Iinv
+A = matrix(H - beta .* complianceTensor.eye,'Voigt');
+A(1:3,1:3,:,:) = A(1:3,1:3,:,:) - shiftdim(gamma,-2);
+B = A;
+for k = 1:numel(A)/36
+  B(:,:,k) = inv(A(:,:,k));
+end
+
+% now perform the averaging of the B matrix to get two numbers (B1 and B2)
+% that are related to the two isotropic moduli - perturbations of the
+% reference values (eq 21 and 22 Watt Peselnick 1980)
+sB1 = sum(B(1:3,1:3,:,:),[1,2]);
+
+sB2 = B(1,1,:,:) + B(2,2,:,:) + B(3,3,:,:) + 2 * (B(4,4,:,:) + B(5,5,:,:) + B(6,6,:,:));
+B1 = reshape(2*sB1 - sB2,size(alpha)) / 15;
+B2 = reshape(3*sB2 - sB1,size(alpha)) / 30;
+
+% The Hashin-Shtrikman moduli are then calculated as changes from the
+% reference isotropic body. (eq 25 & 27 Watt and Peselnick 1980)
+khs = K0 + ( 3 * B1 + 2 * B2 ) ./ ( 3 + alpha .* (3 * B1 + 2 * B2) );
+ghs = G0 + B2 ./ ( 1 + 2 * beta .* B2 );
+
+% The moduli are valid in two limits - minimizing or maximizing the
+% anisotropic difference elastic energy. Think of it as the maximum
+% positive deviations from the reference state or the maximum negative
+% deviations from the reference state. These are either in the positive
+% definite or negative definite regime of the matrix R. Here is a test of
+% the properties of R.  A +1 is returned if positive definite, a -1 is
+% retunred if negative definite.  0 returned otherwise.
+
+minmax = zeros(size(R));
+D = eig(R);
+minmax(all(D>0,1)) = 1;
+minmax(all(D<0,1)) = -1;

TensorAnalysis/@stiffnessTensor/bulkModulus.m

@@ -0,0 +1,19 @@
+function [KV,KR,KVRH] = bulkModulus(C)
+% Voigt-Reuss-Hill elastic bulk moduli
+%
+% Syntax
+%
+%   [KV,KR,KVRH] = bulkModulus(S)
+%
+% Input
+%  S - @complianceTensor
+%
+% Output
+%  KV - bulk modulus Voigt average, upper bound
+%  KR - bulk modulus Reuss average, lower bound
+%  KVRH -bulk modulus Voigt Reuss Hill average, 
+
+
+[KV,KR,KVRH] = bulkModulus(inv(C));
+
+end

tools/math_tools/findJumps.m

@@ -0,0 +1,28 @@
+function x = findJumps(fun,xmin,xmax)
+% find all jumps of a monotonsly increasing function
+%
+% Syntax
+%
+% Input
+%
+% Output
+%
+
+v = fun([xmin,xmax]);
+m = (xmin + xmax)./2;
+
+if v(1) == v(2)
+
+  x = [];
+
+elseif xmax - xmin < 1e-5
+
+  x = m;
+
+else
+
+  x = [findJumps(fun,xmin,m),findJumps(fun,m,xmax)];
+
+end
+
+end