Code used to predict strain-dependent relaxation

Import Data

Relax1=Import["INSERT FILEPATH HERE"];
Relax1Time=Drop[Relax1[[All,1]],2]-Relax1[[3,1]];
Relax1Strain=Drop[Relax1[[All,10]],2];
Relax1Stress=Drop[Relax1[[All,11]],2];
Relax2=Import["INSERT FILEPATH HERE"];
Relax2Time=Drop[Relax2[[All,1]],2]-Relax2[[3,1]];
Relax2Strain=Drop[Relax2[[All,10]],2];
Relax2Stress=Drop[Relax2[[All,11]],2];
Relax3=Import["INSERT FILEPATH HERE"];
Relax3Time=Drop[Relax3[[All,1]],2]-Relax3[[3,1]];
Relax3Strain=Drop[Relax3[[All,10]],2];
Relax3Stress=Drop[Relax3[[All,11]],2];

Calculate positions of peak stresses at 6% strain

Position[Relax1Stress,Max[Relax1Stress]]
Position[Relax2Stress,Max[Relax2Stress]]
Position[Relax3Stress,Max[Relax3Stress]]

Extract part of curve from peak stress onwards

Relax1StressTime6=Transpose[{Relax1Time[[60161;;]]-Relax1Time[[60161]],Relax1Stress[[60161;;]]/Relax1Stress[[60161]]}];
Relax1StrainTime6=Transpose[{Relax1Time[[60161;;]]-Relax1Time[[60161]],Relax1Strain[[60161;;]]/Relax1Stress[[60161]]}];
Relax2StressTime6=Transpose[{Relax2Time[[60159;;]]-Relax2Time[[60159]],Relax2Stress[[60159;;]]/Relax2Stress[[60159]]}];
Relax2StrainTime6=Transpose[{Relax2Time[[60159;;]]-Relax2Time[[60159]],Relax2Strain[[60159;;]]/Relax2Stress[[60159]]}];
Relax3StressTime6=Transpose[{Relax3Time[[60174;;]]-Relax3Time[[60174]],Relax3Stress[[60174;;]]/Relax3Stress[[60174]]}];
Relax3StrainTime6=Transpose[{Relax3Time[[60174;;]]-Relax3Time[[60174]],Relax3Strain[[60174;;]]/Relax3Stress[[60174]]}];

Calculate positions of peak stresses at 4% strain

Position[Relax1Stress[[1;;50000]],Max[Relax1Stress[[1;;50000]]]]
Position[Relax2Stress[[1;;50000]],Max[Relax2Stress[[1;;50000]]]]
Position[Relax3Stress[[1;;50000]],Max[Relax3Stress[[1;;50000]]]]

Extract part of curve from peak stress onwards

Relax1StressTime4=Transpose[{Relax1Time[[30098;;]]-Relax1Time[[30098]],Relax1Stress[[30098;;]]/Relax1Stress[[30098]]}];
Relax1StrainTime4=Transpose[{Relax1Time[[30098;;]]-Relax1Time[[30098]],Relax1Strain[[30098;;]]/Relax1Stress[[30098]]}];
Relax2StressTime4=Transpose[{Relax2Time[[30095;;]]-Relax2Time[[30095]],Relax2Stress[[30095;;]]/Relax2Stress[[30095]]}];
Relax2StrainTime4=Transpose[{Relax2Time[[30095;;]]-Relax2Time[[30095]],Relax2Strain[[30095;;]]/Relax2Stress[[30095]]}];
Relax3StressTime4=Transpose[{Relax3Time[[30106;;]]-Relax3Time[[30106]],Relax3Stress[[30106;;]]/Relax3Stress[[30106]]}];
Relax3StrainTime4=Transpose[{Relax3Time[[30106;;]]-Relax3Time[[30106]],Relax3Strain[[30106;;]]/Relax3Stress[[30106]]}];

Calculate positions of peak stresses at 2% strain

Position[Relax1Stress[[1;;20000]],Max[Relax1Stress[[1;;20000]]]]
Position[Relax2Stress[[1;;20000]],Max[Relax2Stress[[1;;20000]]]]
Position[Relax3Stress[[1;;20000]],Max[Relax3Stress[[1;;20000]]]]

Extract part of curve from peak stress onwards

Relax1StressTime2=Transpose[{Relax1Time[[58;;]]-Relax1Time[[58]],Relax1Stress[[58;;]]/Relax1Stress[[58]]}];
Relax1StrainTime2=Transpose[{Relax1Time[[58;;]]-Relax1Time[[58]],Relax1Strain[[58;;]]/Relax1Stress[[58]]}];
Relax2StressTime2=Transpose[{Relax2Time[[16;;]]-Relax2Time[[16]],Relax2Stress[[16;;]]/Relax2Stress[[16]]}];
Relax2StrainTime2=Transpose[{Relax2Time[[16;;]]-Relax2Time[[16]],Relax2Strain[[16;;]]/Relax2Stress[[16]]}];
Relax3StressTime2=Transpose[{Relax3Time[[6;;]]-Relax3Time[[6]],Relax3Stress[[6;;]]/Relax3Stress[[6]]}];
Relax3StrainTime2=Transpose[{Relax3Time[[6;;]]-Relax3Time[[6]],Relax3Strain[[6;;]]/Relax3Stress[[6]]}];

6% strain curves

ListPlot[{Relax1StressTime6,Relax2StressTime6,Relax3StressTime6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

4% strain curves

ListPlot[{Relax1StressTime4,Relax2StressTime4,Relax3StressTime4},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

2% strain curves

ListPlot[{Relax1StressTime2,Relax2StressTime2,Relax3StressTime2},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

Sample 1 at 2%, 4% and 6% strain

ListPlot[{Relax1StressTime2,Relax1StressTime4,Relax1StressTime6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

Sample 2 at 2%, 4% and 6% strain

ListPlot[{Relax2StressTime2,Relax2StressTime4,Relax2StressTime6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

Sample 3 at 2%, 4% and 6% strain

ListPlot[{Relax3StressTime2,Relax3StressTime4,Relax3StressTime6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

Calculate average responses across all three samples

RelaxStressTime2=Table[{Relax3StressTime2[[i,1]],Mean[{Relax1StressTime2[[i,2]],Relax2StressTime2[[i,2]],Relax3StressTime2[[i,2]]}]},{i,1,30000}];
RelaxStressTime4=Table[{Relax3StressTime4[[i,1]],Mean[{Relax1StressTime4[[i,2]],Relax2StressTime4[[i,2]],Relax3StressTime4[[i,2]]}]},{i,1,30000}];
RelaxStressTime6=Table[{Relax3StressTime6[[i,1]],Mean[{Relax1StressTime6[[i,2]],Relax2StressTime6[[i,2]],Relax3StressTime6[[i,2]]}]},{i,1,30000}];

ListPlot[{RelaxStressTime2,RelaxStressTime4,RelaxStressTime6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s}]

Calculate smoothed version of above

Transpose[{MovingAverage[RelaxStressTime2[[1;;500,1]],100],MovingAverage[RelaxStressTime2[[1;;500,2]],100]}];

Averaged2=Join[Transpose[{MovingAverage[RelaxStressTime2[[1;;500,1]],100],MovingAverage[RelaxStressTime2[[1;;500,2]],100]}],Transpose[{MovingAverage[RelaxStressTime2[[All,1]],1000],MovingAverage[RelaxStressTime2[[All,2]],1000]}]];
Averaged4=Join[Transpose[{MovingAverage[RelaxStressTime4[[1;;500,1]],100],MovingAverage[RelaxStressTime4[[1;;500,2]],100]}],Transpose[{MovingAverage[RelaxStressTime4[[All,1]],1000],MovingAverage[RelaxStressTime4[[All,2]],1000]}]];
Averaged6=Join[Transpose[{MovingAverage[RelaxStressTime6[[1;;500,1]],100],MovingAverage[RelaxStressTime6[[1;;500,2]],100]}],Transpose[{MovingAverage[RelaxStressTime6[[All,1]],1000],MovingAverage[RelaxStressTime6[[All,2]],1000]}]];

experimental=ListPlot[{Averaged2,Averaged4,Averaged6},Joined->True,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},PlotRange->{{0,300},{0,1}},AxesLabel->{t,s},AxesStyle->Thick,LabelStyle->Directive[20],Ticks->{{0,50,100,150,200,250},{0,0.2,0.4,0.6,0.8,1.0}}]

Theoretical relaxation tests

e[t_,emax_]:=If[t<0,0,If[t<0.01,t emax/0.01,emax]]

Plot[{e[t,0.02],e[t,0.04],e[t,0.06]},{t,0,300},AspectRatio->0.6,PlotStyle->{{Thick,Red},{Thick,Green},{Thick,Blue}},AxesLabel->{t,"e"},Ticks->True,AxesStyle->Thick,LabelStyle->Directive[20],PlotRange->{0,0.06}]

Relaxation function and its derivative

EE[t_]:=E0+E1 Exp[-t/τ1]+E2 Exp[-t/τ2]

DE[t_]:=Evaluate[D[EE[t],t]]

Relaxation function parameters

E0=3000;(*MPa*)
E1=950;(*MPa*)
E2=810;(*MPa*)
τ1=1.9;(*s*)
τ2=52;(*s*)

Fibril strain

ef[t_,ec1_,emax_]:=If[(e[t,emax]-ec1)/(1+ec1)<0,0,(e[t,emax]-ec1)/(1+ec1)]

Fibril stress

σf[t_,ec1_,emax_]:=If[(e[t,emax]-ec1)/(1+ec1)<0,0,EE[0](e[t,emax]-ec1)/(1+ec1)+NIntegrate[DE[t-s]ef[s,ec1,emax],{s,0,t}]]

Plot[σf[t,0.03,0.06]/σf[0.01,0.03,0.06],{t,0,300},PlotRange->{0,1}]

Matrix stress

EEm[t_]:=E∞m+(E0m-E∞m) E^(-t/τrm)

DEm[t_]:=Evaluate[D[EEm[t],t]]

Matrix constitutive parameter values

E0m=20;
E∞m=10;
τrm=1;

σm[t_,emax_]:=EEm[0] e[t,emax]+NIntegrate[DEm[t-s] e[s,emax],{s,0,t}]

σ02m=Table[{t,σm[t+0.01,0.02]/σm[0.01,0.02]},{t,0,100}];

σ04m=Table[{t,σm[t+0.01,0.04]/σm[0.01,0.04]},{t,0,100}];

σ06m=Table[{t,σm[t+0.01,0.06]/σm[0.01,0.06]},{t,0,100}];

ListPlot[{σ02m,σ04m,σ06m},Joined->True,PlotStyle->{{Thick,Black},{Thick,Red},{Thick,Blue},{Thick,Purple}},AxesLabel->{t,"s"},PlotRange->{0,1}]

Fascicle stress

Truncated Normal Distribution

NN[x_,μ_,σ_]:=PDF[TruncatedDistribution[{0,Infinity},NormalDistribution[μ,σ]],x](*Define critical strain distribution*)

Plot[NN[x,0.05,0.004],{x,0,0.2},PlotRange->All]

Note-range of integration is equal to range of maximum macroscopic strain

σF[t_,μ_,σ_,ϕ_,emax_]:=(1-ϕ) σm[t,emax]+ϕ NIntegrate[Hold[NN[ec1,μ,σ] σf[t,ec1,emax]],{ec1,0,emax},AccuracyGoal->8]

trange=Table[300/(10^10-1)(10^x-1),{x,0,10,0.1}];

a2=σF[0.01,0.05,0.004,0.8,0.02];
σ2=Quiet[Table[{t,σF[t,0.05,0.004,0.8,0.02]/a2},{t,trange}]];

a4=σF[0.01,0.05,0.004,0.8,0.04];
σ4=Quiet[Table[{t,σF[t,0.05,0.004,0.8,0.04]/a4},{t,trange}]];

a6=σF[0.01,0.05,0.004,0.8,0.06];
σ6=Quiet[Table[{t,σF[t,0.05,0.004,0.8,0.06]/a6},{t,trange}]];

predicted=ListPlot[{σ2,σ4,σ6},Joined->True,PlotStyle->{{Thick,Darker[Red],Dashed},{Thick,Darker[Green],Dashed},{Thick,Darker[Blue],Dashed}},AxesLabel->{t,s},PlotRange->{{0,300},{0,1}},AxesStyle->Thick,LabelStyle->Directive[20],Ticks->{{0,50,100,150,200,250},{0,0.2,0.4,0.6,0.8,1.0}}]

Show[experimental,predicted]