(* Parameters *)
Y = 100 ;(* Young's modulus (GPa) *)
v = 0.22 ;(* Poisson's ratio *)
G = 27 ;(* Shear modulus (GPa) *)
a =10*(10^(-9)); (* Contact Radius (m) *)
p ; (* Mean Contact Pressure (GPa) *)

(* Displacement *)

u[r_,z_]:= 0.5*(((r^2)+(z^2)-(a^2))+(((((((r^2)+(z^2)-(a^2)))^2)+(((2 a z))^2)))^(1/2)));
σr[r_,z_,u_,p_]:=((3p)/2)(((1-2v)/3) ((a^2)/(r^2))(1-((((z/(Sqrt[u]))))^3))+((((z/(Sqrt[u]))))^3)(((a^2)u)/((u^2)+(a^2)(z^2)))+(z/(Sqrt[u]))(u ((1-v)/((a^2)+u))+(1+v)((Sqrt[u])/a)ArcTan[(a/(Sqrt[u]))]-2));
σθ[r_,z_,u_,p_]:=-((3p)/2)(((1-2v)/3) ((a^2)/(r^2))(1-((((z/(Sqrt[u]))))^3))+(z/(Sqrt[u]))(2v+u ((1-v)/((a^2)+u))-(1+v)((Sqrt[u])/a)ArcTan[(a/(Sqrt[u]))]));
σz[r_,z_,u_,p_]:=-((3p)/2)((((z/(Sqrt[u]))))^3)((((a^2)u)/((u^2)+(a^2)(z^2))));
τrz[r_,z_,u_,p_]:=-((3p)/2)((((z^2)r)/((u^2)+(a^2)(z^2))))((((a^(2 ))(Sqrt[u]))/((a^2)+u)));

(* Strain-Stress *)
ϵr[r_,z_,p_] := 1/Y (σr[r,z,u[r,z],p]-v (σθ[r,z,u[r,z],p] + σz[r,z,u[r,z],p]));
ϵθ[r_,z_,p_] := 1/Y (σθ[r,z,u[r,z],p] -v (σz[r,z,u[r,z],p] + σr[r,z,u[r,z],p]));
ϵz[r_,z_,p_] := 1/Y (σz[r,z,u[r,z],p] -v (σr[r,z,u[r,z],p] + σθ[r,z,u[r,z],p]));

VectorPlot[{{σr[r,z,u[r,z],100],σz[r,z,u[r,z],100]},{ϵr[r,z,100],ϵz[r,z,100]}},{r,0.001*(10^(-9)),10*(10^(-9))},{z,0.001*(10^(-9)),10*(10^(-9))}, VectorPoints->Fine,VectorColorFunction->"Rainbow",VectorColorFunctionScaling->False]

(* Strain and strain gradient with repect to z-axis *)
Plot[{Evaluate[ϵz[0.1*(10^(-9)),z,10]],Evaluate[ϵz[0.1*(10^(-9)),z,20]],Evaluate[ϵz[0.1*(10^(-9)),z,50]],Evaluate[ϵz[0.1*(10^(-9)),z,100]]},{z,0,0.1*(10^(-6))}, PlotRange->{-1.0,0}]
Plot[{Evaluate[D[ϵz[0.1*(10^(-9)),z,4],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,8],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,12],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,16],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,20],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,24],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,28],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,32],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,36],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,40],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,44],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,48],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,52],z]]},{z,0,0.06*(10^(-6))}, PlotRange->{-8*(10^6),40*(10^6)}]

(* Strain distribution *)
ContourPlot[Evaluate[ϵr[r,z,50]],{r,0.1*(10^(-9)),0.05*(10^(-6))},{z,0.1*(10^(-9)),0.05*(10^(-6))},PlotRange-> {-0.5,0.1},ContourLabels->True,Contours->30]
ContourPlot[Evaluate[ϵθ[r,z,50]],{r,0.1*(10^(-9)),0.05*(10^(-6))},{z,0.1*(10^(-9)),0.05*(10^(-6))},PlotRange-> {-0.2,0.08},ContourLabels->True,Contours->30]
ContourPlot[Evaluate[ϵz[r,z,50]],{r,0.1*(10^(-9)),0.05*(10^(-6))},{z,0.1*(10^(-9)),0.05*(10^(-6))},PlotRange-> {-0.7,0.01},ContourLabels->True,Contours->30]

(* strain gradients with respect to z *)
ContourPlot[Evaluate[D[ϵr[r,z,50],z]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-200*(10^6),200*(10^6)},ContourLabels->True,Contours->50]
ContourPlot[Evaluate[D[ϵθ[r,z,50],z]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-10*(10^6),200*(10^6)},ContourLabels->True,Contours->50]
ContourPlot[Evaluate[D[ϵz[r,z,50],z]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-100*(10^6),100*(10^6)},ContourLabels->True,Contours->50]

(* strain gradients with respect to r *)
ContourPlot[Evaluate[D[ϵr[r,z,50],r]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-60*(10^6),200*(10^6)},Contours->30]
ContourPlot[Evaluate[D[ϵθ[r,z,50],r]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-20*(10^6),50*(10^6)},Contours->30]
ContourPlot[Evaluate[D[ϵz[r,z,50],r]],{r,0.1*(10^(-9)),0.03*(10^(-6))},{z,0.1*(10^(-9)),0.03*(10^(-6))},PlotRange->{-20*(10^6),500*(10^6)},Contours->30]

Plot[{Evaluate[D[ϵz[0.1*(10^(-9)),z,10],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,20],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,50],z]],Evaluate[D[ϵz[0.1*(10^(-9)),z,100],z]]},{z,0,0.1*(10^(-6))}, PlotRange->{-8*(10^6),1*(10^6)}]

Solve[D[ϵz[0.1*(10^(-9)),z,10],z]==0,z]