Dear all,
I am trying to apply shear wave load (such as F=sinwt) on the top boundary of a square plate.
My code is following:
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np
L = 25.
H = 1.
Nx = 250
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, “crossed”)
File(‘before0.pvd’) << mesh
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = “plane_stress”
mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
if model == “plane_stress”:
lmbda = 2mulmbda/(lmbda+2*mu)
def sigma(v):
return lmbdatr(eps(v))Identity(2) + 2.0mueps(v)
rho_g = 1e-3
f = Constant((0,-rho_g))
V = VectorFunctionSpace(mesh, ‘Lagrange’, degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0],0.)
bc = DirichletBC(V, Constant((0.,0.)), left)
u = Function(V, name=“Displacement”)
solve(a == l, u, bc)
plot(1e3u, mode=“displacement”)
print(“Maximal deflection:”, -u(L,H/2.)[1])
print(“Beam theory deflection:”, float(3rho_g*L4/2/E/H3))
Vsig = TensorFunctionSpace(mesh, “DG”, degree=0)
sig = Function(Vsig, name=“Stress”)
sig.assign(project(sigma(u), Vsig))
print(“Stress at (a,H):”, sig(0, H))
file_results = XDMFFile(“elasticity_results.xdmf”)
file_results.parameters[“flush_output”] = True
file_results.parameters[“functions_share_mesh”] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
vtkfile = File(‘displacement/solution.pvd’)
vtkfile << u
I’m not sure how to handle it properly. Can somebody help me?
Any little help is highly appreciated!