Frequency domain analysis with damping

1

Hi,

Am currently performing frequency domain analysis of the damped structure. I observe that the real and the imaginary parts of the results obtained are the same for both the x and y components of the displacement vector. It would be really appreciable if anyone knows what exactly is happening with the below implementation. I will be looking forward to your suggestions:

from dolfin import *
import numpy
from ufl import indices
import matplotlib.pyplot as plt

frequency=50000
omega=2 * numpy.pi * frequency

L = 0.2
H = 0.01
Nx = 200
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)

#Function space definition over the mesh
Element1 = VectorElement(“CG”, mesh.ufl_cell(), 1)
Element2 = VectorElement(“CG”, mesh.ufl_cell(), 1)

Define the mixed element

V_elem = MixedElement([Element1, Element2])

Define the mixed function space

Vcomplex = FunctionSpace(mesh, V_elem)
#Material properties
E, rho, nu = 300E+8, 8000, 0.3
#Rayleigh Damping Coefficients
alpha = Constant(0)
beta = Constant(1e-6)

Lame coefficient for constitutive relation

mu = E/2./(1+nu)
lmbda = E * nu/(1+nu)/(1-2*nu)

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

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

#Applying the Boundary conditions
zero=Constant((0.0, 0.0))
one=Constant((0.0, -1.0))

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

def rigth(x, on_boundary):
return near(x[0], 0.2)

boundary = [DirichletBC(Vcomplex.sub(0), zero, left), DirichletBC(Vcomplex.sub(1), zero, left)]
source = [DirichletBC(Vcomplex.sub(0), one, rigth), DirichletBC(Vcomplex.sub(1), one, rigth)]

bc=boundary+source

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

M = omega** 2 * rho * (dot(uR,wR) + dot(uI,wI)) *dx
K = (inner(sigma(uR), grad(wR))+inner(sigma(uI), grad(wI)))*dx
C = alpha * M + beta * K

F = M - K + C
a, L = lhs(F), rhs(F)

Set up the PDE

solnU = Function(Vcomplex)
problem = LinearVariationalProblem(a, L, solnU, bcs=bc)
solver = LinearVariationalSolver(problem)
solver.solve()
solnUR, solnUI=solnU.split()

Set up file for exporting results

file_results = XDMFFile(“Amplitude_response.xdmf”)
file_results.parameters[“flush_output”] = True
file_results.parameters[“functions_share_mesh”] = True
file_results.write(solnUR, 0)
file_results.write(solnUI, 0)

2

At a glance it looks like you’re not computing the inner product correctly.

Recall that:

(u, v) = \int_\Omega u \cdot \overline{v} \; \mathrm{d}\vec{x} = \int_\Omega (u_r + j u_i)(v_r - j v_i) \; \mathrm{d} \vec{x} = \int_\Omega \left(u_r v_r + u_i v_i + j u_i v_r - j u_r v_i\right) \; \mathrm{d}x

where \cdot_r and \cdot_i denote the real and imaginary components of their argument, respectively. You’ll notice you have a system like

\begin{pmatrix} (u_r, v_r) & j(u_i, v_r) \\ -j(u_r, v_i) & (u_i, v_i) \end{pmatrix} = \begin{pmatrix} (u_r, v_r) & 0 \\ 0 & (u_i, v_i) \end{pmatrix} + j \begin{pmatrix} 0 & (u_i, v_r) \\ -(u_r, v_i) & 0 \end{pmatrix}

which is the notorious sesquilinear formulation that you’re probably looking for.

3

Hi,

Thanks for the information. I have modified the equations accordingly, but still, the results are not correct. I have made the following changes to the above code:

#Real and Imaginary parts of the trial and test functions
uR, uI = TrialFunctions(Vcomplex)
wR, wI = TestFunctions(Vcomplex)

M = omega**2*rho * (dot(uR,wR) + dot(uI,wI) + dot(uI,wR) - dot(uR,wI)) * dx
K = (inner(sigma(uR), grad(wR)) + inner(sigma(uI), grad(wI)) + inner(sigma(uI), grad(wR)) - inner(sigma(uR), grad(wI)))*dx
C = alpha * M + beta * K

F = M - K + C

a, L = lhs(F), rhs(F)

can you kindly go through the implementation?