Hysteretic model
Hysteretic models are mathematical models capable of simulating the complex nonlinear behavior characterizing mechanical systems and materials used in different fields of engineering, such as aerospace, civil, and mechanical engineering. Some examples of mechanical systems and materials having hysteretic behavior are:
- materials, such as steel, reinforced concrete, wood;
- structural elements, such as steel, reinforced concrete, or wood joints;
- devices, such as seismic isolators[1] and dampers.
Hysteretic models may have a generalized displacement [math]\displaystyle{ u }[/math] as input variable and a generalized force [math]\displaystyle{ f }[/math] as output variable, or vice versa. In particular, in rate-independent hysteretic models, the output variable does not depend on the rate of variation of the input one.[2][3]
Rate-independent hysteretic models can be classified into four different categories depending on the type of equation that needs to be solved to compute the output variable:
- Algebraic models;
- Transcendental models;
- Differential models;
- Integral models.
Algebraic models
In algebraic models, the output variable is computed by solving algebraic equations.
Bilinear model
Model formulation
In the bilinear model formulated by Vaiana et al. (2018),[4] the generalized force at time t, representing the output variable, is evaluated as a function of the generalized displacement as follows:
[math]\displaystyle{ f_t=k_a\left(u_t-u_j\right)+k_b u_j+s_t f_0 \quad \text{when} \quad u_t s_t\lt u_j s_t, }[/math]
[math]\displaystyle{ f_t=k_b u_t+s_t f_0 \quad \text{when} \quad u_t s_t\gt u_j s_t, }[/math]
where [math]\displaystyle{ k_a, k_b, }[/math] and [math]\displaystyle{ f_0 }[/math] are the three model parameters to be calibrated from experimental or numerical tests, whereas [math]\displaystyle{ s_t }[/math] is the sign of the generalized velocity at time [math]\displaystyle{ t }[/math], that is, [math]\displaystyle{ s_t = \operatorname{sign}(\dot{u}_t) }[/math]. Furthermore, [math]\displaystyle{ u_0 }[/math]is an internal model parameter evaluated as:
[math]\displaystyle{ u_0=\frac{f_0}{k_a-k_b}, }[/math]
whereas [math]\displaystyle{ u_j }[/math] is the internal variable:
[math]\displaystyle{ u_j=\frac{k_a u_{t-\Delta t}+s_t f_0-f_{t-\Delta t}}{k_a-k_b} }[/math].
Hysteresis loop shapes
Figure 1.1 shows two different hysteresis loop shapes obtained by applying a sinusoidal generalized displacement having unit amplitude and frequency and simulated by adopting the Bilinear Model (BM) parameters listed in Table 1.1. File:Bilinear Model.tif
| [math]\displaystyle{ k_a }[/math] | [math]\displaystyle{ k_b }[/math] | [math]\displaystyle{ f_0 }[/math] | |
|---|---|---|---|
| (a) | 10.0 | 1.0 | 0.5 |
| (b) | 10.0 | -1.0 | 0.5 |
Matlab code
% =========================================================================================
% June 2020
% Bilinear Model Algorithm
% Nicolò Vaiana, Research Fellow in Structural Mechanics and Dynamics, PhD
% Department of Structures for Engineering and Architecture
% University of Naples Federico II
% via Claudio, 21 - 80124, Napoli
% =========================================================================================
clc; clear all; close all;
%% APPLIED DISPLACEMENT TIME HISTORY
dt = 0.001; %time step
t = 0:dt:1.50; %time interval
a0 = 1; %applied displacement amplitude
fr = 1; %applied displacement frequency
u = a0*sin((2*pi*fr)*t(1:length(t))); %applied displacement vector
v = 2*pi*fr*a0*cos((2*pi*fr)*t(1:length(t))); %applied velocity vector
n = length(u); %applied displacement vector length
%% 1. INITIAL SETTINGS
%1.1 Set the three model parameters
ka = 10.0; %model parameter
kb = 1.0; %model parameter
f0 = 0.5; %model parameter
%1.2 Compute the internal model parameters
u0 = f0/(ka-kb); %internal model parameter
%1.3 Initialize the generalized force vector
f = zeros(1,n);
%% 2. CALCULATIONS AT EACH TIME STEP
for i = 2:n
%2.1 Update the history variable
uj = (ka*u(i-1)+sign(v(i))*f0-f(i-1))/(ka-kb);
%2.2 Evaluate the generalized force at time t
if (sign(v(i))*uj)-2*u0 < sign(v(i))*u(i) && sign(v(i))*u(i) < sign(v(i))*uj
f(i) = ka*(u(i)-uj)+kb*uj+sign(v(i))*f0;
else
f(i) = kb*u(i)+sign(v(i))*f0;
end
end
%% PLOT
figure
plot(u,f,'k','linewidth',4)
set(gca,'FontSize',28)
set(gca,'FontName','Times New Roman')
xlabel('generalized displacement'), ylabel('generalized force')
grid
zoom off
Asymmetric bilinear model
Model formulation
In the asymmetric bilinear model formulated by Vaiana et al. (2020),[3] the generalized force at time t, representing the output variable, is evaluated as a function of the generalized displacement as follows:
[math]\displaystyle{ f_t=k_a\left(u_t-u_j\right)+k_b u_j+s_t f_0 \quad \text{when} \quad u_t s_t\lt u_j s_t, }[/math]
[math]\displaystyle{ f_t=k_b u_t+s_t f_0 \quad \text{when} \quad u_t s_t\gt u_j s_t, }[/math]
where [math]\displaystyle{ k_a=k_a^+(k_a^-), k_b=k_b^+(k_b^-), }[/math] and [math]\displaystyle{ f_0 = f_0^+ (f_0^-) }[/math] are the three model parameters of the generic loading (unloading) case to be calibrated from experimental or numerical tests, whereas [math]\displaystyle{ s_t }[/math] is the sign of the generalized velocity at time [math]\displaystyle{ t }[/math], that is, [math]\displaystyle{ s_t = \operatorname{sign}(\dot{u}_t) }[/math]. Furthermore, [math]\displaystyle{ u_0 }[/math] is an internal model parameter evaluated as:
[math]\displaystyle{ u_0=\frac{(k_b^+-k_b^-)u_j+f_0^++f_0^-}{2(k_a^+-k_b^-)} \quad\Biggl({\displaystyle {\frac {(k_{b}^{+}-k_{b}^{-})u_{j}+f_{0}^{+}+f_{0}^{-}}{2(k_{a}^{-}-k_{b}^{+})}}}\Biggr) \quad\text{when}\quad s_t\gt 0 \quad (s_t\lt 0), }[/math]
whereas [math]\displaystyle{ u_j }[/math] is the internal variable:
[math]\displaystyle{ u_j=\frac{k_a^+ u_{t-\Delta t}+s_t f_0^+-f_{t-\Delta t}}{k_a^+-k_b^+} \quad\Biggl({\displaystyle {\frac {k_a^- u_{t-\Delta t}+s_t f_0^--f_{t-\Delta t}}{k_a^--k_b^-}}}\Biggr) \quad\text{when}\quad s_t\gt 0 \quad (s_t\lt 0). }[/math]
Matlab code
% =========================================================================================
% February 2021
% Asymmetric Bilinear Model Algorithm
% Nicolò Vaiana, Research Fellow in Structural Mechanics and Dynamics, PhD
% Department of Structures for Engineering and Architecture
% University of Naples Federico II
% via Claudio, 21 - 80124, Napoli
% =========================================================================================
clc; clear all; close all;
%% APPLIED DISPLACEMENT TIME HISTORY
dt = 0.001; %time step
t = 0:dt:1.50; %time interval
a0 = 1; %applied displacement amplitude
fr = 1; %applied displacement frequency
u = a0*sin((2*pi*fr)*t(1:length(t))); %applied displacement vector
v = 2*pi*fr*a0*cos((2*pi*fr)*t(1:length(t))); %applied velocity vector
n = length(u); %applied displacement vector length
%% 1. INITIAL SETTINGS
%1.1 Set the six model parameters
kap = 5.0; %model parameter
kbp = 0.5; %model parameter
f0p = 0.75; %model parameter
kam = 15.0; %model parameter
kbm = 0.1; %model parameter
f0m = 0.25; %model parameter
%1.2 Initialize the generalized force vector
f = zeros(1,n);
%% 2. CALCULATIONS AT EACH TIME STEP
for i = 2:n
%2.1 Update the model parameters, the history variable, and the internal model parameter
if v(i)>0
ka = kap; kb = kbp; f0 = f0p;
else
ka = kam; kb = kbm; f0 = f0m;
end
uj = (ka*u(i-1)+sign(v(i))*f0-f(i-1))/(ka-kb);
if v(i)>0
u0 = ((kbp-kbm)*uj+f0p+f0m)/(2*(kap-kbm));
else
u0 = ((kbp-kbm)*uj+f0p+f0m)/(2*(kam-kbp));
end
%2.2 Evaluate the generalized force at time t
if (sign(v(i))*uj)-2*u0 < sign(v(i))*u(i) && sign(v(i))*u(i) < sign(v(i))*uj
f(i) = ka*(u(i)-uj)+kb*uj+sign(v(i))*f0;
else
f(i) = kb*u(i)+sign(v(i))*f0;
end
end
%% PLOT
figure
plot(u,f,'k','linewidth',4)
set(gca,'FontSize',28)
set(gca,'FontName','Times New Roman')
xlabel('generalized displacement'), ylabel('generalized force')
grid
zoom off
Animation
The following gif shows the nonlinear response of a single-degree-of-freedom (SDOF) mechanical system, with unit mass and asymmetric rate-independent hysteretic behavior, subjected to an external random force. To simulate its response, the following ABM parameters have been used: [math]\displaystyle{ k_a^+=50 \,\, (k_a^-=50), k_b^+=10 \,\, (k_b^-=1), f_0^+=1 \,\, (f_0^-=2) }[/math].
Algebraic model by Vaiana et al. (2019)
Model formulation
In the algebraic model developed by Vaiana et al. (2019),[5] the generalized force at time [math]\displaystyle{ t }[/math], representing the output variable, is evaluated as a function of the generalized displacement as follows:
[math]\displaystyle{ f_t=\beta_1 u_t^3 + \beta_2 u_t^5+k_b u_t+\left(k_a-k_b\right)\left[\frac{(1+s_t u_t-s_t u_j+2u_0)^{1-\alpha}}{s_t (1-\alpha)}-\frac{(1+2u_0)^{1-\alpha}}{s_t (1-\alpha)}\right]+s_t f_0 \quad \text{when} \quad u_t s_t\lt u_j s_t, }[/math]
[math]\displaystyle{ f_t=\beta_1 u_t^3+\beta_2 u_t^5 + k_b u_t+s_t f_0 \quad \text{when} \quad u_t s_t\gt u_j s_t, }[/math]
where [math]\displaystyle{ k_a, k_b, \alpha }[/math], [math]\displaystyle{ \beta_1 }[/math] and [math]\displaystyle{ \beta_2 }[/math] are the five model parameters to be calibrated from experimental or numerical tests, whereas [math]\displaystyle{ s_t }[/math] is the sign of the generalized velocity at time [math]\displaystyle{ t }[/math], that is, [math]\displaystyle{ s_t = \operatorname{sign}(\dot{u}_t) }[/math]. Furthermore, [math]\displaystyle{ u_0 }[/math]and [math]\displaystyle{ f_0 }[/math] are two internal model parameters evaluated as:
[math]\displaystyle{ u_0=\frac{1}{2}\left[\left(\frac{k_a-k_b}{10^{-20}}\right)^{1/\alpha}-1\right], }[/math]
[math]\displaystyle{ f_0=\frac{k_a-k_b}{2}\left[\frac{\left(1+2u_0\right)^{1-\alpha}-1}{1-\alpha}\right], }[/math]
whereas [math]\displaystyle{ u_j }[/math] is the internal variable:
[math]\displaystyle{ u_j=u_{t-\Delta t}+s_t(1+2u_0)-s_t\left\lbrace \frac{s_t(1-\alpha)}{k_a-k_b} \left[f_{t-\Delta t}-\beta_1 u_{t-\Delta t}^3 - \beta_2 u_{t-\Delta t}^5 - k_b u_{t-\Delta t}-s_t f_0+(k_a-k_b) \frac{(1+2u_0)^{1-\alpha}}{s_t(1-\alpha)}\right]\right\rbrace^{1/(1-\alpha)}. }[/math]
Hysteresis loop shapes
Figure 1.2 shows four different hysteresis loop shapes obtained by applying a sinusoidal generalized displacement having unit amplitude and frequency and simulated by adopting the Algebraic Model (AM) parameters listed in Table 1.2.File:Figure NV 3.tif
| [math]\displaystyle{ k_a }[/math] | [math]\displaystyle{ k_b }[/math] | [math]\displaystyle{ \alpha }[/math] | [math]\displaystyle{ \beta_1 }[/math] | [math]\displaystyle{ \beta_2 }[/math] | |
|---|---|---|---|---|---|
| (a) | 10.0 | 1.0 | 10.0 | 0.0 | 0.0 |
| (b) | 10.0 | 1.0 | 10.0 | 0.2 | 0.2 |
| (c) | 10.0 | 1.0 | 10.0 | −0.2 | −0.2 |
| (d) | 10.0 | 1.0 | 10.0 | −1.2 | 1.2 |
Matlab code
% =========================================================================================
% September 2019
% Algebraic Model Algorithm
% Nicolò Vaiana, Post-Doctoral Researcher, PhD
% Department of Structures for Engineering and Architecture
% University of Naples Federico II
% via Claudio, 21 - 80125, Napoli
% =========================================================================================
clc; clear all; close all;
%% APPLIED DISPLACEMENT TIME HISTORY
dt = 0.001; %time step
t = 0:dt:1.50; %time interval
a0 = 1; %applied displacement amplitude
fr = 1; %applied displacement frequency
u = a0*sin((2*pi*fr)*t(1:length(t))); %applied displacement vector
v = 2*pi*fr*a0*cos((2*pi*fr)*t(1:length(t))); %applied velocity vector
n = length(u); %applied displacement vector length
%% 1. INITIAL SETTINGS
%1.1 Set the five model parameters
ka = 10.0; %model parameter
kb = 1.0; %model parameter
alfa = 10.0; %model parameter
beta1 = 0.0; %model parameter
beta2 = 0.0; %model parameter
%1.2 Compute the internal model parameters
u0 = (1/2)*((((ka-kb)/10^-20)^(1/alfa))-1); %internal model parameter
f0 = ((ka-kb)/2)*((((1+2*u0)^(1-alfa))-1)/(1-alfa)); %internal model parameter
%1.3 Initialize the generalized force vector
f = zeros(1, n);
%% 2. CALCULATIONS AT EACH TIME STEP
for i = 2:n
%2.1 Update the history variable
uj = u(i-1)+sign(v(i))*(1+2*u0)-sign(v(i))*((((sign(v(i))*(1-alfa))/(ka-kb))*(f(i-1)-beta1*u(i-1)^3-beta2*u(i-1)^5-kb*u(i-1)-sign(v(i))*f0+(ka-kb)*(((1+2*u0)^(1-alfa))/(sign(v(i))*(1-alfa)))))^(1/(1-alfa)));
%2.2 Evaluate the generalized force at time t
if (sign(v(i))*uj)-2*u0 < sign(v(i))*u(i) || sign(v(i))*u(i) < sign(v(i))*uj
f(i) = beta1*u(i)^3+beta2*u(i)^5+kb*u(i)+(ka-kb)*((((1+2*u0+sign(v(i))*(u(i)-uj))^(1-alfa))/(sign(v(i))*(1-alfa)))-(((1+2*u0)^(1-alfa))/(sign(v(i))*(1-alfa))))+sign(v(i))*f0;
else
f(i) = beta1*u(i)^3+beta2*u(i)^5+kb*u(i)+sign(v(i))*f0;
end
end
%% PLOT
figure
plot(u, f, 'k', 'linewidth', 4)
set(gca, 'FontSize', 28)
set(gca, 'FontName', 'Times New Roman')
xlabel('generalized displacement'), ylabel('generalized force')
grid
zoom off
Transcendental models
In transcendental models, the output variable is computed by solving transcendental equations, namely equations involving trigonometric, inverse trigonometric, exponential, logarithmic, and/or hyperbolic functions.
Exponential models
Exponential model by Vaiana et al. (2018)
Model formulation
In the exponential model developed by Vaiana et al. (2018),[4] the generalized force at time [math]\displaystyle{ t }[/math], representing the output variable, is evaluated as a function of the generalized displacement as follows:
[math]\displaystyle{ f_t=-2\beta u_t+e^{\beta u_t}-e^{-\beta u_t} +k_b u_t-s_t\frac{k_a-k_b}{\alpha}\left[e^{-\alpha(u_t s_t-u_j s_t+2u_0)}-e^{-2\alpha u_0}\right]+f_0 s_t \quad \text{when} \quad u_t s_t \lt u_j s_t, }[/math]
[math]\displaystyle{ f_t=-2\beta u_t+e^{\beta u_t}-e^{-\beta u_t}+k_b u_t+f_0 s_t \quad \text{when} \quad u_t s_t\gt u_j s_t, }[/math]
where [math]\displaystyle{ k_a, k_b, \alpha }[/math] and [math]\displaystyle{ \beta }[/math] are the four model parameters to be calibrated from experimental or numerical tests, whereas [math]\displaystyle{ s_t }[/math] is the sign of the generalized velocity at time [math]\displaystyle{ t }[/math], that is, [math]\displaystyle{ s_t = \operatorname{sign}(\dot{u}_t) }[/math]. Furthermore, [math]\displaystyle{ u_0 }[/math]and [math]\displaystyle{ f_0 }[/math] are two internal model parameters evaluated as:
[math]\displaystyle{ u_0=-\frac{1}{2\alpha} \ln\left(\frac{10^{-20}}{k_a-k_b}\right), }[/math]
[math]\displaystyle{ f_0=\frac{k_a-k_b}{2\alpha}\left(1-e^{-2\alpha u_0}\right), }[/math]
whereas [math]\displaystyle{ u_j }[/math] is the internal variable:
[math]\displaystyle{ u_j=u_{t-\Delta t}+2u_0 s_t + \frac{s_t}{\alpha} \ln\left[\frac{\alpha s_t}{k_a-k_b} \left(-2\beta u_{t-\Delta t}+e^{\beta u_{t-\Delta t}}-e^{-\beta u_{t-\Delta t}}+k_b u_{t-\Delta t}+\frac{k_a-k_b}{\alpha} s_t e^{-2\alpha u_0} + f_0 s_t-f_{t-\Delta t}\right)\right]. }[/math]
Hysteresis loop shapes
Figure 2.1 shows four different hysteresis loop shapes obtained by applying a sinusoidal generalized displacement having unit amplitude and frequency and simulated by adopting the Exponential Model (EM) parameters listed in Table 2.1.File:Figure NV 2.tif
| [math]\displaystyle{ k_a }[/math] | [math]\displaystyle{ k_b }[/math] | [math]\displaystyle{ \alpha }[/math] | [math]\displaystyle{ \beta }[/math] | |
|---|---|---|---|---|
| (a) | 5.0 | 0.5 | 5.0 | 0.0 |
| (b) | 5.0 | −0.5 | 5.0 | 0.0 |
| (c) | 5.0 | 0.5 | 5.0 | 1.0 |
| (d) | 5.0 | 0.5 | 5.0 | −1.0 |
Matlab code
% =========================================================================================
% September 2019
% Exponential Model Algorithm
% Nicolò Vaiana, Post-Doctoral Researcher, PhD
% Department of Structures for Engineering and Architecture
% University of Naples Federico II
% via Claudio, 21 - 80125, Napoli
% =========================================================================================
clc; clear all; close all;
%% APPLIED DISPLACEMENT TIME HISTORY
dt = 0.001; %time step
t = 0:dt:1.50; %time interval
a0 = 1; %applied displacement amplitude
fr = 1; %applied displacement frequency
u = a0*sin((2*pi*fr)*t(1:length(t))); %applied displacement vector
v = 2*pi*fr*a0*cos((2*pi*fr)*t(1:length(t))); %applied velocity vector
n = length(u); %applied displacement vector length
%% 1. INITIAL SETTINGS
%1.1 Set the four model parameters
ka = 5.0; %model parameter
kb = 0.5; %model parameter
alfa = 5.0; %model parameter
beta = 1.0; %model parameter
%1.2 Compute the internal model parameters
u0 = -(1/(2*alfa))*log(10^-20/(ka-kb)); %internal model parameter
f0 = ((ka-kb)/(2*alfa))*(1-exp(-2*alfa*u0)); %internal model parameter
%1.3 Initialize the generalized force vector
f = zeros(1, n);
%% 2. CALCULATIONS AT EACH TIME STEP
for i = 2:n
%2.1 Update the history variable
uj = u(i-1)+2*u0*sign(v(i))+sign(v(i))*(1/alfa)*log(sign(v(i))*(alfa/(ka-kb))*(-2*beta*u(i-1)+exp(beta*u(i-1))-exp(-beta*u(i-1))+kb*u(i-1)+sign(v(i))*((ka-kb)/alfa)*exp(-2*alfa*u0)+sign(v(i))*f0-f(i-1)));
%2.2 Evaluate the generalized force at time t
if (sign(v(i))*uj)-2*u0 < sign(v(i))*u(i) || sign(v(i))*u(i) < sign(v(i))*uj
f(i) = -2*beta*u(i)+exp(beta*u(i))-exp(-beta*u(i))+kb*u(i)-sign(v(i))*((ka-kb)/alfa)*(exp(-alfa*(sign(v(i))*(u(i)-uj)+2*u0))-exp(-2*alfa*u0))+sign(v(i))*f0;
else
f(i) = -2*beta*u(i)+exp(beta*u(i))-exp(-beta*u(i))+kb*u(i)+sign(v(i))*f0;
end
end
%% PLOT
figure
plot(u ,f, 'k', 'linewidth', 4)
set(gca, 'FontSize', 28)
set(gca, 'FontName', 'Times New Roman')
xlabel('generalized displacement'), ylabel('generalized force')
grid
zoom off
References
- ↑ Vaiana, Nicolò; Spizzuoco, Mariacristina; Serino, Giorgio (June 2017). "Wire rope isolators for seismically base-isolated lightweight structures: Experimental characterization and mathematical modeling" (in en). Engineering Structures 140: 498–514. doi:10.1016/j.engstruct.2017.02.057. https://linkinghub.elsevier.com/retrieve/pii/S0141029617306272.
- ↑ Dimian, Mihai; Andrei, Petru (4 November 2013). Noise-driven phenomena in hysteretic systems. ISBN 9781461413745.
- ↑ 3.0 3.1 Vaiana, Nicolò; Sessa, Salvatore; Rosati, Luciano (January 2021). "A generalized class of uniaxial rate-independent models for simulating asymmetric mechanical hysteresis phenomena" (in en). Mechanical Systems and Signal Processing 146: 106984. doi:10.1016/j.ymssp.2020.106984. Bibcode: 2021MSSP..14606984V.
- ↑ 4.0 4.1 Vaiana, Nicolò; Sessa, Salvatore; Marmo, Francesco; Rosati, Luciano (26 April 2018). "A class of uniaxial phenomenological models for simulating hysteretic phenomena in rate-independent mechanical systems and materials". Nonlinear Dynamics 93 (3): 1647–1669. doi:10.1007/s11071-018-4282-2.
- ↑ Vaiana, Nicolò; Sessa, Salvatore; Marmo, Francesco; Rosati, Luciano (March 2019). "An accurate and computationally efficient uniaxial phenomenological model for steel and fiber reinforced elastomeric bearings". Composite Structures 211: 196–212. doi:10.1016/j.compstruct.2018.12.017.
 |
