1D Oscillator with External Eigensolver¶
Here, we will analyze the oscillations of a square plate by performing an eigenvalue analysis in different configurations using an external eigensolver (Python script):
- Undamped 1D Oscillator (click here for the xml file)
- Damped 1D Oscillator (click here for the xml file)
Here you can download the python script, mesh input file, final mesh and material file.
Setup¶
We consider a square mass of dimension 1mx1m, which is constrained in x direction on the left and is connected to a spring and a damper in the lower left corner:
The system constitutes a 1D oscillator in the y-direction, described by the following ordinary differential equation:
From the normalized form of \eqref{eq:momentumbalance} \begin{equation} \ddot{y} + 2\xi\omega_0 \dot{y}+\omega_0^2 y=\frac{f(t)}{m} \end{equation}
with \omega_0 = \sqrt{\frac{k}{m}} and \xi = \frac{c}{2\omega_0m}, we can find the following analytical expression for the eigenvalues: \begin{equation} \lambda_{1,2} = -\xi\omega_0 \pm \omega_0 \sqrt{\xi^2-1} \end{equation} Note that the eigenvalue for the undamped case (\xi=0) is simply \omega_0^2.
The corresponding eigenvalue problem in FE-formulation takes the form: \begin{equation} (\boldsymbol{K}+\lambda\boldsymbol{C}+\lambda^2\boldsymbol{M})\boldsymbol{x}=\boldsymbol{0} \label{eq:quadraticEVP} \end{equation} For the undamped case, \eqref{eq:quadraticEVP} simplifies to a generalized eigenvalue problem: \begin{equation} \boldsymbol{K}\boldsymbol{x}=\lambda\boldsymbol{M}\boldsymbol{x} \end{equation} The quadratic eigenvalue problem is also converted to a generalized one, similar to the transformation of second order ODEs to first order ODEs. The resulting generalized eigenvalue problem with increased system size is solved with the external eigensolver.
Meshing using Gmsh¶
We use Gmsh as a preprocessor.
ASCII *.geo
input files are used to create geometry and mesh.
In the following, the workflow for creating the input file for the square mass is illustrated:
-
Define geometrical constants and choose element type.
-
Create the square plate by first defining the corner points, connecting them with lines and then generating a closed surface.
-
Define a structured mesh; in this case 1 element for the entire domain is sufficient, since we are not interested in the mechanical deformation of the plate.
-
Define physical groups for relevant nodes, edges and the surface to define boundary conditions and the computational domain in the simulation.
-
Export the mesh.
-
The .msh file can be directly generated from the
UnitSquare.geo
input file. Note that only Gmsh version 2 is supported by openCFS, so we need to export the mesh in the legacy format by the following command: - Final mesh:
External Eigensolver¶
As already mentioned, we use a Python script to solve the generalized eigenvalue problem for our oscillator:
- Import the required libraries for linear algebra functions and the
sys
module for taking arguments from the terminal: - Read in the arguments from the terminal call:
-
Import the system matrices from
.mtx
format: -
Solve the generalized eigenvalue problem with the
scipy.linalg.sparse
package for sparse linear algebra: -
Reshape the arrays holding eigenvalues and eigenvectors since the openCFS input routine only supports vector-like structures:
-
Export the resulting arrays to the Matrix Market coordinate format:
Undamped Oscillator¶
We start by setting up the XML file for the undamped case with m=1kg and k=1 \frac{N}{m}. Thus, we expect an eigenvalue of \lambda=1 \frac{1}{s^2}.
-
Define the
.msh
file as input. The results are stored in a*.hdf5
file. The material parameters are taken from themat.xml
file. -
Define the computational domain as plane and assign the material:
-
Define the eigenvalue analysis step: we have a generalized eigenvalue problem and specify a shift point around which the solver should look for the eigenvalues. We only expect one eigenvalue, so we set the number of computed eigenvalues to 1.
<sequenceStep> <sequenceStep> <analysis> <eigenValue> <valuesAround> <shiftPoint> <Real>0.95</Real> <Imag>0</Imag> </shiftPoint> <number>1</number> </valuesAround> <eigenVectors normalization="norm" side="right"/> <problemType> <Generalized> <aMatrix>stiffness</aMatrix> <bMatrix>mass</bMatrix> </Generalized> </problemType> </eigenValue> </analysis>
-
Define plane strain solid mechanics as the governing partial differential equation (PDE) on our computational domain:
-
Constrain the plate in x direction on the left side. Define a spring element acting in y-direction at the lower left node with the stiffness of k=1 \frac{N}{m}.
-
Define the stored results:
-
Choose the external eigensolver for our problem and indicate what parameters should be included in the terminal call:
<linearSystems> <system> <solutionStrategy> <standard> <eigenSolver id="ext"/> </standard> </solutionStrategy> <eigenSolverList> <external id="ext"> <logging>yes</logging> <cmd>python3 EigenSolver.py</cmd> <arguments> <AFileName/> <BFileName/> <shiftPoint/> <number/> <tolerance>1e-9</tolerance> <EigenValuesFileName/> <EigenVectorsFileName/> </arguments> </external> </eigenSolverList> </system> </linearSystems> </sequenceStep> </cfsSimulation>
-
The material data is found in the
mat.xml
file:<?xml version='1.0' encoding='utf-8'?> <cfsMaterialDataBase xmlns="http://www.cfs++.org/material"> <material name="UnitMaterial"> <mechanical> <density> <linear> <real>1.0</real> </linear> </density> <elasticity> <linear> <isotropic> <elasticityModulus> <real>1.0e+9</real> </elasticityModulus> <poissonNumber> <real>0.25</real> </poissonNumber> </isotropic> </linear> </elasticity> </mechanical> </material> </cfsMaterialDataBase>
-
Running the simulation by the command
cfs UndampedOscillator
, gives us the following results:++ Reading mesh ... "UnitSquare.msh" -> regular ++ Creating PDE 'mechanic' for analysis 'eigenValue' ++ External Eigensolver - Executing the command 'python3 EigenSolver.py UndampedOscillator_AMatrix.mtx UndampedOscillator_BMatrix.mtx 0.950000+0.000000j 1 1.000000e-09 UndampedOscillator_EigenValues.mtx UndampedOscillator_EigenVectors.mtx' Import Olivia Import Olivia Mode | Real Part | Error Bounds | Imaginary Part | 1 | 1 | 1e-09 | 2.23709e-18 |
The call to the Python script is printed with all arguments and can be used for debugging.
Damped Oscillator¶
We set up the XML file for the damped case with the same mass and stiffness as before (m=1kg and k=1 \frac{N}{m}) and introduce a damping c=\frac{6}{5}\frac{Ns}{m} . This gives us \lambda_{1,2}= -0.6 \pm 0.8j \frac{1}{s} for the analytical eigenvalues. Only few changes are necessary to achieve this setup starting from the undamped case:
-
We define a quadratic eigenvalue problem and compute 2 instead of one eigenvalues around a chosen shift point:
<analysis> <eigenValue> <valuesAround> <shiftPoint> <Real>-0.6</Real> <Imag>0</Imag> </shiftPoint> <number>2</number> </valuesAround> <eigenVectors normalization="norm" side="right"/> <problemType> <Quadratic> <quadratic>mass</quadratic> <linear>damping</linear> <constant>stiffness</constant> </Quadratic> </problemType> </eigenValue> </analysis>
-
In addition to the stiffness, we define a damping at the lower left node
-
Define the solvers for the quadratic eigenvalue problem: the system is converted to a generalized eigenvalue problem (
quadratic
eigensolver), which is then solved with the external eigensolver:<linearSystems> <system> <solutionStrategy> <standard> <eigenSolver id="myquadratic"/> </standard> </solutionStrategy> <eigenSolverList> <quadratic id="myquadratic"> <generalisedEigenSolver id="ext"/> <linearisation>firstCompanion</linearisation> </quadratic> <external id="ext"> <logging>yes</logging> <cmd>python3 EigenSolver.py</cmd> <arguments> <AFileName/> <BFileName/> <shiftPoint/> <number/> <tolerance>1e-9</tolerance> <EigenValuesFileName/> <EigenVectorsFileName/> </arguments> </external> </eigenSolverList> </system> </linearSystems>
-
Running the simulation by the command
cfs DampedOscillator
, gives us the following results:++ Reading mesh ... "UnitSquare.msh" -> regular ++ Creating PDE 'mechanic' for analysis 'eigenValue' Time needed for the linearisation: 0.00065 sec ++ External Eigensolver - Executing the command 'python3 EigenSolver.py DampedOscillator_AMatrix.mtx DampedOscillator_BMatrix.mtx -0.600000+0.000000j 2 1.000000e-09 DampedOscillator_EigenValues.mtx DampedOscillator_EigenVectors.mtx' Import Olivia Import Olivia Mode | Real Part | Error Bounds | Imaginary Part | 1 | -0.6 | 1e-09 | -0.8 | 2 | -0.6 | 1e-09 | 0.8 |