20.1. V2.08.011: FORMA11 - Modal dynamic analysis

The tests in this tutorial provide further practice with modal dynamics in Code_Aster. Three geometries, made entirely of a linear elastic material, are explored. We perform several modal analyses which involves the search for eigen frequencies and associated modes.

The models are:

  • Model A: Search of the multiple modes of a beam (POU_D_E elements)

  • Model B: Search of the rigid body modes of a sphere (3D elements)

  • Model C: Modal analysis of a cooling tower (DKT elements)

The aims of this tutorial are to:

  • Compute the multiple modes of a cantilevered beam

  • Compute the rigid body modes on an unconstrained sphere

  • Use of the macro-command CALC_MODES

20.1.1. Model A : Modal analysis of a cantilever beam

20.1.1.1. Problem description

20.1.1.1.1. Objective

The objective of this model is to determine the first 11 eigen frequencies of the cantilever beam using four methods:

  • SORENSEN method

  • TRI_DIAG (Lanczos) method

  • JACOBI (Bathe and Wilson) method

  • Inverse iteration method

20.1.1.1.2. Geometry

../../_images/fig1_v2.08.011.png

Fig. 20.1.1 Schematic of cantilever beam

We consider a beam of rectangular cross-section (0.05 m \(\times\) 0.05 m) and length 0.5 m fixed at point A. We will analyze the free vibration of the beam.

20.1.1.1.3. Material properties

The material is linear elastic:

  • Young’s modulus: \(E\) = 210 GPa

  • Poisson’s ratio: \(\nu\) = 0.3

  • Mass density: \(\rho\) = 7000 kg/m \(^3\)

20.1.1.1.4. Boundary conditions and loads

The beam is fixed at point A

20.1.1.2. Creating and running the tutorial

20.1.1.2.1. Geometry and mesh with Salome-Meca

The line mesh can be constructed interactively using Salomé-Meca.

  • Define the points A, B and then the line AB. The size of the elements is the constant. The mesh consists of 10 SEG2 elements and 11 nodes

  • In the Mesh module of Salomé define a node group (point_A) for the vertex A and an element group (bar) for the line AB.

  • Save the mesh in MED format.

Straight Euler beam elements (POU_D_E) will be used for the model.

You can also create the model and the mesh using a Python script (download from here). The annotated script is listed below.

import sys
import salome
import GEOM
import SMESH
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math

ExportPATH="/home/banerjee/Salome/forma11/"

salome.salome_init()

#-----------
# Geometry
#-----------
geompy = geomBuilder.New()

# Create points and line
point_A = geompy.MakeVertex(0, 0, 0)
point_B = geompy.MakeVertex(0.5, 0, 0)
bar = geompy.MakeLineTwoPnt(point_A, point_B)

# Create and assign groups
g_point_A = geompy.CreateGroup(bar, geompy.ShapeType["VERTEX"])
g_bar = geompy.CreateGroup(bar, geompy.ShapeType["EDGE"])
geompy.UnionIDs(g_point_A, [2])
geompy.UnionIDs(g_bar, [1])

# Add the geometries to the study
geompy.addToStudy(point_A, 'point_A')
geompy.addToStudy(point_B, 'point_B')
geompy.addToStudy(bar, 'bar')

# Add the groups to the study
geompy.addToStudyInFather(bar, g_point_A, 'point_A')
geompy.addToStudyInFather(bar, g_bar, 'bar')

#-----------
# Mesh
#-----------
smesh = smeshBuilder.New()

# Define mesh algorithms
mesh = smesh.Mesh(bar)
Regular_1D = mesh.Segment()
Number_of_Segments_1 = Regular_1D.NumberOfSegments(10)
smesh.SetName(mesh.GetMesh(), 'mesh')
smesh.SetName(Regular_1D.GetAlgorithm(), 'Regular_1D')
smesh.SetName(Number_of_Segments_1, 'Number of Segments_1')

# Compute the mesh
isDone = mesh.Compute()

# Map geometry groups to mesh
bar_2 = mesh.GroupOnGeom(g_bar, 'bar', SMESH.EDGE)
point_A_2 = mesh.GroupOnGeom(g_point_A, 'point_A', SMESH.NODE)

# Export mesh
mesh.ExportMED( r''+ExportPATH+'forma11a.mmed'+'')

# Update GUI tree
if salome.sg.hasDesktop():
  salome.sg.updateObjBrowser()

20.1.1.2.2. Creation of the command file

  • Read of the MED mesh (LIRE_MAILLAGE)

  • Define the finite element type (AFFE_MODELE) and use beam elements POU_D_E (EULER_BEAM).

  • Assign the characteristics of the beam elements (AFFE_CARA_ELEM). The section of all elements of the beam is the same. Use SECTION = ‘RECTANGLE’ and CARA = (‘HY’, ‘HZ’).

  • Define and assign the material (DEFI_MATERIAU and AFFE_MATERIAU). You will need to define RHO so that the mass matrix can be computed.

  • Assign the boundary conditions (AFFE_CHAR_MECA): * Constrain the fixed end (LIAISON = ‘ENCASTRE’)

  • Compute the element stiffness and mass matrices with the macro ASSEMBLAGE. Name the stiffness matrix stif_mat, the mass matrix mass_mat, and the DOF vector num_dof. The macro will create the assembled matrices.

  • Solve the modal analysis problem with CALC_MODES [U4.52.02]

You can try the following:

  • Calculate the 11 smallest eigen frequencies and the associated displacement modes

  • Write the eigen modes (IMPR_RESU) in MED format MED for visualization

For additional exploration:

  • Calculate the eigen frequencies and the associated modes present in the band of frequencies 0 Hz and 6000 Hz, with the four methods (OPTION) SORENSEN, LANCZOS, TRI_DIAG (Bathe and Wilson), and QZ (keyword factor SOLVEUR_MODAL = _F (METHODE = …)) and then using inverse iterations (OPTION = ‘AJUSTE’).

  • Print the eigen modes in MED format

The annotated Python command file for Model A is forma11a.comm. The contents of the file are listed below for convenience.

# Modal analysis of a straight bar made of beam elements
# POU_D_E
DEBUT(LANG='EN')

# Read the mesh
mesh = LIRE_MAILLAGE(FORMAT = 'MED',
                     UNITE = 20)

# Check mesh connectivity
MACR_INFO_MAIL(MAILLAGE = mesh,
               QUALITE = 'OUI',
               INTERPENETRATION = 'OUI',
               CONNEXITE = 'OUI',
               TAILLE = 'OUI',
               PROP_CALCUL = 'OUI')

# Assign model
model = AFFE_MODELE(AFFE = _F(MODELISATION=('POU_D_E'),
                              PHENOMENE='MECANIQUE',
                              TOUT='OUI'),
                    MAILLAGE = mesh)

# Beam element properties
section = AFFE_CARA_ELEM(MODELE = model,
                         POUTRE = _F(GROUP_MA = ('bar'),
                                     SECTION = 'RECTANGLE',
                                     CARA = ('HY', 'HZ'),
                                     VALE = (0.05, 0.05)))

# Assign material
steel = DEFI_MATERIAU(ELAS=_F(E = 210000000000.0,
                              NU = 0.3,
                              RHO = 7800.0))

mater = AFFE_MATERIAU(AFFE = _F(MATER = (steel),
                                TOUT = 'OUI'),
                      MODELE = model)

# Assign BCs
disp_bc = AFFE_CHAR_MECA(DDL_IMPO = _F(GROUP_NO = ('point_A'),
                                       LIAISON = 'ENCASTRE'),
                         MODELE = model)

# Assemble stiffness and mass matrices
ASSEMBLAGE(CARA_ELEM = section,
           CHAM_MATER = mater,
           CHARGE = (disp_bc),
           MATR_ASSE = (_F(MATRICE = CO('stif_mat'),
                           OPTION = 'RIGI_MECA'),
                        _F(MATRICE = CO('mass_mat'),
                           OPTION = 'MASS_MECA')),
           MODELE = model,
           NUME_DDL = CO('num_dof'))

# Modal analysis (0): PLUS_PETITE/TRI_DIAG
mode_0 = CALC_MODES(CALC_FREQ = _F(NMAX_FREQ = 11),
                    INFO = 2,
                    MATR_MASS = mass_mat,
                    MATR_RIGI = stif_mat,
                    OPTION = 'PLUS_PETITE',
                    SOLVEUR_MODAL = _F(METHODE = 'TRI_DIAG'))

# Modal analysis (1): PLUS_PETITE/QZ
mode_1 = CALC_MODES(CALC_FREQ = _F(NMAX_FREQ = 11),
                    INFO = 2,
                    MATR_MASS = mass_mat,
                    MATR_RIGI = stif_mat,
                    OPTION = 'PLUS_PETITE',
                    SOLVEUR_MODAL = _F(METHODE = 'QZ'))

# Modal analysis (2): CENTRE/QZ
mode_2 = CALC_MODES(CALC_FREQ = _F(FREQ = 2597.04,
                                   NMAX_FREQ = 11),
                    INFO = 2,
                    MATR_MASS = mass_mat,
                    MATR_RIGI = stif_mat,
                    OPTION = 'CENTRE',
                    SOLVEUR_MODAL = _F(METHODE='QZ'))

# Modal analysis (3): ALL/QZ
mode_3 = CALC_MODES(MATR_MASS = mass_mat,
                    MATR_RIGI = stif_mat,
                    OPTION = 'TOUT',
                    SOLVEUR_MODAL = _F(METHODE = 'QZ'))

# Modal analysis (4): BANDE/SORENSEN
mode_sor = CALC_MODES(CALC_FREQ = _F(FREQ = (0.0, 6000.0)),
                      MATR_MASS = mass_mat,
                      MATR_RIGI = stif_mat,
                      NORM_MODE = _F(NORME = 'EUCL'),
                      OPTION = 'BANDE',
                      SOLVEUR_MODAL = _F(METHODE='SORENSEN'),
                      TITRE = 'Sorensen')

# Modal analysis (5): BANDE/Lanzcos (TRI_DIAG)
mode_lan = CALC_MODES(CALC_FREQ = _F(FREQ = (0.0, 6000.0)),
                      MATR_MASS = mass_mat,
                      MATR_RIGI = stif_mat,
                      NORM_MODE = _F(NORME = 'EUCL'),
                      OPTION ='BANDE',
                      SOLVEUR_MODAL = _F(METHODE = 'TRI_DIAG'),
                      TITRE ='Lanzcos')

# Modal analysis (6): BANDE/Bathe (JACOBI)
mode_jac = CALC_MODES(CALC_FREQ = _F(FREQ = (0.0, 6000.0)),
                      MATR_MASS = mass_mat,
                      MATR_RIGI = stif_mat,
                      NORM_MODE = _F(NORME = 'EUCL'),
                      OPTION = 'BANDE',
                      SOLVEUR_MODAL = _F(METHODE = 'JACOBI'),
                      TITRE = 'Bathe',
                      VERI_MODE = _F(SEUIL = 1e-05))

# Modal analysis (7): BANDE/QZ
mode_qz = CALC_MODES(CALC_FREQ = _F(FREQ = (0.0, 6000.0)),
                     MATR_MASS = mass_mat,
                     MATR_RIGI = stif_mat,
                     NORM_MODE = _F(NORME = 'EUCL'),
                     OPTION = 'BANDE',
                     SOLVEUR_MODAL =_F(METHODE = 'QZ'),
                     TITRE = 'QZ',
                     VERI_MODE = _F(SEUIL = 1e-05))

# Modal analysis (8): AJUSTE/MULT_FRONT
mode_inv = CALC_MODES(CALC_FREQ = _F(FREQ = (0.0, 6000.0)),
                      MATR_MASS = mass_mat,
                      MATR_RIGI = stif_mat,
                      NORM_MODE = _F(NORME = 'EUCL'),
                      OPTION = 'AJUSTE',
                      SOLVEUR = _F(METHODE = 'MULT_FRONT'),
                      TITRE = 'inverse')

# Write out material and section properties
mat_sec = POST_ELEM(MODELE = model,
                    CHAM_MATER = mater,
                    CARA_ELEM = section,
                    MASS_INER = _F(TOUT = 'OUI'))
IMPR_TABLE(TABLE = mat_sec,
           SEPARATEUR = ',')

# Write results in tables
IMPR_RESU(FORMAT = 'RESULTAT',
          RESU = _F(RESULTAT = mode_0,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT = 'RESULTAT',
          RESU = _F(RESULTAT = mode_1,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_2,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_3,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_sor,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_lan,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_jac,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_qz,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

IMPR_RESU(FORMAT='RESULTAT',
          RESU = _F(RESULTAT = mode_inv,
                    TOUT_PARA = 'OUI',
                    TOUT_CHAM = 'NON'),
          UNITE = 80)

# Write displacement results in tables
IMPR_RESU(FORMAT = 'RESULTAT',
          RESU = (_F(IMPR_COOR = 'OUI',
                     NOM_CHAM = 'DEPL',
                     NOM_CMP = ('DY', 'DZ'),
                     NUME_ORDRE = (1, 2),
                     RESULTAT = mode_sor),
                  _F(IMPR_COOR='OUI',
                     NOM_CHAM = 'DEPL',
                     NOM_CMP = ('DY', 'DZ'),
                     NUME_ORDRE = (1, 2),
                     RESULTAT = mode_lan),
                  _F(IMPR_COOR='OUI',
                     NOM_CHAM = 'DEPL',
                     NOM_CMP = ('DY', 'DZ'),
                     NUME_ORDRE = (1, 2),
                     RESULTAT = mode_jac),
                  _F(IMPR_COOR='OUI',
                     NOM_CHAM = 'DEPL',
                     NOM_CMP = ('DY', 'DZ'),
                     NUME_ORDRE = (1, 2),
                     RESULTAT = mode_qz),
                  _F(IMPR_COOR='OUI',
                     NOM_CHAM = 'DEPL',
                     NOM_CMP = ('DY', 'DZ'),
                     NUME_ORDRE = (1, 2),
                     RESULTAT = mode_inv)),
          UNITE = 81)

# Write out results in MED format for viz
IMPR_RESU(FORMAT = 'MED',
          RESU = (_F(MAILLAGE = mesh),
                  _F(RESULTAT = mode_0),
                  _F(RESULTAT = mode_1),
                  _F(RESULTAT = mode_2),
                  _F(RESULTAT = mode_3),
                  _F(RESULTAT = mode_sor),
                  _F(RESULTAT = mode_lan),
                  _F(RESULTAT = mode_jac),
                  _F(RESULTAT = mode_qz),
                  _F(RESULTAT = mode_inv)),
          UNITE = 82)

# Verification tests
TEST_RESU(RESU = _F(NUME_ORDRE = 6,
                    PARA = 'FREQ',
                    RESULTAT = mode_1,
                    VALE_CALC = 2597.040657041,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 11,
                    PARA = 'FREQ',
                    RESULTAT = mode_2,
                    VALE_CALC = 5769.9054933878997,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 1,
                    PARA = 'FREQ',
                    RESULTAT = mode_3,
                    VALE_CALC = 167.63819422541999,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 1,
                    PARA = 'FREQ',
                    RESULTAT = mode_sor,
                    VALE_CALC = 167.63819422438999,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 4,
                    PARA = 'FREQ',
                    RESULTAT = mode_lan,
                    VALE_CALC = 1050.6045040198001,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 7,
                    PARA = 'FREQ',
                    RESULTAT = mode_jac,
                    VALE_CALC = 2942.3746290629001,
                    CRITERE = 'RELATIF'))

TEST_RESU(RESU = _F(NUME_ORDRE = 11,
                    PARA = 'FREQ',
                    RESULTAT = mode_qz,
                    VALE_CALC = 5769.9054933878997,
                    CRITERE = 'RELATIF'))

# End
FIN()

20.1.1.2.3. Postprocessing with Salome-Meca

To visualize the modal deformations with Salomé-Meca, it is convenient to use the Results tab in AsterStudy. Some of the modes from the SORENSEN algorithm are shown below.

../../_images/forma11a_mode0_167.png

Fig. 20.1.2 Mode 1 - 167 Hz

../../_images/forma11a_mode0_1050.png

Fig. 20.1.3 Mode 3 - 1050 Hz

../../_images/forma11a_mode0_2940.png

Fig. 20.1.4 Mode 7 - 2942 Hz

../../_images/forma11a_mode0_5769.png

Fig. 20.1.5 Mode 10 - 5770 Hz

20.1.1.3. Remarks

  • A good practice consists in using the option BANDE rather than PLUS_PETITE or CENTRE.

  • The number of frequencies is determined automatically by the Sturm test and thus all multiplicities are taken into account.

  • The method of SORENSEN is the fastest (in spite of its iterative character which makes it possible to obtain a better precision).

  • We can project on a smaller sized space (compared to TRI_DIAG and JACOBI) and we do not need a factorization (compared to the algorithms of the inverse powers activated when one uses OPTION = ‘AJUSTE’/’SEPARE’/’PROCHE’ (OPTION = ‘ADJUSTED’ / ‘SEPARATED’ / ‘NEAR’).

  • In the presence of multiple modes, small modifications of the parameters of the modal method (or of the associated mesh DOFs NUME_DDL) produces different eigenvectors. But they always generate the same eigen space. You can convince yourself by visualizing the modal deformations or adding the modal effective masses.

    These quantities can be saved in text form using:

    IMPR_RESU(FORMAT = 'RESULTAT',
              RESU = _F(RESULTAT = mode_0,
                        TOUT_PARA = 'OUI',
                        TOUT_CHAM = 'NON'),
              UNITE = 80)
    

TRI_DIAG

SORENSEN

MASS_EFFE_DY

MASS_EFFE_DZ

MASS_EFFE_DY

MASS_EFFE_DZ

Mode 1

3.16920

2.80830

1.13915

4.83834

Mode 2

2.80672

3.17077

4.83834

1.13915

TOTAL

5.97592

5.97907

5.97749

5.97749

20.1.1.4. Verification tests

The results obtained with the SORENSEN method are presented in the table below

Mode

Frequency in Hz (SORENSEN method)

1

1.67638E+02

2

1.67638E+02

3

1.05060E+03

4

1.05060E+03

5

1.48054E+03

6

2.59704E+03

7

2.94237E+03

8

2.94237E+03

9

4.47822E+03

10

5.76991E+03

11

5.76991E+03

20.1.2. Model B: Modal analysis of an unconstrained sphere

20.1.2.1. Description of the problem

20.1.2.1.1. Objective

The objective of this model is to determine the modes for a spherical structure in free-free mode (with the presence of multiple rigid body modes of rigid body). In particular, we will determine:

  • The natural frequencies in a given frequency band with the method of SORENSEN

  • The natural frequencies in a frequency band with the LANCZOS method, with or without rigid mode (OPTION = MODE_RIGIDE)

  • The 16 smallest natural frequencies with the method of SORENSEN

20.1.2.1.2. Geometry

../../_images/fig2_v2.08.011.png

20.1.2.1.3. Material properties

The material is linear elastic:

  • Young’s modulus: \(E\) = 100 MPa

  • Poisson’s ratio: \(\nu\) = 0.3

  • Mass density: \(\rho\) = 10,000 kg/m \(^3\)

20.1.2.1.4. Boundary conditions and loads

None.

20.1.2.2. Creating and running the model

20.1.2.2.1. Geometry and mesh creation with gmsh

../../_images/fig3_v2.08.011.png

Fig. 20.1.6 Sphere mesh created with gmsh

The geometry and the mesh in this example have been created with gmsh (https://gmsh.info/doc/texinfo/gmsh.html). The Python script for creating the model is in forma11b-geo.py (listed below).

#!/usr/bin/env python
# coding: utf-8

import gmsh
import sys
import math

gmsh.initialize()

# Choose Built-in kernel
model = gmsh.model
geo = model.geo
mesh = model.mesh
model.add("forma11b")


# Sphere size
el_size = 1
rad = 0.01
rad2 = rad/math.sqrt(3.0)
num_nodes = 5

# Add points
p0 = geo.addPoint(0, 0, 0, el_size)
p102 = geo.addPoint(rad2, rad2, -rad2, el_size)
p103 = geo.addPoint(-rad2, rad2, -rad2, el_size)
p104 = geo.addPoint(-rad2, -rad2, -rad2, el_size)
p105 = geo.addPoint(rad2, -rad2, -rad2, el_size)

# Add circles
l1 = geo.addCircleArc(p103, p0, p102)
l2 = geo.addCircleArc(p102, p0, p105)
l3 = geo.addCircleArc(p105, p0, p104)
l4 = geo.addCircleArc(p104, p0, p103)

# Create curve loop
loop1 = geo.addCurveLoop([1, 2, 3, 4])

# Create surface
surf1 = geo.addSurfaceFilling([loop1], sphereCenterTag = 1)

# Create and rotate copies (top surf)
surf2 = geo.copy([(2, surf1)])
geo.rotate(surf2, 0, 0, 0, 1, 0, 0, math.pi/2)

# Create and rotate copies (opposite surf)
surf3 = geo.copy([(2, surf1)])
geo.rotate(surf3, 0, 0, 0, 1, 0, 0, math.pi)

# Create and rotate copies (bot surf)
surf4 = geo.copy([(2, surf1)])
geo.rotate(surf4, 0, 0, 0, 1, 0, 0, 3*math.pi/2)

# Create and rotate copies (left surf)
surf5 = geo.copy([(2, surf1)])
geo.rotate(surf5, 0, 0, 0, 0, 1, 0, math.pi/2)

# Create and rotate copies (right surf)
surf6 = geo.copy([(2, surf1)])
geo.rotate(surf6, 0, 0, 0, 0, 1, 0, -math.pi/2)

# Create surface loop
sloop1 = geo.addSurfaceLoop([surf1, surf2[0][1], surf3[0][1], surf4[0][1], surf5[0][1], surf6[0][1]])

# Create volume
sphere = geo.addVolume([sloop1])

# Synchronize
geo.synchronize()

curve_list = [1, 2, 3, 4, 6, 7, 9, 11, 12, 14, 17, 19]
surface_list = [1, 5, 10, 15, 20, 21]
vol_list = [1]

# Set transfinite interpolation
for curve in curve_list:
    mesh.setTransfiniteCurve(curve, num_nodes)

for surf in surface_list:
    mesh.setTransfiniteSurface(surf)

for vol in vol_list:
    mesh.setTransfiniteVolume(vol)

# Set meshing options
gmsh.option.setNumber('Mesh.RecombineAll', 1)
gmsh.option.setNumber('Mesh.RecombinationAlgorithm', 1)
gmsh.option.setNumber('Mesh.Recombine3DLevel', 2)
gmsh.option.setNumber('Mesh.ElementOrder', 2)
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)

# Generate mesh
mesh.generate(3)

# Save mesh
gmsh.write("forma11b.msh")

# Activate GUI
gmsh.fltk.run()

# Clean up
gmsh.finalize()

The mesh contains 64 HEXA27 elements, and 730 nodes. Only the volume elements (3D) will be used in the model and the other elements can be deleted if you wish (see https://bbanerjee.github.io/ParSim/fem/cracks/python/salome-meca/code-aster/inspecting-and-manipulating-meshes/ for an example).

The mesh can then be converted into MED format using AsterStudy:

DEBUT(LANG='EN')

mesh = LIRE_MAILLAGE(FORMAT='GMSH',
                     INFO=2,
                     UNITE=20,
                     VERI_MAIL=_F(VERIF='OUI'))

IMPR_RESU(FORMAT='MED',
          RESU=_F(MAILLAGE=mesh),
          UNITE=80)

IMPR_RESU(FORMAT='ASTER',
          RESU=_F(MAILLAGE=mesh),
          UNITE=2)

FIN()

20.1.2.2.2. Creation of the command file in AsterStudy

  • Read of the MED mesh (LIRE_MAILLAGE)

  • Define the finite element type (AFFE_MODELE) and use 3D elements (3D)

  • Define and assign the material (DEFI_MATERIAU and AFFE_MATERIAU). You will need to define RHO so that the mass matrix can be computed.

  • Compute the element stiffness and mass matrices with the macro ASSEMBLAGE. Name the stiffness matrix stif_mat, the mass matrix mass_mat, and the DOF vector num_dof. The macro will create the assembled matrices.

  • Solve the modal analysis problem with CALC_MODES [U4.52.02]

You can try the following:

  • Calculate with the SORENSEN method (keyword SOLVEUR_MODAL) the frequencies in the frequency band (0 Hz, 2880 Hz) as well as the associated modes

  • If the calculation fails, the frequency band can be extended to have a small negative value for the lower bound

  • Print the eigen modes (IMPR_RESU) in MED format MED for visualization

  • Calculate with the TRI_DIAG (Lanczos) method (with or without option MODE_RIGIDE) the frequencies in the band (0 Hz, 2880 Hz) as well as the associated modes

  • Calculate with the SORENSEN method the 16 smallest frequencies as well as the associated modal shapes. The PREC_SHIFT parameter can be used (under the keyword factor CALC_FREQ) to get around the problem of zero frequencies.

The annotated Python command file for Model B is forma11b.comm. The contents of the file are listed below for convenience.

# Modal analysis of a solid sphere
DEBUT(LANG = 'EN')

# Read mesh
mesh = LIRE_MAILLAGE(FORMAT = 'MED')

# Define and assign material
elastic = DEFI_MATERIAU(ELAS= _F(E = 1.E8,
                                 NU = 0.3,
                                 RHO = 1.E4))

mater = AFFE_MATERIAU(MAILLAGE = mesh,
                      AFFE = _F(TOUT = 'OUI',
                                MATER = elastic))

# Assign 3D elements to mesh
model = AFFE_MODELE(MAILLAGE = mesh,
                    AFFE = _F(TOUT = 'OUI',
                              PHENOMENE = 'MECANIQUE',
                              MODELISATION = '3D'))

# Create stiffness and mass matrices
ASSEMBLAGE(CHAM_MATER = mater,
           MATR_ASSE = (_F(MATRICE = CO('stif_mat'),
                           OPTION = 'RIGI_MECA'),
                        _F(MATRICE = CO('mass_mat'),
                           OPTION = 'MASS_MECA')),
           MODELE = model,
           NUME_DDL = CO('num_dof'))

# Modal analysis 1 (BANDE/SORENSEN)
mode_sob = CALC_MODES(OPTION='BANDE',
                      SOLVEUR_MODAL = _F(METHODE = 'SORENSEN',
                                         COEF_DIM_ESPACE = 8),
                      SOLVEUR = _F(NPREC = 9),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (-0.3, 2700.0)),
                      TITRE = '1 - SORENSEN/BANDE')

# Modal analysis 2 (BANDE/TRI_DIAG) (Lanzcos)
mode_lan = CALC_MODES(OPTION = 'BANDE',
                      SOLVEUR_MODAL = _F(MODE_RIGIDE = 'OUI',
                                         METHODE = 'TRI_DIAG'),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (-0.0, 2800.0)),
                      TITRE = '2 - LANCZOS/BANDE')

# Modal analysis 3 (PLUS_PETITE/SORENSEN) 
# 16 smallest modes
# The threshold frqeuency (SEUIL_FREQ) is used to avoid
# rigid body modes
mode_sop = CALC_MODES(OPTION = 'PLUS_PETITE',
                      SOLVEUR_MODAL = _F(METHODE = 'SORENSEN'),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(NMAX_FREQ = 16,
                                     SEUIL_FREQ = 1.0),
                      TITRE = '3 - SORENSEN/PLUs_PETITE')

# Modal analysis 4 (BANDE/TRI_DIAG) (Lanzcos - no rigid body modes)
mode_las = CALC_MODES(OPTION = 'BANDE',
                      SOLVEUR_MODAL = _F(METHODE = 'TRI_DIAG',
                                         MODE_RIGIDE = 'NON'),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (-0.0, 2800.0)),
                      VERI_MODE = _F(STOP_ERREUR = 'NON'),
                      TITRE = '4 - LANCZOS/BANDE with MODE_RIGIDE')

# Modal analysis 5 (BANDE/QZ) 
mode_qz = CALC_MODES(OPTION = 'BANDE',
                     SOLVEUR_MODAL = _F(METHODE = 'QZ',
                                        TYPE_QZ = 'QZ_QR'),
                     MATR_RIGI = stif_mat,
                     MATR_MASS = mass_mat,
                     CALC_FREQ = _F(FREQ = (0.0, 2700.0)),
                     TITRE = '5 - QZ/BANDE')

# Write results in table format 
IMPR_RESU(FORMAT = "RESULTAT",
          RESU = _F(RESULTAT = mode_sob,
                    TOUT_CHAM = 'NON',
                    TOUT_PARA = 'OUI'))

IMPR_RESU(FORMAT = "RESULTAT",
          RESU = _F(RESULTAT = mode_lan,
                    TOUT_CHAM = 'NON',
                    TOUT_PARA = 'OUI'))

IMPR_RESU(FORMAT = "RESULTAT",
          RESU = _F(RESULTAT = mode_sop,
                    TOUT_CHAM = 'NON',
                    TOUT_PARA = 'OUI'))

IMPR_RESU(FORMAT = "RESULTAT",
          RESU = _F(RESULTAT = mode_las,
                    TOUT_CHAM = 'NON',
                    TOUT_PARA = 'OUI'))

IMPR_RESU(FORMAT = "RESULTAT",
          RESU = _F(RESULTAT = mode_qz,
                    TOUT_CHAM = 'NON',
                    TOUT_PARA = 'OUI'))

# Write results in MED format 
IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(MAILLAGE = mesh)))

IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(RESULTAT = mode_sob)))

IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(RESULTAT = mode_lan)))

IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(RESULTAT = mode_sop)))

IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(RESULTAT = mode_las)))

IMPR_RESU(FORMAT = 'MED',
          UNITE = 80,
          RESU = (_F(RESULTAT = mode_qz)))

# Verification tests
# The SORENSEN zero mode: may issue warning.  First non-zero mode is 2613 Hz.
TEST_RESU(RESU = _F(NUME_ORDRE = 1,
                    PARA = 'FREQ',
                    RESULTAT = mode_sob,
                    VALE_CALC = 0.,
                    VALE_REFE = 0.,
                    PRECISION = 5.E-03,
                    REFERENCE = 'NON_DEFINI',
                    CRITERE = 'ABSOLU'))

TEST_RESU(RESU = _F(NUME_ORDRE = 7,
                    PARA = 'FREQ',
                    RESULTAT = mode_sob,
                    VALE_CALC = 2469.64060726,
                    CRITERE = 'RELATIF'))

# The QZ zero mode: may issue warning.  First non-zero mode is 2613 Hz.
TEST_RESU(RESU = _F(NUME_ORDRE = 1,
                    PARA = 'FREQ',
                    REFERENCE = 'NON_DEFINI',
                    RESULTAT = mode_qz,
                    VALE_CALC = 0.0,
                    VALE_REFE = 0.0,
                    CRITERE = 'ABSOLU',
                    PRECISION = 3.5E-3))

TEST_RESU(RESU = _F(NUME_ORDRE = 16,
                    PARA = 'FREQ',
                    RESULTAT = mode_qz,
                    VALE_CALC = 2613.64123957,
                    CRITERE = 'RELATIF'))


TEST_RESU(RESU = _F(NUME_ORDRE = 11,
                    PARA = 'FREQ',
                    RESULTAT = mode_lan,
                    VALE_CALC = 2470.84079534))


# Finish
FIN();

20.1.2.2.3. Postprocessing with Salome-Meca

To visualize the modal deformations with Salomé-Meca, it is convenient to use the Results tab in AsterStudy. Some of the modes from the BANDE/SORENSEN algorithm are shown below.

../../_images/forma11b_sob_1.png

Fig. 20.1.7 Mode 1 - negative

../../_images/forma11b_sob_7.png

Fig. 20.1.8 Mode 7 - 2470 Hz

../../_images/forma11b_sob_10.png

Fig. 20.1.9 Mode 10 - 2473 Hz

../../_images/forma11b_sob_14.png

Fig. 20.1.10 Mode 14 - 2616 Hz

20.1.2.3. Remarks

Factoring of the shifted matrices \((\mathbf{K} - \sigma \mathbf{I}) = \mathbf{L} \mathbf{D} \mathbf{L}^T\) governed by the settings: NMAX_ITER_SHIFT, SOLVEUR / NPREC [U4.52.02] and PREC_SHIFT, can be carried out only if these are regular.

This often poses a problem when the magnitude of components of \(\mathbf{K}\) is large compared to that of the components of \(\mathbf{M}\) and \(\sigma\) is a good approximate eigenvalue. The Code_Aster policy is to emit an ALARM when it is about a Sturm test and an ERREUR_FATALE when it is about the working matrix. In the event of a problem, you can always change the option or the shift.

The detection of rigid body modes is often problematic for classical modal solvers. It is good practice to use the SORENSEN method with a band of frequencies whose lower bound is zero, or even slightly negative. The solver normally captures the rigid body modes without any problem (it’s just a special multiple mode case). Otherwise, with TRI_DIAG, the MODE_RIGIDE option is recommended.

20.1.2.4. Verification tests

The results obtained with the SORENSEN method and the TRI_DIAG (LANCZOS) method with and without the option OPTION = MODE_RIGIDE are presented in the table below.

Warning

These values are obtained with the mesh that is provide in the default Code-Aster installation and are slightly different from those computed with the gmsh HEAX27 element mesh.

Sorensen method (METHOD = ‘SORENSEN’)

Lanczos method (with MODE_RIGIDE)

Lanczos method (without MODE_RIGIDE)

Mode

Frequency

Mode

Frequency

Mode

Frequency

5

-3.88562E-05

6

0.00000E+00

4

1.22874E-04

4

-3.88562E-05

5

0.00000E+00

5

2.70602E-04

3

-6.73009E-05

4

0.00000E+00

3

-5.03634E-04

2

1.09902E-04

3

0.00000E+00

2

-5.69081E-04

6

1.25908E-04

2

0.00000E+00

6

4.46762E-03

1

1.73770E-04

1

0.00000E+00

7

2.46964E+03

7

2.46964E+03

7

2.46964E+03

8

2.46964E+03

8

2.46964E+03

8

2.46964E+03

9

2.47084E+03

9

2.47084E+03

9

2.47084E+03

10

2.47084E+03

10

2.47084E+03

10

2.47084E+03

11

2.47084E+03

11

2.47084E+03

11

2.47084E+03

12

2.61345E+03

12

2.61345E+03

12

2.61345E+03

13

2.61345E+03

13

2.61345E+03

13

2.61345E+03

14

2.61345E+03

14

2.61345E+03

14

2.61345E+03

15

2.61364E+03

15

2.61364E+03

15

2.61364E+03

16

2.61364E+03

16

2.61364E+03

16

2.61364E+03

17

3.48770E+03

20.1.3. Model C: Modal analysis of a cooling tower

../../_images/fig4_v2.08.011.png

20.1.3.1. Description of the problem

In this tutorial we will compute the fundamental modes of a large cooling tower made of a concrete-like material. We will model the tower using DKT shell elements.

The following actions will be performed:

  • Compute modes in specified frequency bands by simultaneous iteration

  • Apply of a filtering criterion for the modes (based on checking if a modal parameter value is greater than a threshold).

  • Concatenate the calculated data structures into one.

The modes are calculated with the command CALC_MODES [U4.52.02] with the option ‘BANDE’. They are then normalized with the command NORM_MODE [U4.52.11]. The calculated modes are filtered and concatenated with the command EXTR_MODE [U4.52.12].

20.1.3.1.1. Objective

The objective of this model is:

  • To extract geometric/physical quantities from the finite element model

  • To evaluate the number of modes

  • To use macro CALC_MODES with the option ‘BANDE’ (with the input frequency range divided into several sub-bands) to calculate the eigen frequencies.

20.1.3.1.2. Geometry and mesh

A mesh of the problem is provided in the test case directory in the Code_Aster installation (.../share/aster/tests). We will not use that mesh directly but create our own version of the geometry based on that model.

The thickness of the shell of the cooling tower is 0.305 m.

Extracting the coordinates of the tower profile

We will use the Code_Aster mesh file forma11c.mmed to extract the coordinates of the profile and save them in JSON format (forma11c_profile.json). The Python script for this process (forma11c_profile.py) can be run in Salome-Meca and is listed below.

#!/usr/bin/env python
import sys
import salome
import GEOM
import SMESH
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import math
import json

ExportPATH="/home/banerjee/Salome/forma11/"

salome.salome_init()
geompy = geomBuilder.New()
smesh = smeshBuilder.New()

# Read mesh
([mesh], status) = smesh.CreateMeshesFromMED(r'/home/banerjee/Salome/forma11/forma11c.mmed')

# Create a profile by selecting nodes by coordinate
profile = mesh.CreateEmptyGroup( SMESH.NODE, 'profile' )

# Create selection box
eps = 1.0e-3
mesh_bb = mesh.BoundingBox()
x_min = 0
x_max = mesh_bb[3] + eps
z_min = mesh_bb[2] - eps
z_max = mesh_bb[5] + eps
y_min = -eps
y_max =  eps

box_min = geompy.MakeVertex(x_min, y_min, z_min)
box_max = geompy.MakeVertex(x_max, y_max, z_max)
box = geompy.MakeBoxTwoPnt(box_min, box_max)
geompy.addToStudy(box, "selection_box")

# Get the nodes that are in the selection box
filter = smesh.GetFilter(SMESH.NODE,  SMESH.FT_BelongToGeom, box)
ids = mesh.GetIdsFromFilter(filter)

# Get the nodal coordinates
coords = []
for id in ids:
   coords.append(mesh.GetNodeXYZ(id))
   nbAdd = profile.Add([id])
  
# Save the coordinates (string)
file_id = open( r''+ExportPATH+'forma11c_profile.coords'+'', 'w')
for coord in coords:
    file_id.write(str(coord) + '\n')

file_id.close()

# Save the coordinates (json)
file_id = open( r''+ExportPATH+'forma11c_profile.json'+'', 'w')
json.dump(coords, file_id)
file_id.close()


# Update GUI tree
if salome.sg.hasDesktop():
  salome.sg.updateObjBrowser()

One the profile coordinates are available, the geometry and mesh can be generated either with Salome-Meca or with gmsh.

Geometry and mesh with Salome-Meca

The Python script for generating the geometry and mesh in MED format with Salome-Meca is forma11c_salome.py and is listed below. The model is meshed with the medial axis algorithm to produce a fully quadrilateral mesh.

Note

Note that the bottom edge of the cooling tower (which is fixed) is assigned the edge geometry group GM4 and the faces of the tower are assigned the geometry group GM1. These names have been used to match those produce when converting the gmsh mesh to MED format with Code_Aster.

#!/usr/bin/env python
import salome
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import GEOM, SMESH
import math
import json

ExportPATH="/home/banerjee/Salome/forma11/"

# Read the profile coordinates
file_id = open( r''+ExportPATH+'forma11c_profile.json'+'', 'r')
coords = json.load(file_id)
file_id.close()

salome.salome_init()

#-------------------
# Geometry
geompy = geomBuilder.New()

O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)

v_profile = []
for coord in coords:
    v = geompy.MakeVertex(coord[0], coord[1], coord[2])
    v_profile.append(v)

curve = geompy.MakeInterpol(v_profile, False, True)
surf = geompy.MakeRevolution2Ways(curve, OZ, math.pi)
geompy.addToStudy(surf, 'tower' )

# Geometry groups
surf_ID = geompy.GetSubShapeID(surf, surf)
g_surf = geompy.CreateGroup(surf, geompy.ShapeType["FACE"])
geompy.AddObject(g_surf, surf_ID)
id_g_surf = geompy.addToStudy(g_surf, 'GM1')

surf_bb = geompy.BoundingBox(surf)
x_max = surf_bb[1]
z_min = surf_bb[4]
pt = geompy.MakeVertex(x_max, 0, z_min)
bot_edge = geompy.GetEdgeNearPoint(surf, pt)
bot_edge_ID = geompy.GetSubShapeID(surf, bot_edge)
g_bot_edge = geompy.CreateGroup(surf, geompy.ShapeType["EDGE"])
geompy.AddObject(g_bot_edge, bot_edge_ID)
id_g_bot_edge = geompy.addToStudy(g_bot_edge, 'GM4')


#-------------------
# Mesh
smesh = smeshBuilder.New()

# Initialize mesh
mesh = smesh.Mesh(surf)
smesh.SetName(mesh.GetMesh(), 'mesh')

# Algorithms
Regular_1D = mesh.Segment()
QuadFromMedialAxis_1D2D = mesh.Quadrangle(algo = smeshBuilder.QUAD_MA_PROJ)
Number_of_Segments = Regular_1D.NumberOfSegments(60)
Number_of_Layers = QuadFromMedialAxis_1D2D.NumberOfLayers(31)

smesh.SetName(Number_of_Layers, 'Number of Layers')
smesh.SetName(Number_of_Segments, 'Number of Segments')

# Compute mesh
isDone = mesh.Compute()

# Map geometry groups to mesh
base = mesh.GroupOnGeom(g_bot_edge, 'GM4', SMESH.EDGE)
tower = mesh.GroupOnGeom(g_surf, 'GM1', SMESH.FACE)

# Export mesh
mesh.ExportMED( r''+ExportPATH+'forma11c_salome.mmed'+'')

# Update GUI
if salome.sg.hasDesktop():
  salome.sg.updateObjBrowser()

Geometry and mesh with gmsh

Alternatively, we can create the geometric model and the mesh with gmsh. The Python script for this process is at forma11c_gmsh.py. A listing is provided below for convenience. The model is mesh with transfinite interpolation followed by element recombination.

#!/usr/bin/env python
# coding: utf-8

import gmsh
import sys
import math
import json

ExportPATH="/home/banerjee/Salome/forma11/"

gmsh.initialize()


# Choose OpenCascade kernel
model = gmsh.model
occ = model.occ
mesh = model.mesh
model.add("forma11c_gmsh")

# Default size
el_size = 0.7

# Read the profile coordinates
file_id = open( r''+ExportPATH+'forma11c_profile.json'+'', 'r')
coords = json.load(file_id)
file_id.close()

# Create profile points
v_profile = []
for coord in coords:
    v = occ.addPoint(coord[0], coord[1], coord[2], el_size)
    v_profile.append(v)

# Create spline going through profile points
l1 = occ.addBSpline(v_profile)

# Create copies and rotate
l2 = occ.copy([(1,l1)])
l3 = occ.copy([(1,l1)])
l4 = occ.copy([(1,l1)])

# Rotate the copies
occ.rotate(l2, 0, 0, 0, 0, 0, 1, math.pi/2)
occ.rotate(l3, 0, 0, 0, 0, 0, 1, math.pi)
occ.rotate(l4, 0, 0, 0, 0, 0, 1, 3*math.pi/2)

# Sweep the lines
surf1 = occ.revolve([(1, l1)], 0, 0, 0, 0, 0, 1, math.pi/2)
surf2 = occ.revolve(l2, 0, 0, 0, 0, 0, 1, math.pi/2)
surf3 = occ.revolve(l3, 0, 0, 0, 0, 0, 1, math.pi/2)
surf4 = occ.revolve(l4, 0, 0, 0, 0, 0, 1, math.pi/2)

# Join the surfaces
surf5 = occ.fragment(surf1, surf2)
surf6 = occ.fragment(surf3, surf4)
surf7 = occ.fragment(surf5[0], surf6[0])

# Update
occ.synchronize()


# Get all geometric entities
faces = [face[1] for face in occ.getEntities(2)]
edges = [curve[1] for curve in occ.getEntities(1)]         
vertices = [point[1] for point in occ.getEntities(0)]

# Add physical groups for the mesh
chimney_faces = model.addPhysicalGroup(2, faces)
model.setPhysicalName(2, chimney_faces, "chimney_faces")
chimney_edges = model.addPhysicalGroup(1, edges)
model.setPhysicalName(1, chimney_edges, "chimney_edges")
chimney_verts = model.addPhysicalGroup(0, vertices)
model.setPhysicalName(0, chimney_verts, "chimney_verts")

# Add physical group for the base edges
bb1 = model.getBoundingBox(2, 1)
bb2 = model.getBoundingBox(2, 2)
bb3 = model.getBoundingBox(2, 3)
bb4 = model.getBoundingBox(2, 4)
xmax = bb1[3]
xmin = bb2[0]
ymax = xmax
ymin = xmin
zmin = bb1[2]
zmax = zmin

eps = 1.0e-3
base_edges = model.getEntitiesInBoundingBox(xmin-eps, ymin-eps, zmin-eps, xmax+eps, ymax+eps, zmax+eps, dim=1)
edges = [curve[1] for curve in base_edges]
chimney_base = model.addPhysicalGroup(1, edges)
model.setPhysicalName(1, chimney_base, "chimney_base")

# Assign transfinite interpolation to enetities
num_nodes_circ = 15
for curve in occ.getEntities(1):
    mesh.setTransfiniteCurve(curve[1], num_nodes_circ)

num_nodes_vert = 32
vertical_curves = [7, 10, 13, 17]
for curve in vertical_curves:
    mesh.setTransfiniteCurve(curve, num_nodes_vert)

for surf in occ.getEntities(2):
    mesh.setTransfiniteSurface(surf[1])

# Set up meshing options
gmsh.option.setNumber('Mesh.RecombineAll', 1)
gmsh.option.setNumber('Mesh.RecombinationAlgorithm', 1)
gmsh.option.setNumber('Mesh.Recombine3DLevel', 2)
gmsh.option.setNumber('Mesh.ElementOrder', 1)
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)

# Generate mesh
mesh.generate(2)

# Save mesh
gmsh.write(r''+ExportPATH+'forma11c_gmsh.msh'+'')

# Check in GUI
gmsh.fltk.run()

# Clean up
gmsh.finalize()

The mesh can be converted into MED format using AsterStudy using the commands:

DEBUT(LANG='EN')

mesh = LIRE_MAILLAGE(FORMAT='GMSH',
                     INFO=2,
                     UNITE=20,
                     VERI_MAIL=_F(VERIF='OUI'))

IMPR_RESU(FORMAT='MED',
          RESU=_F(MAILLAGE=mesh),
          UNITE=80)

IMPR_RESU(FORMAT='ASTER',
          RESU=_F(MAILLAGE=mesh),
          UNITE=2)

FIN()

20.1.3.1.3. Material Properties

The material is linear elastic:

  • Young’s modulus: \(E\) = 27.6 GPa

  • Poisson’s ratio: \(\nu\) = 0.166

  • Mass density: \(\rho\) = 2244 kg/m \(^3\)

20.1.3.1.4. Boundary conditions

The base of the tower is fixed.

20.1.3.2. Characteristics of model

20.1.3.2.1. Characteristics of the mesh

The mesh produced by Salome-Meca is composed of 1860 QUAD4 elements and 151 edges. There are 1920 nodes.

The mesh generated by gmsh contains 1736 QUAD4 element, 292 edges, and 37 vertices. There are 1821 nodes in the mesh.

The tower is modeled with DKT shell elements

20.1.3.2.2. Aster command file

The main commands needed for the Aster computation are:

  • Read of the MED mesh (LIRE_MAILLAGE)

  • Define the finite element type (AFFE_MODELE) and use shell elements DKT

  • Assign the characteristics of the shell elements (AFFE_CARA_ELEM). The thickness of all shell elements is the same. Use EPAIS = 0.305 (THICKNESS = 0.305).

  • Define and assign the material (DEFI_MATERIAU and AFFE_MATERIAU). You will need to define RHO so that the mass matrix can be computed.

  • Assign the boundary conditions (AFFE_CHAR_MECA):

    • Constrain the base (GM4) with DDL_IMPO (fix the displacement and rotation DOFs)

  • Compute the element stiffness and mass matrices with the macro ASSEMBLAGE. Name the stiffness matrix stif_mat (option RIGI_MECA), the mass matrix mass_mat (option MASS_MECA), and the DOF vector num_dof (option NUME_DDL). The macro will create the assembled matrices.

  • Write the model parameters (POST_ELEM) to a table and write it to disk with IMPR-TABLE

  • Extract the modes between 0 and 4 Hz with INFO_MODE.

  • Solve the modal analysis problem with CALC_MODES [U4.52.02]

For example, you can explore the following:

  • Calculate the eigen frequencies and the associated modes present in the band of frequency 0 Hz to 4 Hz (CALC_MODES)

  • Normalize with the infinity norm (NORM_MODE), on all the DOFs (TRAN_ROTA) of the physical nodes

  • Extract modes based on the unit effective mass criterion (EXTR_MODE / MASS_EFFE_UN).

  • Print the eigen modes (IMPR_RESU) in MED format for visualization.

The first three commands above can be combined into a single command.

  • Recalculate the eigen frequencies with option BANDE after combining the three operations CALC_MODES, NORM_MODE and EXTR_MODE in the CALC_MODES command. Divide CALC_FREQ into 2 intervals.

  • Redo the above step with 4 subdivisions of CALC_FREQ

  • Redo the previous step but use the option TRI_DIAG (Lanczos method) this time.

The command file for these studies can be downloaded from forma11c.comm. It is listed below for convenience.

# Modal analysis of chimney shell
DEBUT(LANG = 'EN')

# Read mesh
mesh = LIRE_MAILLAGE(FORMAT = 'MED')

# Assign model
model = AFFE_MODELE(MAILLAGE = mesh,
                    AFFE = _F(GROUP_MA = 'GM1',
                              PHENOMENE = 'MECANIQUE',
                              MODELISATION = 'DKT'))

# Define and assign material
concrete = DEFI_MATERIAU(ELAS = _F(E = 2.76E10,
                                   NU = 0.166,
                                   RHO = 2244.0))

mater = AFFE_MATERIAU(MAILLAGE = mesh,
                      AFFE = _F(GROUP_MA = 'GM1',
                                MATER = concrete))

# Assign shell properties
shell = AFFE_CARA_ELEM(MODELE = model,
                       COQUE = _F(GROUP_MA = 'GM1',
                                  EPAIS = 0.305))

# Displcamenet BC
disp_bc = AFFE_CHAR_MECA(MODELE = model,
                         DDL_IMPO = _F(GROUP_MA = 'GM4',
                                       DX = 0.0,
                                       DY = 0.0,
                                       DZ = 0.0,
                                       DRX = 0.0,
                                       DRY = 0.0,
                                       DRZ = 0.0))

# Create stiffness and mass matrices
ASSEMBLAGE(CHAM_MATER = mater,
           CARA_ELEM = shell,
           CHARGE = disp_bc,
           MATR_ASSE = (_F(MATRICE = CO('stif_mat'),
                           OPTION = 'RIGI_MECA'),
                        _F(MATRICE = CO('mass_mat'),
                           OPTION = 'MASS_MECA')),
           MODELE = model,
           NUME_DDL = CO('num_dof'))

# Write model parameters to a table
model_p = POST_ELEM(MODELE = model,
                    CHAM_MATER = mater,
                    CARA_ELEM = shell,
                    MASS_INER = _F(TOUT = 'OUI'))
IMPR_TABLE(TABLE = model_p,
           SEPARATEUR = ',')

# Compute number of modes btween 0 and 4 Hz
num_mode = INFO_MODE(MATR_RIGI = stif_mat,
                     MATR_MASS = mass_mat,
                     FREQ = (0.,4.0,))

# Modal analysis 1 : BANDE/SORENSEN
mode_sor = CALC_MODES(OPTION = 'BANDE',
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (0.0, 4.0)))

# Normalize modes with L-infinity norm on all nodal components
# required for the computation of effective unit masses
mode_nrm = NORM_MODE(MODE = mode_sor,
                     NORME = 'TRAN_ROTA')

# Sort the modes
mode_srt = EXTR_MODE(FILTRE_MODE = _F(MODE = mode_nrm,
                                      CRIT_EXTR = 'MASS_EFFE_UN',
                                      SEUIL = 0.001),
                     IMPRESSION = _F(CUMUL = 'OUI'))

# Model analysis 2: Perform the three operations above in two steps
# using a single command
mode_sr2 = CALC_MODES(OPTION = 'BANDE',
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (0.0, 2.0, 4.0)),
                      STOP_BANDE_VIDE = 'OUI',
                      NORM_MODE = _F(NORME = 'TRAN_ROTA'),
                      FILTRE_MODE =_F(CRIT_EXTR = 'MASS_EFFE_UN'),
                      IMPRESSION = _F(CUMUL = 'OUI'))

# Model analysis 3: Perform the three operations above in four steps
# using a single command with a MULT_FRONT solver
mode_sr4 = CALC_MODES(OPTION = 'BANDE',
                      SOLVEUR = _F(METHODE = 'MULT_FRONT'),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (0.0, 1.0, 2.0, 3.0, 4.0)),
                      STOP_BANDE_VIDE = 'OUI',
                      NORM_MODE = _F(NORME = 'TRAN_ROTA'),
                      FILTRE_MODE =_F(CRIT_EXTR = 'MASS_EFFE_UN'),
                      IMPRESSION = _F(CUMUL = 'OUI'))

# Model analysis 4: Perform the three operations above in four steps
# using a single command with a Lanzcos method
mode_ln4 = CALC_MODES(OPTION = 'BANDE',
                      SOLVEUR_MODAL = _F(METHODE = 'TRI_DIAG',
                                         DIM_SOUS_ESPACE = 180),
                      VERI_MODE = _F(SEUIL = 1e-05),
                      MATR_RIGI = stif_mat,
                      MATR_MASS = mass_mat,
                      CALC_FREQ = _F(FREQ = (0.0, 1.0, 2.0, 3.0, 4.0)),
                      STOP_BANDE_VIDE = 'OUI',
                      NORM_MODE = _F(NORME = 'TRAN_ROTA'),
                      FILTRE_MODE =_F(CRIT_EXTR = 'MASS_EFFE_UN'),
                      IMPRESSION = _F(CUMUL = 'OUI'))
                          
# Write results in MED format
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(MAILLAGE = mesh,
                    RESULTAT = mode_sor))
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(RESULTAT = mode_nrm))
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(RESULTAT = mode_srt))
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(RESULTAT = mode_sr2))
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(RESULTAT = mode_sr4))
IMPR_RESU(FORMAT = 'MED',
          RESU = _F(RESULTAT = mode_ln4))

# Verification tests
# Check the number of modes
TEST_TABLE(VALE_CALC_I = 86,
           NOM_PARA = 'NB_MODE',
           TABLE = num_mode)

# Check the first mode
TEST_RESU(RESU = _F(NUME_ORDRE = 1,
                  PARA = 'FREQ',
                  RESULTAT = mode_nrm,
                  VALE_CALC = 1.1385562506602))

# Finish
FIN();

20.1.3.2.3. Postprocessing

The easiest way to explore the results is through the AsterStudy Results tab. Some of the mode shapes are shown below.

../../_images/forma11c_sor_1_13.png

Fig. 20.1.11 Mode 1 - 1.13 Hz

../../_images/forma11c_sor_1.png

Fig. 20.1.12 Mode 47 - 2.80 Hz

../../_images/forma11c_sor_2.png

Fig. 20.1.13 Mode 48 - 2.80 Hz

../../_images/forma11c_sor_3_84.png

Fig. 20.1.14 Mode 60 - 3.84 Hz

20.1.3.3. Remarks

Good practices:

  • Assess the model with POST_ELEM

  • Evaluate the spectrum with INFO_MODE.

  • Use CALC_MODES, with the option BANDE subdivided into several sub-bands. This approach is more economical if you want explore a broad spectrum.

  • Check the concordance of the identified ends of spectrum

  • Explore is there is any time saved by combining operations: normalization, filtering, concatenation of data structures etc.

20.1.3.4. Verification

The two modes that are filtered out correspond to the case of unit effective mass. The reference values of both these frequencies (multiplicity) is 12.80065 Hz.