%% 线性粘弹性元件模型
% 粘弹性元件包括虎克弹簧和牛顿粘壶。它们分别代表着以下两个本构关系（以简单剪切的表观应力\sigma\left(t\right)和应变\gamma\left(t\right)为例）：
% 虎克弹簧：\sigma\left(t\right)=E\gamma\left(t\right)；
% 牛顿粘壶：\sigma\left(t\right)=\eta\dot{\gamma}\left(t\right)
% 它们的线性组合方式包括“串联”和“并联”。串联的两个元件，应力相等且等于总应力，应变相加等于总应变；并联的两个元件，应力相加等于总应力，应变相等且等于总应变。根据这个原则，我们可以写出任何线性粘弹性元件模型的本构关系。
%
% 图为各类线性粘弹性元件模型示意图
% 例如Maxwell模型。设其虎克弹簧的本构关系是\sigma_1\left(t\right)=E\gamma_1\left(t\right)，其牛顿粘壶的本构关系是\sigma_2\left(t\right)=\eta\dot{\gamma}_2\left(t\right)。则有\sigma\left(t\right)=\sigma_1\left(t\right)=\sigma_2\left(t\right)，\gamma\left(t\right)=\gamma_1\left(t\right)+\gamma_2\left(t\right)。由此有\dot{\gamma}\left(t\right)=\dot{\gamma}_1\left(t\right)+\dot{\gamma}_2\left(t\right)，把相应的本构关系（该求导的求导）代入该式即得到Maxwell模型的本构关系
% \dot{\gamma}\left(t\right)=E^{-1}\dot{\sigma}\left(t\right)+\eta\sigma\left(t\right)
% 一般令\tau\equiv\eta/E，称为Maxwell模型的松弛时间，则Maxwell模型的本构关系整理成更常见形式：
% Maxwell模型：\eta\dot{\gamma}\left(t\right)=\tau\dot{\sigma}\left(t\right)+\sigma\left(t\right)
% 用相同的方法可以推导出Kelvin-Voigt模型的本构关系：
% Kelvin-Voigt模型：E\gamma\left(t\right)+\tau E\dot{\gamma}\left(t\right)=\sigma\left(t\right)
% 其中的E和\tau\equiv\eta/E都是Kelvin-Voigt模型中的相应元件参数。
%% 1. 蠕变
% 蠕变（creep）测试是迫使材料形成如下的阶跃应力
% \sigma\left(t\right)=\sigma_0H\left(t\right)
% 考察材料的应变\gamma\left(t\right)，其中H\left(t\right) 是单位阶跃函数。我们还定义蠕变柔量（creep compliance）：
% J\left(t\right)\equiv\frac{\gamma\left(t\right)}{\sigma_0}
% 虎克弹簧的蠕变响应是
% \gamma\left(t\right)=\sigma_0E^{-1}H\left(t\right)
% 蠕变柔量是
% J\left(t\right)=E^{-1}H\left(t\right)
% 牛顿粘壶的蠕变响应是
% \gamma\left(t\right)=\int_{-\infty}^t\dot{\gamma}\left(t^\prime\right)\mathrm{d}t^\prime=\sigma_0\eta^{-1}\int_{-\infty}^tH\left(t^\prime\right)\mathrm{d}t^\prime=\sigma_0E^{-1}\frac{t}{\tau}H\left(t\right)
% Maxwell模型的蠕变响应：
% \gamma\left(t\right)=\int_{-\infty}^t\dot{\gamma}\left(t^\prime\right)\mathrm{d}t^\prime=\sigma_0E^{-1}\left(1+\frac{t}{\tau}\right)H\left(t\right)
% Kelvin-Voigt模型的蠕变响应，需要把\sigma\left(t\right)=\sigma_0H\left(t\right)代入其本构方程中，然后解关于\gamma\left(t\right)的常微分方程，其初条件是\gamma\left(0\right)=0，结果是：
% \gamma\left(t\right)=\sigma_0E^{-1}\left(1-e^{-t/\tau}\right)H\left(t\right)

log_sigma_0 =-0.55; % exponent for sigma_0
log_E =-0.48; % exponent for E
log_tau = -0.6; % expoenent for tau
blnCompliance=false;

% recover the values of the parameters from their exponents
sigma_0=10.^log_sigma_0;
E=10.^log_E;
tau=10.^log_tau;

eta=E.*tau;

% time vector
t=logspace(-3,2,1000)';

% strain of each model
gamma_Hookian=sigma_0.*E.^(-1).*heaviside(t);
gamma_Newtonian=sigma_0.*eta.^(-1).*t.*heaviside(t);
gamma_Maxwell=sigma_0.*(E.^(-1)+eta.^(-1).*t).*heaviside(t);
gamma_Kelvin_Voigt=sigma_0.*E.^(-1).*(1-exp(-t./tau));

% prepare data to plot
plot_data=[gamma_Hookian,gamma_Newtonian,gamma_Maxwell,gamma_Kelvin_Voigt];
if blnCompliance
    plot_data=plot_data./sigma_0;
end

% plot handle and properties 
h=loglog(t,plot_data);

h(1).LineWidth=1;
h(2).LineWidth=1;
h(3).LineWidth=1;
h(4).LineWidth=1;

% axes handle and properties

a=gca;
a.LineWidth=1;
a.XLim=[1e-3,1e2];
a.YLim=[1e-4,1e2];
a.XLabel.String="\itt\rm (s)";

if blnCompliance
    a.Title.String="剪切柔量";
    a.YLabel.String="\itJ\rm(\itt\rm)";
else
    a.Title.String="剪切应变";
    a.YLabel.String="\it\gamma\rm(\itt\rm)";
end

% legend handle and properties
l=legend;
l.String={'虎克弹簧','牛顿粘壶','Maxwell模型','Kelvin-Voigt模型'};
l.Location="best";

%
%% 2. 应力松弛
% 应力松弛（stress relaxation）测试是迫使材料形成如下的阶跃应变
% \gamma\left(t\right)=\gamma_0H\left(t\right)
% 考察材料的应力\sigma\left(t\right)。我们还定义应力松弛模量（stress relaxatxion modulus）
% G\left(t\right)\equiv\frac{\sigma\left(t\right)}{\gamma_0}
% 虎克弹簧的应力松弛响应是
% \sigma\left(t\right)=\gamma_0EH\left(t\right)
% 牛顿粘壶的应力松弛响应是
% \sigma\left(t\right)=\gamma_0\eta\delta\left(t\right)
% 其中\delta\left(t\right)是单位脉冲函数。
% Maxwell模型的应力松弛响应，需要将\gamma\left(t\right)=\gamma_0H\left(t\right)代入其本构方程，然后解关于\sigma\left(t\right)的常微分方程，其初条件是\sigma\left(0\right)=0，结果是：
% \sigma\left(t\right)=\gamma_0\frac{\eta}{\tau}e^{-t/\tau}\left[H\left(t\right)-H\left(0\right)\right]
% Kevin-Voigt模型的应力松弛响应是
% \sigma\left(t\right)=E\gamma_0H\left(t\right)+\eta\gamma_0\delta\left(t\right)

log_gamma_0 =-0.58; % exponent for sigma_0
log_E =0.42; % exponent for E
log_tau = 1.3; % expoenent for tau
blnModulus=false;

% recover the values of the parameters from their exponents
gamma_0=10.^log_gamma_0;
E=10.^log_E;
tau=10.^log_tau;

eta=E.*tau;

% time vector
t=logspace(-3,2,1000)';

% strain of each model
sigma_Hookian=gamma_0.*E.*heaviside(t);
sigma_Newtonian=gamma_0.*eta.*dirac(t);
sigma_Maxwell=gamma_0.*(eta./tau).*exp(-t./tau).*(heaviside(t)-heaviside(0));
sigma_Kelvin_Voigt=gamma_0.*E.*heaviside(t)+gamma_0.*eta.*dirac(t);


% prepare data to plot
plot_data=[sigma_Hookian,sigma_Newtonian,sigma_Maxwell,sigma_Kelvin_Voigt];
if blnModulus
    plot_data=plot_data./gamma_0;
end

% plot handle and properties 
h=loglog(t,plot_data);

h(1).LineWidth=1;
h(2).LineWidth=1;
h(3).LineWidth=1;
h(4).LineWidth=1;

% axes handle and properties

a=gca;
a.LineWidth=1;
a.XLim=[1e-3,1e2];
a.YLim=[1e-4,1e2];
a.XLabel.String="\itt\rm (s)";

if blnModulus
    a.Title.String="剪切模量";
    a.YLabel.String="\itG\rm(\itt\rm)";
else
    a.Title.String="剪切应力";
    a.YLabel.String="\it\sigma\rm(\itt\rm)";
end

% legend handle and properties
l=legend;
l.String={'虎克弹簧','牛顿粘壶','Maxwell模型','Kelvin-Voigt模型'};
l.Location="best";

%
%