linear_elasticity/modal_analysis_declarative.py

Description

Modal analysis of a wheel set.

The first six modes are the rigid body modes because no boundary conditions are applied.

Running the simulation:

sfepy-run sfepy/examples/linear_elasticity/modal_analysis_declarative.py

The eigenvalues are saved to wheelset_eigs.txt and the eigenvectros to wheelset.vtk. View the results using:

sfepy-view wheelset.vtk -f u003:wu003:f30%:p0 1:vw:p0

The first six frequencies calculated by SfePy:

[11.272, 11.322, 34.432, 80.711, 80.895, 93.149]

The results of modal analysis performed in Ansys:

[11.306, 11.316, 34.486, 80.901, 81.139, 93.472]

source code

"""
Modal analysis of a wheel set.

The first six modes are the rigid body modes because no boundary
conditions are applied.

Running the simulation::

  sfepy-run sfepy/examples/linear_elasticity/modal_analysis_declarative.py

The eigenvalues are saved to wheelset_eigs.txt and the eigenvectros to
wheelset.vtk. View the results using::

  sfepy-view wheelset.vtk -f u003:wu003:f30%:p0 1:vw:p0

The first six frequencies calculated by SfePy::

  [11.272, 11.322, 34.432, 80.711, 80.895, 93.149]

The results of modal analysis performed in Ansys::

  [11.306, 11.316, 34.486, 80.901, 81.139, 93.472]
"""
import numpy as nm
from sfepy.base.base import output
from sfepy.mechanics.matcoefs import stiffness_from_youngpoisson
from sfepy import data_dir


def report_eigs(pb, evp):
    eigs = evp.eigs
    n_rbm = evp.pb.conf.n_rbm

    output('%d eigenvalues converged (%d ignored as rigid body modes)' %
           (len(eigs), n_rbm))

    eigs = eigs[n_rbm:]

    omegas = nm.sqrt(eigs)
    freqs = omegas / (2 * nm.pi)

    output('number |         eigenvalue |  angular frequency '
           '|          frequency')
    for ii, eig in enumerate(eigs):
        output('%6d | %17.12e | %17.12e | %17.12e'
               % (ii + 1, eig, omegas[ii], freqs[ii]))


def define(n_eigs=6, approx_order=1, density=7850., young=210e9, poisson=0.3):
    filename_mesh, dim = data_dir + '/meshes/3d/wheelset.vtk', 3

    n_rbm = 0

    options = {
        'n_eigs': n_eigs + n_rbm,
        'eigs_only': False,
        'post_process_hook_final': 'report_eigs',
        'evps': 'eig',
    }

    regions = {
        'Omega': 'all',
        'Fix': ('vertices in (x < -1.08999)', 'vertex'),
    }

    materials = {
        'm': ({
            'D': stiffness_from_youngpoisson(dim, young, poisson),
            'rho': density,
        },),
    }

    fields = {
        'displacement': ('real', 'vector', 'Omega', approx_order),
    }

    variables = {
        'u': ('unknown field', 'displacement'),
        'v': ('test field', 'displacement', 'u'),
    }

    integrals = {
        'i': 2 * approx_order,
    }

    equations = {
        'lhs': 'dw_lin_elastic.i.Omega(m.D, v, u)',
        'rhs': 'dw_dot.i.Omega(m.rho, v, u)',
    }

    ebcs = {
        'fix': ('Fix', {'u.all': 0.0})  # fix rigid body modes
    }

    solvers = {
        # 'eig': ('eig.matlab', {
        #     'method': 'eigs',
        #     'which': 'sm',
        #     'eps': 1e-6,
        # }),
        #'eig': ('eig.primme', {
        #    'which': 'SM',
        #    'tol': 1e-8,
        #}),
        'eig': ('eig.scipy', {
            'method': 'eigsh',
            'which': 'LM',
            'sigma': 0,
            'tol': 1e-8,
            'linear_solver': ('ls.scipy_superlu', {}),
            # 'linear_solver': ('ls.cholesky', {}),
            # 'linear_solver': ('ls.mumps', {}),
            # 'ls.cholesky' and 'ls.mumps' linear solvers are much faster,
            # but scikit-sparse package and MUMPS library are required
            'maxiter': 1000,
        }),
    }

    return locals()