Animate mass-spring-damper TikZ

animationscalculationstikz-pgf

I've create a TikZ of a mass-spring-damper system that is described by the following ODEs

dx1/dt = x2
dx2/dt = -w0^2*x1 -2*zeta*w0*x2 + u
y = K*w0*x1 + 2*zeta*w0*x2

wherein u is the input and y is the output. Let's assumed w0=K=1 and zeta=0.1 and the initial conditions x1(0)=x2(0)=0.

The TikZ is

\documentclass{article}
\usepackage{tikz}
\usetikzlibrary{arrows, shapes, patterns, calc, decorations.pathmorphing, decorations.markings, positioning, animations}

\begin{document}
    \begin{tikzpicture}
        \tikzstyle{wheel} = [draw, circle, minimum size = 0.75cm]
        \tikzstyle{mass} = [draw, rectangle, minimum height = 1cm, minimum width = 2cm]
        \tikzstyle{spring} = [decorate, decoration = {zigzag, pre length = 0.3cm, post length = 0.3cm, segment length = 6}]
        \tikzstyle{damper} = [decoration = {markings, mark connection node = dmp, mark = at position 0.5 with 
        {
            \node (dmp) [inner sep = 0pt, transform shape, rotate = -90, minimum width = 5pt, minimum height = 10pt, draw=none] {};
            \draw ($(dmp.north east)+(1.5pt,0)$) -- (dmp.south east) -- (dmp.south west) -- ($(dmp.north west)+(1.5pt,0)$);
            \draw ($(dmp.north)+(0,-1.5pt)$) -- ($(dmp.north)+(0,1.5pt)$);
        }
        }, decorate]
        \tikzstyle{groundflat} = [fill, pattern = north east lines, draw = none, minimum width = 0.75cm, minimum height = 0.3cm]
        \draw[fill, pattern = north east lines] (-.5,-.25) cos (0,0) sin (.5,.25) coordinate (top) cos (1,0) sin (1.5,-.25) cos (2,0) sin (2.5,.25) cos (3,0) sin (3.5,-.25) |- (-.5,-1)--cycle;
        \node[wheel, above = 0pt of top] (u) {$u$};
        \node[mass, above = 2cm of u] (m) {$y$};
        \draw [-] (u.north) |- ++(0.5,0.25cm)coordinate (uright);
        \draw [-] (u.north) |- ++(-0.5,0.25cm)coordinate (uleft);
        \draw [spring] (uleft) -- ($(m.south) +(-0.5,0)$);
        \draw [damper] (uright) -- ($(m.south) +(0.5,0)$);
    \end{tikzpicture}
\end{document}

Now, how can I animate this s.t. u=sin(a*t) with different a along a certain time t. y should, of course, follow the solution of the ODEs. And the u wheel should "ride" along the ground sinusoid.

Best Answer

Update: Wheel "riding" on the road and moving observer, as requested in OP and comments.

(Click on the image to run the interactive SVG [4.2 MiB] in the browser.)

The governing ODE system has been solved with \pstODEsolve from package pst-ode.

For the wheel to "ride" on the road profile u(x) at constant speed, a third differential equation, dx/dt, needs to be solved. It describes the movement of the x coordinate of the wheel's contact point with the road. The system of differential equations is now

dx/dt = ωwheel rwheel / sqrt(1+(u'(x))2)

dymass/dt = vmass

dvmass/dt = −ω02(ymassyoffsetywheel) − 2ζω0(vmassvwheel)

The integration parameter t represents time, ωwheel the wheel's angular speed and rwheel its radius. (Hence, ωwheel rwheel t is the distance covered along the curvy road profile.) ω0 is the undamped natural frequency and ζ is the damping ratio. Based on dx/dt, x(t), wheel radius rwheel, local road elevation u(x), slope u'(x) and 2nd derivative u''(x) the wheel axle coordinates xwheel, ywheel and its vert. velocity vwheel are calculated. ywheel and vwheel now serve as input for the mass-spring-damper system.

The road profile u(x) is modelled as a waveform burst:

u(x) = cos(xxoffset) eλ(xxoffset)2

Notice the not-so-harmonic trajectory of the wheel where it bumps into the "potholes" as well as the resulting response curve.

The moving window of the animated plot is achieved by dynamically setting the xmin and xmax options of the axis environment.

The ODE solver (RKF45 method) in pkg pst-ode was coded in PostScript. Recently, Marcel Krüger implemented a PostScript interpreter in Lua such that ps2pdf is not required anymore.

Typeset as PDF by running lualatex 3 times. For GIF and SVG output, uncomment one of the lines at the beginning of the source and follow the instructions therein.

%\PassOptionsToPackage{export}{animate} % lualatex <file> ; convert -density 300 -alpha remove -delay 4 <file>.pdf <file>.gif
%\PassOptionsToPackage{dvisvgm,autoplay}{animate} % dvilualatex <file> ; dvisvgm --zoom=-1 --font-format=woff2 --bbox=papersize <file>.dvi
\documentclass[margin=3pt]{standalone}

\usepackage{pst-ode}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\usepgfplotslibrary{fillbetween}
\usetikzlibrary{arrows.meta, calc, decorations}
\usepackage{listofitems} % read space separated items
\usepackage{animate}
\usepackage{xsavebox} %\xsbox

% adjustable parameters & definitions
\def\tEnd{80} % max. time (integration parameter)
\def\animationframes{201} % number of output points for animated objects
\def\outputpoints{1001} % number of output points for curves
\def\omegaWheel{1.0} % wheel angular velocity
\def\rWheel{1.0} % wheel radius
\def\rCenterOfMass{0.16em} % center-of-mass symbol size
\def\springlength{6} % length of spring in neutral state (and of damper cylinder)
\def\springturns{10} % (integer!) number of spring windings
\def\windowsizeright{5} % moving window x-size around moving object
\def\windowsizeleft{25}
\def\ymin{-1} % y axis limits
\def\ymax{17}
\pstVerb{
  tx@Dict begin
  % mass-spring-damper properties
  /w0 1.0 def
  /zeta 0.1 def
  % initial conditions
  /x0     0.0 def % x(t0)=0 starting x coordinate
  /yMass0 0.0 def % y coord of mass centre (/wo vert. offset)
  /vMass0 0.0 def % y velocity of mass centre
  % road profile U(x) with wave-form burst, and its first and second derivatives U'(x), U''(x)
  /A 1 def        % amplitude
  /a 1 def        % mode
  /lambda (1/40) AlgParser cvx exec def % burst growth/decay rate
  /xOff 20.0 def  % burst offset
  /U (A*Euler^(-lambda*(x[0]-xOff)^2)*cos(a*(x[0]-xOff))) AlgParser cvx bind def
  /dUdx (-A*Euler^(-lambda*(x[0]-xOff)^2)*a*sin(a*(x[0]-xOff))-2*lambda*(x[0]-xOff)*U) AlgParser cvx bind def
  /d2Udx2 (-(a^2+2*lambda+4*lambda^2*(x[0]-xOff)^2)*U-4*lambda*(x[0]-xOff)*dUdx) AlgParser cvx bind def
  % do not change anything below
  /rWheel \rWheel\space def
  /omegaWheel \omegaWheel\space def
  /massWheelOff (rWheel+1.5*\springlength+2)  AlgParser cvx exec def
  /dxdt (rWheel*omegaWheel/sqrt(1+dUdx^2)) AlgParser cvx bind def
  % local wheel hub coordinates
  /xWheel (x[0]-rWheel*dUdx/sqrt(1+dUdx^2)) AlgParser cvx bind def
  /yWheel (U+rWheel/sqrt(1+dUdx^2)) AlgParser cvx bind def
  /dyWheeldx (dUdx*(1-d2Udx2/(1+dUdx^2)*(yWheel-U))) AlgParser cvx bind def
  % angular position of wheel in terms of t
  /thetaWheel {omegaWheel t mul 2 Pi mul div dup cvi sub neg 360 mul} bind def
  end
}

% solve equations of motion with many output points for smooth curves
\pstODEsolve[algebraicAll,saveData]{curves}{% PS variable and file that take result table
                           % table format to be saved in `curves' and `curves.dat':
  x[0] |                   %    x coordinate of contact point
  U |                      %    y coordinate U(x) of contact point
  xWheel | yWheel |        %    wheel hub x (same as that of mass) and y coordinates
  x[1]+massWheelOff+rWheel %    y coordinate of mass (response with const. vertical offset)
}{0}{\tEnd}{\outputpoints}{% t_0, t_end, number of t steps plus 1
  x0 | yMass0 | vMass0     % initial conditions
}{                         % ODE system's RHS:
  dxdt |                   %   dx/dt, movement of contact point (x coordinate)
  x[2] |                   %   centre of mass vert. velocity
  -w0^2*(x[1]-yWheel+rWheel)%  centre of mass vert. acceleration
    -2*zeta*w0*(x[2]-dyWheeldx*dxdt)
}

% solve equations with fewer output points for mass-spring-damper animation; save solution
% table (wheel and mass centre coordinates, valve angular position) in `objects' and `objects.dat'
\pstODEsolve[algebraicAll,saveData]{objects}{
  xWheel | yWheel | x[1]+massWheelOff+rWheel | thetaWheel-90
}{0}{\tEnd}{\animationframes}{x0 | yMass0 | vMass0}{
  dxdt | x[2] | -w0^2*(x[1]-yWheel+rWheel)-2*zeta*w0*(x[2]-dyWheeldx*dxdt)
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \fileopenr{<file stream>}{<file name>}, opens file for reading
\newcommand\fileopenr[2]{%
  \newread#1%
  \immediate\openin#1=#2%
}
% \readtolist[<sep char>]{<file stream>}{\list}
% reads a line from file stream and splits at <sep char> into \list[1], \list[2], ...
\newcommand\readtolist[3][,]{{%
  \setsepchar{#1}%
  \immediate\read#2 to \inputline%
    \ifeof#2
      \immediate\closein#2%
      \ifdefined\multiframebreak\multiframebreak\fi%
    \else%
      \greadlist*#3\inputline%
    \fi%
}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\tikzset{spring/.style 2 args= {decorate, decoration = {zigzag, segment length = \dimexpr(#1-1pt)/\springturns\relax, amplitude=#2}}}%

\begin{document}
\IfFileExists{curves.dat}{}{dummy text\end{document}}%
\xsbox{centreofmass}{%
  \tikz[radius=\rCenterOfMass] {% from https://tex.stackexchange.com/questions/23453
    \fill (0,0) -- ++(\rCenterOfMass,0) arc [start angle=0,end angle=90] -- ++(0,{-2*\rCenterOfMass}) arc [start angle=270, end angle=180];%
    \draw [thin] (0,0) circle;%
  }%
}%
\begin{animateinline}[controls]{25}
  \fileopenr{\data}{objects.dat}
  \readtolist[ ]{\data}{\table}
  \multiframe{10000}{}{
    \begin{tikzpicture}[line cap=rect, line join=bevel]
      \begin{axis}[
        unit vector ratio = 1, axis on top,
        xmin={\table[1]-\windowsizeleft}, xmax={\table[1]+\windowsizeright},
        ymin=\ymin, ymax=\ymax,
        xtick distance=5,
        ytick distance=2,
        xlabel={$x$}, ylabel={$y$}, ylabel style={rotate=-90},
        extra x tick labels={\phantom{\strut00}},extra x ticks={\pgfkeysvalueof{/pgfplots/xmax}}, % prevents wobbling of axis
      ]
        % road profile; off-window coordinates filtered away for smaller file size
        \addplot [name path=road,no markers,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\pgfkeysvalueof{/pgfplots/xmax}+0.1 ? nan : x}]
          table [x index=0, y index=1] {curves.dat} -- (\pgfkeysvalueof{/pgfplots/xmax},0);
        \path [name path=bottom] (0,\pgfkeysvalueof{/pgfplots/ymin}) -- (\pgfkeysvalueof{/pgfplots/xmax},\pgfkeysvalueof{/pgfplots/ymin});
        \addplot [fill=black!20!white] fill between[of=road and bottom];
        % wheel hub trajectory
        \addplot [no markers,blue,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\table[1] ? nan : x}]
          table [x index=2, y index=3] {curves.dat};
        % centre of mass trajectory
        \addplot [no markers,red,unbounded coords=jump,x filter/.expression={x<\pgfkeysvalueof{/pgfplots/xmin}-0.1 || x>\table[1] ? nan : x}]
          table [x index=2, y index=4] {curves.dat};
        \coordinate (w) at (\table[1],\table[2]);  % wheel hub
        \coordinate (m) at (\table[1],\table[3]);  % centre of mass
        % wheel with valve
        \fill[even odd rule] (w) circle [radius={0.7*\rWheel}] (w) circle [radius=\rWheel];
        \draw[thick] ($(w)+(axis direction cs:{0.6*\rWheel*cos(\table[4])},{0.6*\rWheel*sin(\table[4])})$)
          -- ++(axis direction cs:{0.1*\rWheel*cos(\table[4])},{0.1*\rWheel*sin(\table[4])});
        % mass
        \node at (m) {\thecentreofmass};
        \draw (m) ++(axis direction cs:-1,-1) rectangle ++(axis direction cs:2,2);
        % spring/damper support
        \draw (m) ++(axis direction cs:-1,-1) -- ++(axis direction cs:0,{-0.25-0.25*\springlength}) coordinate (sprtop);
        \path (m) ++(axis direction cs:1,-1) coordinate (dmptop);
        \draw (w) ++(axis direction cs:0,{\rWheel+0.5}) -| ++(axis direction cs:-1,{0.25+0.25*\springlength}) coordinate (sprbot);
        \draw let \p1=(axis direction cs:0.35,0) in [{Circle[open, length=\x1]}-, shorten <=-0.5*\x1] (w) -- ++(axis direction cs:0,{\rWheel+0.5}) -| ++(axis direction cs:1,0.5) coordinate (dmpbot);
        % spring
        \path let \p1=($(sprtop)-(sprbot)$), \p2=(axis direction cs:0.6,0) % length and amplitude (half width)
              in [draw,spring={\y1}{\x2}] (sprbot) -- (sprtop);
        % damper lower part
        \draw (dmpbot) -| +(axis direction cs:0.4,\springlength) (dmpbot) -| +(axis direction cs:-0.4,\springlength);
        % damper upper part
        \draw (dmptop) -- ++(axis direction cs:0,-\springlength) ++(axis direction cs:-0.3,0) -- ++(axis direction cs:0.6,0);
      \end{axis}
    \end{tikzpicture}%
    \readtolist[ ]{\data}{\table}
  }
\end{animateinline}
\end{document}