使用FEniCS求解悬臂板固有频率结果异常求助

我正尝试用FEniCS求解简单悬臂板的固有频率，目标是匹配某论文的结果，但目前得到的结果不正确，排查许久仍未找到问题。以下是代码运行输出及完整代码：

运行输出 #

/miniconda3/envs/BASim/lib/python3.10/site-packages/ufl/__init__.py:250: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources
[DESKTOP-PVK4R9R:28050] mca_base_component_repository_open: unable to open mca_btl_openib: librdmacm.so.1: cannot open shared object file: No such file or directory (ignored)
Computing 10 first eigenvalues...

Eigenfrequencies [Hz]:
  Mode 1: 12.11532 Hz
  Mode 2: 44.17645 Hz
  Mode 3: 75.86820 Hz
  Mode 4: 212.21367 Hz
  Mode 5: 273.42441 Hz
  Mode 6: 414.87689 Hz
  Mode 7: 437.26501 Hz
  Mode 8: 685.17764 Hz
  Mode 9: 751.14498 Hz
  Mode 10: 1021.38196 Hz

完整代码 #

from dolfin import *
import numpy as np

# =============================================================================
# GEOMETRY AND MATERIAL PROPERTIES FROM THE PAPER
# =============================================================================
L = 1.0  # Beam length [m]
B = 0.050  # Beam width [m]  
H = 0.005  # Beam height [m]

E_val = 70e9  # Young's modulus [Pa] = 70 GPa
nu_val = 0.35   # Poisson's ratio
rho_val = 2700  # Density [kg/m³]

E = Constant(E_val)
nu = Constant(nu_val)
rho = Constant(rho_val)

# =============================================================================
# MESH GENERATION
# =============================================================================
Nx = 100
Ny = int(B/L*Nx)+1
Nz = int(H/L*Nx)+1

mesh = BoxMesh(Point(0., 0., 0.), Point(L, B, H), Nx, Ny, Nz)

# =============================================================================
# CONSTITUTIVE RELATIONS
# =============================================================================
mu = E / 2.0 / (1.0 + nu)
lmbda = E * nu / (1.0 + nu) / (1.0 - 2.0 * nu)

def eps(v):
    return sym(grad(v))

def sigma(v):
    dim = v.geometric_dimension()
    return 2.0 * mu * eps(v) + lmbda * tr(eps(v)) * Identity(dim)

# =============================================================================
# FUNCTION SPACE AND BOUNDARY CONDITIONS
# =============================================================================
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u_ = TrialFunction(V)
du = TestFunction(V)

def left(x, on_boundary):
    return near(x[0], 0.0) and on_boundary

bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), left)

# =============================================================================
# ASSEMBLE STIFFNESS AND MASS MATRICES
# =============================================================================
k_form = inner(sigma(du), eps(u_)) * dx
m_form = rho * dot(du, u_) * dx

K = PETScMatrix()
M = PETScMatrix()

assemble(k_form, tensor=K)
assemble(m_form, tensor=M)

bc.apply(K)
bc.zero(M)

# =============================================================================
# EIGENVALUE SOLVER
# =============================================================================
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.0

N_eig = 10
print(f"Computing {N_eig} first eigenvalues...")
eigensolver.solve(N_eig)

# =============================================================================
# PRINT EIGENFREQUENCIES
# =============================================================================
print("\nEigenfrequencies [Hz]:")
for i in range(N_eig):
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    
    if r > 0:
        freq = np.sqrt(r) / (2.0 * np.pi)
        print(f"  Mode {i+1}: {freq:.5f} Hz")
    else:
        print(f"  Mode {i+1}: negative eigenvalue ({r:.2e}) - skipped")

内容的提问来源于stack exchange，提问作者FN-2187