openfoam/tutorials/mesh/moveDynamicMesh/simpleHarmonicMotion/shm.gnuplot
2011-01-05 18:25:32 +00:00

77 lines
2.3 KiB
Gnuplot

#------------------------------------------------------------------------------
# ========= |
# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
# \\ / O peration |
# \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
# \\/ M anipulation |
#-------------------------------------------------------------------------------
# License
# This file is part of OpenFOAM.
#
# OpenFOAM is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License
# along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
#
# Script
# shm.gnuplot
#
# Description
# Creates an .eps graph of OpenFOAM results vs analytical solution
# for a simple harmonic motion dynamicMesh case
#
#------------------------------------------------------------------------------
reset
set samples 2000
k = 5000.0
m = 9.6
c = 50.0
a = -0.1
omega = sqrt(k/m)
zeta = c/(2.0*m*omega)
phi = atan((sqrt(1.0 - zeta**2))/zeta)
A = a/sin(phi)
pos(A, t, omega, phi, zeta) = A*exp(-zeta*omega*t)*sin(sqrt(1-zeta**2)*omega*t + phi)
vel(A, t, omega, phi, zeta) = \
A*exp(-zeta*omega*t)*\
( \
sqrt(1-zeta**2)*omega*cos(sqrt(1-zeta**2)*omega*t + phi) \
- zeta*omega*sin(sqrt(1-zeta**2)*omega*t + phi) \
)
set xlabel "Time/[s]"
set ylabel "Position"
set ytics nomirror
set y2tics
set yrange [-0.1:0.1]
set y2range [-2:2]
set xzeroaxis
set terminal postscript eps color enhanced solid
set output "shm.eps"
plot \
"t_vs_cm" u 1:($2 - 0.6) w l t "Simulation, centre of mass relative to start", \
pos(A, x, omega, phi, zeta) w l t "Analytical solution, centre of mass", \
"t_vs_lv" u 1:2 w l axes x1y2 t "Simulation, vertical velocity", \
vel(A, x, omega, phi, zeta) w l axes x1y2 t "Analytical solution, vertical velocity"
#------------------------------------------------------------------------------