テスト済み有限要素コード¶

CalculiX
Code_Aster
Elmer
MoFEM
GetFEM

CalculiX ¶

CalculiXは、剛体のボックス外への移動をサポートします。以下は入力ファイルの説明です。

** Mesh ++++++++++++++++++++++++++++++++++++++++++++++++

*INCLUDE, INPUT=Mesh/TET-0-5.inp # Path to mesh for ccx solver

** Material ++++++++++++++++++++++++++++++++++++++++++++++++

*Material, Name=Material-1	# Definition of material
*Density			# Definition of density
7.829E-09
*Elastic			# Definition of Young modulus and Poisson's ratio
207000, 0.33

** Sections ++++++++++++++++++++++++++++++++++++++++++++++++
** # Assigning material and solid elements to the elements sets in mesh
*Solid section, Elset=Fork-El_Section, Material=Material-1	

** Steps +++++++++++++++++++++++++++++++++++++++++++++++++++
** # Definition of frequency step with 12 eigenmodes
*Step
*Frequency
12

** Field outputs +++++++++++++++++++++++++++++++++++++++++++
** # Commands responsible for saving results
*Node file
RF, U
*El file
S, E

** End step ++++++++++++++++++++++++++++++++++++++++++++++++

*End step

この解析で使用するシミュレーション入力ファイルは、 GitHub にあります!

Code_Aster ¶

残念ながら、Code_Asterはほぼ0の周波数を得るための理想的なソルバではないので、BANDEオプションを使用して特定の帯域幅の周波数を計算しました。以下に、Code_Asterの入力ファイルについて説明します。

DEBUT()


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

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

mater = DEFI_MATERIAU(ELAS=_F(E=207000.0,                   # Defining the material
                              NU=0.33,
                              RHO=7.829e-09))

fieldmat = AFFE_MATERIAU(AFFE=_F(MATER=(mater, ),           # Assigning the material to the mesh
                                 TOUT='OUI'),
                         MAILLAGE=mesh,
                         MODELE=model)

ASSEMBLAGE(CHAM_MATER=fieldmat,                             # Setting the right matrixes for simulation
           MATR_ASSE=(_F(MATRICE=CO('MASS'),
                         OPTION='MASS_MECA'),
                      _F(MATRICE=CO('RIGI'),
                         OPTION='RIGI_MECA')),
           MODELE=model,
           NUME_DDL=CO('ndl'))




modes = CALC_MODES(CALC_FREQ=_F(FREQ=(100.0, 3000.0)),      # Setting the solver BANDE option between 100 and 3000 Hz.
                   MATR_MASS=MASS,
                   MATR_RIGI=RIGI,
                   OPTION='BANDE',
                   SOLVEUR_MODAL=_F(METHODE='TRI_DIAG',
                                    ),
                   TYPE_RESU='DYNAMIQUE')


IMPR_RESU(FORMAT='MED',                                     # Saving the result in 3D format
          RESU=_F(INFO_MAILLAGE='OUI',
                  RESULTAT=modes),
          UNITE=80)

IMPR_RESU(FORMAT='RESULTAT',                                # Saving the result in .dat file
          RESU=_F(RESULTAT=modes),
          UNITE=8)


FIN()

この解析で使用するシミュレーション入力ファイルは、 GitHub にあります!

Elmer ¶

Elmerは、剛体のボックス外へのモードをサポートします。以下はソルバー入力ファイルの説明です。

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "Mesh/QUAD-HEX-0-5" # Path to the mesh
  Include Path ""
  Results Directory "Results"     # Path to results directory
End

Simulation                        # Settings and constants for simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals = 1
  Timestepping Method = BDF
  BDF Order = 1
  Solver Input File = case.sif
  Post File = case.ep
End

Constants
  Gravity(4) = 0 -1 0 9.82
  Stefan Boltzmann = 5.67e-08
  Permittivity of Vacuum = 8.8542e-12
  Boltzmann Constant = 1.3807e-23
  Unit Charge = 1.602e-19
End

Body 1                            # Assigning the material and equations to the mesh
  Target Bodies(1) = 1
  Name = "Body 1"
  Equation = 1
  Material = 1
End

Solver 1                          # Solver settings
  Equation = Linear elasticity
  Procedure = "StressSolve" "StressSolver"
  Eigen System Select = Smallest magnitude
  Eigen System Values = 10
  Variable = -dofs 3 Displacement
  Eigen Analysis = True
  Exec Solver = Always
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Linear System Convergence Tolerance = 1.0e-3
  Steady State Convergence Tolerance = 1.0e-5
  Nonlinear System Convergence Tolerance = 1.0e-7
  Nonlinear System Max Iterations = 1
  Nonlinear System Newton After Iterations = 3
  Nonlinear System Newton After Tolerance = 1.0e-3
  Nonlinear System Relaxation Factor = 1
  Linear System Solver = Direct
  Linear System Direct Method = MUMPS
End

Solver 2                          # Extracting the result to the .dat file
  Exec Solver = After Saving
  Equation ="Equation 1"
  Procedure = "SaveData" "SaveScalars"
  Filename = "QUAD-HEX-0-5.dat"
  Variable 1 = Displacement
  Save Eigenfrequencies = Logical True
End

Equation 1                        # Setting the active solvers
  Name = "Equation 1"
  Active Solvers(2) = 1 2
End

Material 1                        # Defining the material
  Name = "Steel"
  Mesh Poisson ratio = 0.33
  Youngs modulus = 207000
  Poisson ratio = 0.33
  Porosity Model = Always saturated
  Density = 7.829e-9
End

この解析で使用するシミュレーション入力ファイルは、 GitHub にあります!

MoFEM ¶

MoFEMモデルの作成方法を確認したい場合は、MoFEMチュートリアルの VEC-1: Eigen elastic セクションを参照することをお勧めします。

GetFEM ¶

"""

"""
###############################################################################

import time
import getfem as gf
import numpy as np
import scipy
from scipy import io
from scipy.sparse import linalg

###############################################################################
# Let us now define the different physical and numerical parameters of the problem.
#

rho = 7.829e-09  # Young Modulus (kg/mm^3)
E = 207000  # Young Modulus (N/mm^2)
nu = 0.33  # Poisson ratio
clambda = E * nu / ((1 + nu) * (1 - 2 * nu))  # First Lame coefficient (N/mm^2)
cmu = E / (2 * (1 + nu))  # Second Lame coefficient (N/mm^2)

###############################################################################
# Set the numerical parameter of each case.
#

elements_degrees = []
msh_file_names = []

# Linear-Tet 2.0mm
elements_degrees.append(1)
msh_file_names.append("./LIN-TET/TET-2-0.msh")

# Linear-Tet 1.0mm
elements_degrees.append(1)
msh_file_names.append("./LIN-TET/TET-1-0.msh")

# Linear-Tet 0.5mm
elements_degrees.append(1)
msh_file_names.append("./LIN-TET/TET-0-5.msh")

# Quad-Tet 2.0mm
elements_degrees.append(2)
msh_file_names.append("./LIN-TET/TET-2-0.msh")

# Quad-Tet 1.0mm
elements_degrees.append(2)
msh_file_names.append("./LIN-TET/TET-1-0.msh")

# Quad-Tet 0.5mm
elements_degrees.append(2)
msh_file_names.append("./LIN-TET/TET-0-5.msh")

# Linear-Hex 2.0mm
elements_degrees.append(1)
msh_file_names.append("./LIN-HEX/HEX-2-0.msh")

# Linear-Hex 1.0mm
elements_degrees.append(1)
msh_file_names.append("./LIN-HEX/HEX-1-0.msh")

# Linear-Hex 0.5mm
elements_degrees.append(1)
msh_file_names.append("./LIN-HEX/HEX-0-5.msh")

# Quad-Hex 2.0mm
elements_degrees.append(2)
msh_file_names.append("./LIN-HEX/HEX-2-0.msh")

# Quad-Hex 1.0mm
elements_degrees.append(2)
msh_file_names.append("./LIN-HEX/HEX-1-0.msh")

# Quad-Hex 0.5mm
elements_degrees.append(2)
msh_file_names.append("./LIN-HEX/HEX-0-5.msh")

t = time.process_time()

###############################################################################
# Importing the mesh from gmsh format.
#

meshs = []

for i, msh_file_name in enumerate(msh_file_names):
    mesh = gf.Mesh("import", "gmsh", msh_file_name)
    meshs.append(mesh)

print("Time for import mesh", time.process_time() - t)
t = time.process_time()

###############################################################################
# Definition of finite elements methods and integration method
#

mfus = []
mfds = []
mims = []

for elements_degree, mesh in zip(elements_degrees, meshs):

    mfu = gf.MeshFem(mesh, 3)
    mfu.set_classical_fem(elements_degree)
    mfus.append(mfu)

    mfd = gf.MeshFem(mesh, 1)
    mfd.set_classical_fem(elements_degree)
    mfds.append(mfd)

    mim = gf.MeshIm(mesh, elements_degree * 2)
    mims.append(mim)

###############################################################################
# We get the mass and stiffness matrices using asm function.
#

mass_matrixs = []
linear_elasticitys = []

for i, (mfu, mfd, mim) in enumerate(zip(mfus, mfds, mims)):

    mass_matrix = gf.asm_mass_matrix(mim, mfu, mfu)
    mass_matrix.scale(rho)
    mass_matrixs.append(mass_matrix)

    lambda_d = np.repeat(clambda, mfd.nbdof())
    mu_d = np.repeat(cmu, mfd.nbdof())
    linear_elasticity = gf.asm_linear_elasticity(mim, mfu, mfd, lambda_d, mu_d)
    linear_elasticitys.append(linear_elasticity)

    mass_matrix.save("hb", "M" + "{:02}".format(i) + ".hb")
    linear_elasticity.save("hb", "K" + "{:02}".format(i) + ".hb")

print("Time for assemble matrix", time.process_time() - t)
t = time.process_time()

###############################################################################
# Solve the eigenproblem.
#
# Compute the residual error for SciPy.
#
# :math:`Err=\frac{||(K-\lambda.M).\phi||_2}{||K.\phi||_2}`
#
# Convert Lambda values to Frequency values:
# :math:`freq = \frac{\sqrt(\lambda)}{2.\pi}`
#

for i, (mass_matrix, linear_elasticity, mfu) in enumerate(
    zip(mass_matrixs, linear_elasticitys, mfus)
):

    M = io.hb_read("M" + "{:02}".format(i) + ".hb")
    K = io.hb_read("K" + "{:02}".format(i) + ".hb")
    vals, vecs = linalg.eigsh(A=K, k=6, M=M, sigma=400.0, which="LA")

    omegas = np.sqrt(np.abs(vals))
    freqs = omegas / (2.0 * np.pi)

    nev = 6

    scipy_acc = np.zeros(nev)

    print(f"case{i}")

    for j in range(nev):
        lam = vals[j]  # j-th eigenvalue
        phi = vecs.T[j]  # j-th eigenshape

        kphi = K * phi.T  # K.Phi
        mphi = M * phi.T  # M.Phi

        kphi_nrm = np.linalg.norm(kphi, 2)  # Normalization scalar value

        mphi *= lam  # (K-\lambda.M).Phi
        kphi -= mphi

        scipy_acc[j] = np.linalg.norm(kphi, 2) / kphi_nrm  # compute the residual
        print(f"[{j}] : Freq = {freqs[j]:8.2f} Hz\t Residual = {scipy_acc[j]:.5}")

    np.savetxt("freqs" + "{:0=3}".format(i) + ".txt", freqs)
    mfu.export_to_vtk(
        "vecs" + "{:0=3}".format(i) + ".vtk",
        mfu,
        vecs[:, 0],
        "Mode_1",
        mfu,
        vecs[:, 1],
        "Mode_2",
        mfu,
        vecs[:, 2],
        "Mode_3",
        mfu,
        vecs[:, 3],
        "Mode_4",
        mfu,
        vecs[:, 4],
        "Mode_5",
        mfu,
        vecs[:, 5],
        "Mode_6",
    )

print("Time for solve eigen value problem", time.process_time() - t)
t = time.process_time()

CoFEA Initiative

テスト済み有限要素コード¶

CalculiX¶

Code_Aster¶

Elmer¶

MoFEM¶

GetFEM¶

CalculiX ¶

Code_Aster ¶

Elmer ¶

MoFEM ¶

GetFEM ¶