动画质量弹簧阻尼器 TikZ

动画质量弹簧阻尼器 TikZ

我创建了一个质量弹簧阻尼器系统的 TikZ,该系统由以下 ODE 描述

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

其中u是输入,y是输出。我们假设w0=K=1zeta=0.1以及初始条件x1(0)=x2(0)=0

TikZ 是

\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}

现在,我该如何为这段时间u=sin(a*t)的不同时间a点的 st 制作动画ty当然,应该遵循 ODE 的解。车轮u应该沿着地面正弦曲线“行驶”。

答案1

更新:正如 OP 和评论中所要求的那样,车轮在道路上“行驶”并移动观察者。

(点击图片运行交互式 SVG [4.2 MB]。

控制 ODE 系统已使用\pstODEsolve包解决pst-ode

为了让车轮“行驶”在路面上X)以恒定速度,第三个微分方程,dX/天,需要解决。它描述了X车轮与道路接触点的坐标。微分方程组现在为

dX/天= ω r轮子/ sqrt(1+('(X))2

d质量/d=大量的

d质量/d= −ω 0 2质量偏移量) − 2ζω 0 (质量车轮

积分参数表示时间,ω wheel表示车轮的角速度,r车轮的半径。(因此,ω车轮 r车轮 是沿弯曲道路轮廓行进的距离。)ω 0是无阻尼固有频率,ζ 是阻尼比。基于 dX/天X)、车轮半径r. 轮, 地方 道路 高程 ;X)、斜率'(X)和二阶导数''(X)轮轴坐标X车轮车轮及其垂直速度车轮经过计算。车轮车轮现在作为质量弹簧阻尼器的输入。

道路概况X)被建模为波形突发:

X) = 余弦(XX抵消λXX偏移2

请注意车轮撞到“坑洼”时不太和谐的轨迹以及由此产生的响应曲线。

动画图的移动窗口是通过动态设置环境的xminxmax选项来实现的axis

ODE 求解器(RKF45pkg 中的算法)pst-ode是用 PostScript 编码的。最近,Marcel Krüger 实现了一个PostScript 解释器在 Lua 中,ps2pdf不再需要这样做。

运行lualatex3 次即可排版为 PDF。对于 GIF 和 SVG 输出,请取消注释源代码开头的其中一行,然后按照其中的说明进行操作。

%\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}

相关内容