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:

$$\begin{equation} m\ddot{y}+c\dot{y}+ky=f(t) \label{eq:momentumbalance} \end{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:

gmsh -3 -format msh2 UnitSquare.geo

• 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:

from scipy.io import mmread, mmwrite
import scipy.sparse
import numpy as np
import sys

• Read in the arguments from the terminal call:

AFileName = sys.argv[1]
BFileName = sys.argv[2]
arg_sigma = complex(sys.argv[3])
arg_k = int(sys.argv[4])
arg_tol = sys.argv[5]

• Import the system matrices from .mtx format:

A = mmread(AFileName)
B = mmread(BFileName)

• Solve the generalized eigenvalue problem with the scipy.linalg.sparse package for sparse linear algebra:

eigenvalues, eigenvectors = scipy.sparse.linalg.eigs(A,k=arg_k,M=B,sigma=arg_sigma,tol=arg_tol)

• Reshape the arrays holding eigenvalues and eigenvectors since the openCFS input routine only supports vector-like structures:

eigenvalues = eigenvalues[np.newaxis].T
eigenvectors = eigenvectors.reshape((-1,1), order='F')

• Export the resulting arrays to the Matrix Market coordinate format:

eigenvalues = scipy.sparse.coo_matrix(eigenvalues)
eigenvectors = scipy.sparse.coo_matrix(eigenvectors)
mmwrite(sys.argv[6], eigenvalues)
mmwrite(sys.argv[7], eigenvectors)

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 the mat.xml file.

  <fileFormats>
    <input>
      <gmsh fileName="UnitSquare.msh"/>
    </input>
    <output>
      <hdf5/>
    </output>
    <materialData file="mat.xml" format="xml"/>
  </fileFormats>

• Define the computational domain as plane and assign the material:

  <domain geometryType='plane'  printGridInfo="no">
    <regionList>
      <region name="V_square" material="UnitMaterial" />
    </regionList>
  </domain>

• 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:

 <pdeList>
      <mechanic subType='planeStrain'>
        <regionList>
          <region name="V_square"/>
        </regionList>      

• 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}$.

         <bcsAndLoads>
          <fix name="S_fixed">
            <comp dof="x"/>
          </fix>
          <concentratedElem name="N_concentrated-element" dof="y" stiffnessValue="1.0"/>
        </bcsAndLoads>      

• Define the stored results:

        <storeResults>
          <nodeResult type="mechDisplacement">
            <allRegions/>
          </nodeResult>
        </storeResults>
      </mechanic>
    </pdeList>

• 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

 <concentratedElem name="N_concentrated-element" dof="y" stiffnessValue="1.0" dampingValue="3/5*2"/>

• 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 |