%%% magnetic composite constitutive equation and its response to the constant magnetic field at z direction 
% due to the placement of the magnets in relation to the sample, it is
% assumed that the magnetic intencity change through the thickness of
% sample by this relation H=H0*(1-(delta*(z/t))) in which t is half of
% thickness (because of the symmetry of the problem) and delta (dissipation
% of magnetic field through the thickness)is considered to be 80% according
% to the measurment 
clear
clc

%% Iron particles material properties
Ei=211*(10^9); mioi=0.29;

%% PLA epoxy mechanical property
% pla mechanical properties are highly affected by the tempreture
Ep=3.1*(10^9); % at room degree of celsius
miop=0.33;  

%% calculation of the magnetic field induced stress 
m0=4*pi*(10^(-7));%magnetic permeability
x1=200000; %magnetic volume susceptibility of iron
mr1=x1+1;% magnetic relative permeability of iron
mn1=m0*mr1; %magnetic permeability of iron
a=input("please enter magnetic flux density (mT)\")

B0= a*(10^(-3));% Magnetic flux density at the surface (measured by the equipment= 65 mT)
H0= B0/mn1; %magnetic intensity at the surfac
vf=0.11; % iron particles volume fraction 
d=10*(10^(-6)); % averaged particle size
delta=0; % magnetic penetration through the thicknss
t= 7*(10^(-3)); %half of the thickness
A=(50*15)*(10^(-6)); % area exposed to magnetic field
V_particle=pi*(d^3)/6;
V_t=A*t;
Nt=(vf*V_t)/V_particle;
%% distribution of iron particles inside the composite 
% Considering spherical iron particles with the same diameter, 
% The sample is divided into smaller volumes, each with a thickness of dz. 
% Within each of these volumes, we make the assumption that every iron 
% particle is situated at a distance of z from the exposed surface.
dz=50*(10^(-6));
nz=round(t/dz);
Nuc=(vf*dz*A)/V_particle; %number of particles at each seperated volume

z1=dz/2;
f=zeros(nz,1);
for n=1:nz
    if n==1
        z=z1;
    else
    z=z+dz;
    end
H=H0*(1-(delta*(z/t)));
M=x1*H; %magnitude of magnetization
I=M*(pi*d);
f(n)=Nuc*I*pi*d*mn1*H;%Lorentz force acting on each iron particle
end
f;
ft=sum(f);
zigma_t=ft/A;
%% iron and pla compliant matrix
Si1=[1 -mioi -mioi;-mioi 1 -mioi;-mioi -mioi 1];
Si2=[1+mioi 0 0;0 1+mioi 0;0 0 1+mioi];
Si=(1/Ei)*[Si1 zeros(3,3);zeros(3,3) Si2]; %iron

Sp1=[1 -miop -miop;-miop 1 -miop;-miop -miop 1];
Sp2=[1+miop 0 0;0 1+miop 0;0 0 1+miop];
Sp=(1/Ep)*[Sp1 zeros(3,3);zeros(3,3) Sp2]; %pla

mio_eff=vf*mioi+(1-vf)*miop+(vf*(mioi-miop))/(1-vf*(mioi-miop));

%% micromechanical model---Mori-Tanaka
% Esheby tensor for spherical inhemogenities and perfect interface ---equivalent eigenstrain problem
ES=zeros(3,3,3,3);
for i=1:3
    for j=1:3
        for k=1:3
            for l=1:3
ES(i,j,k,l)=((5*miop-1)/(15*(1-miop)))*kroneckerDelta(i,j)*kroneckerDelta(k,l)+((4-5*miop)/(15*(1-miop)))*(kroneckerDelta(i,k)*kroneckerDelta(j,l)+kroneckerDelta(i,l)*kroneckerDelta(j,k));
            end
        end
    end
end
%Mandel notation for transforminng the fourth order tensor into a 6*6 matrix
ESm=[ES(1,1,1,1) ES(1,1,2,2) ES(1,1,3,3) sqrt(2)*ES(1,1,2,3) sqrt(2)*ES(1,1,3,1) sqrt(2)*ES(1,1,1,2);
    ES(2,2,1,1) ES(2,2,2,2) ES(2,2,3,3) sqrt(2)*ES(2,2,2,3) sqrt(2)*ES(2,2,3,1) sqrt(2)*ES(2,2,1,2);
    ES(3,3,1,1) ES(3,3,2,2) ES(3,3,3,3) sqrt(2)*ES(3,3,2,3) sqrt(2)*ES(3,3,3,1) sqrt(2)*ES(3,3,1,2);
    sqrt(2)*ES(2,3,1,1) sqrt(2)*ES(2,3,2,2) sqrt(2)*ES(2,3,3,3) 2*ES(2,3,2,3) 2*ES(2,3,3,1) 2*ES(2,3,1,2);
    sqrt(2)*ES(3,1,1,1) sqrt(2)*ES(3,1,2,2) sqrt(2)*ES(3,1,3,3) 2*ES(3,1,2,3) 2*ES(3,1,3,1) 2*ES(3,1,1,2);
    sqrt(2)*ES(1,2,1,1) sqrt(2)*ES(1,2,2,2) sqrt(2)*ES(1,2,3,3) 2*ES(1,2,2,3) 2*ES(1,2,3,1) 2*ES(1,2,1,2)];
cm=0.83; cr=0.17;
Br_dil=inv(eye(6)+ESm*inv(Sp)*(Si-Sp));
Br_MT= Br_dil*inv(cm*eye(6)+cr*Br_dil);
Bm_MT=inv(cm*eye(6)+cr*Br_dil);
St=cr*Si*Br_MT+cm*Sp*Bm_MT;
b=input("please enter tensile stress (Mpa)\")
tzigma=b*(10^6);%Tensile stress

%stress tensor
tzigma_n=tzigma;
zigma_tn=zigma_t-mio_eff*tzigma;
zigma=[tzigma_n; 0; zigma_tn; 0; 0; 0]; 
epsi=St*zigma   %strain tensor
epsi1=epsi(1)*100 %tensile strain %
u3=epsi(3)*t *(10^3); % deformation at z-direction (in mm)



    

 