Eigenvalue of Euler-Bernoulli Beam Equation with axial load

Hi,
I hope all are fine. I need help to solve the eigenvalue problem to compute critical load for the beam using Euler-Bernoulli beam equation and apply axial load. For this I solved 1D poisson equation by using Arnoldi method my computed value is similar to the analytical value. I also want to compute eigenvector of 1D poisson equation. Can anyone tell me how I can compute the eigenvector? Its only for benchmark. But I want to compute critical load for the beam.I used mixed saddle point method reduce power of derivative. I am using the same Arnoldi method but computed value did not matched with analytical value. Please tell me about my problem. I will be very thankful to you.
Best,
Rauf.

#Code for poisson equation:
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.meshes import Make1DMesh
import math
mesh = Make1DMesh(100)

fes = H1(mesh, order=2, dirichlet=“left|right”)
u = fes.TrialFunction()
v = fes.TestFunction()

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx

m = BilinearForm(fes)
m += uvdx

a.Assemble()
m.Assemble()

u = GridFunction(fes,multidim=3)

with TaskManager():
lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),list(u.vecs), shift=1)
print ("lam: ", lam)

#Code for Beam

from ngsolve import *
from ngsolve.meshes import Make1DMesh
from ngsolve.la import EigenValues_Preconditioner
from ngsolve.webgui import Draw
import numpy as np
import math
import scipy.linalg
from scipy import random
import scipy.sparse as sp

mesh = Make1DMesh(100)

w = H1(mesh, order=2, dirichlet=“left|right”)
sigma = H1(mesh, order=2, dirichlet=" ")
fesm = w*sigma

w, sigma = fesm.TrialFunction()
v , tau = fesm.TestFunction()

a = BilinearForm(fesm,symmetric=True)
a += (grad(sigma)grad(v) + sigmatau + grad(w)*grad(tau))*dx

m = BilinearForm(fesm,symmetric=True)
m += grad(w)*grad(v)*dx

a.Assemble()
m.Assemble()

u = GridFunction(fesm)

with TaskManager():
lam = ArnoldiSolver(a.mat, m.mat, fesm.FreeDofs(),list(u.vecs), shift=1)

computed value is -262.11 but analytical value is 4pi^2

Hi Rauf,

you have to give the dimension of the Krylov-space by the number of vectors. With

u = GridFunction(fesm, multidim=10)

I get
(-39.4784,-1.99848e-16)
(-80.7629,-1.5946e-15)
(-157.914,2.09357e-16)

where the first eigenvalue is very close to (the negative of) 4 \pi^2

Joachim

here an alternative solver using the LOBPCG solver for generalized symmetric and positive definite EVPs.

bernoulli_evp2.py (635 Bytes)

in case you want to continue with mixed methods for plates you can try out the HHJ method:

HHJ - method

Hi Joachim,
Bundle of thanks for your reply. I got the same results by giving the dimension of the Krylov-space by the number of vectors. Can I compute the eigenvector using the Arnoldi method?

Best,
Rauf.

Hi Rauf,

the eigenvectors should be saved in the multidim GridFunction u.
When drawing you can select the multidim index, corresponding to the computed eigenvector.

Best,
Michael

Hi Michael!
Thanks for your reply.
Best,
Rauf.