SDoF Bouc-Wen model with version checking

Integrating the differential equation of a Bouc-Wen model of hysteresis is very easy using Mathematica. In Mathematica 9, some additional checks prevent NDSolve[ ] from integrating the equation; the error message is NDSolve::smpf: Failure to project onto the discontinuity surface when computing Filippov continuation at time 0.  We can create a version-aware code to take into account this; if the Mathematica version is greater than 8, then remove DiscontinuityProcessing method.

The equation of motion of a single-degree-of-freedom Bouc-Wen system with external linear viscous damping is expressed as:

m{\rm{ }}\ddot u\left( t \right) + c{\rm{ }}\dot u\left( t \right) + F\left( t \right) = f\left( t \right) 

where, u\left( t \right) is the displacement, F(t) is the restoring force, f(t) is the excitation force, c is the damping coefficient and the overdot denotes the derivative with respect to time. The restoring force is expressed as:

F\left( t \right) = a\frac{{{F_y}}}{{{u_y}}}u\left( t \right) + (1 - a){F_y}{\rm{ }}z\left( t \right)

where, {F_y} is the yield force, {u_y} is the yield displacement, a is the ratio of post-yield to pre-yield (elastic) stiffness and z\left( t \right) a dimensionless hysteretic parameter that obeys the following non-linear differential equation with zero initial conditions: 

\dot z(t) = \frac{1}{{{u_y}}}\left[ {A - {{\left| {z(t)} \right|}^n}(\beta + {{\rm sgn}} (\dot u(t){\rm{ }}z(t)){\rm{ }}\gamma )} \right]\dot u(t) 

The following example evaluates the response of a single-degree-of-freedom Bouc-Wen model under a sinusoidal excitation, with the addition of sliders and other controls to allow the user to fiddle with the parameters:

Manipulate[
ParametricPlot[
Evaluate[{x1[t], a YP[[2]]/YP[[1]] x1[t] + (1 - a) YP[[2]] x3[t]} /.
Quiet@NDSolve[
{x1'[t] == x2[t],
x2'[t] == -1/Mass (c x2[t] + a YP[[2]]/YP[[1]] x1[t] + (1 - a) YP[[2]] x3[t] - Fmax Sin[(2 π)/T t]),
x3'[t] == x2[t]/YP[[1]] (1 - Abs[x3[t]]^n (γ Sign[x2[t] x3[t]] + (1 - γ))),
x1[0] == 0,
x2[0] == 0,
x3[0] == 0},
{x1[t], x2[t], x3[t]},
{t, 0, tTotal},
Method -> If[$VersionNumber > 8, {"DiscontinuityProcessing" -> False}, Automatic]]],
{t, 0, tTotal},
ImageSize -> {450, 450}, PlotRange -> 10, AxesLabel -> {"u", "F"}],
{{tTotal, 20, "Total time"}, 0.5, 100, Appearance -> "Labeled"},
{{Mass, 2.86, "m"}, 0.1, 10, 0.01, Appearance -> "Labeled"},
{{T, 4.0, "T"}, 0.1, 10, 0.01, Appearance -> "Labeled"},
{{Fmax, 8.0, "Fmax"}, 0.1, 10, 0.01, Appearance -> "Labeled"},
{{n, 2.0, "n"}, 0.1, 10, 0.01, Appearance -> "Labeled"},
{{c, 0.0, "c"}, 0.0, 10, 0.01, Appearance -> "Labeled"},
{{a, 0.05, "a"}, 0.0, 1, 0.01, Appearance -> "Labeled"},
{{γ, 0.5, "γ"}, 0.01, 1, 0.01, Appearance -> "Labeled"},
{{YP, {0.111, 2.86}}, {0, 0}, {10, 10}, Locator}]