from future import print_function
from dolfin import *
-----------------------------------------------------------------------------------------------
External load
class Traction(Expression):
#def init(self, end, **kwargs):
#Expression.__init__(self)
#self.t = 0.0
#self.end = end
#def eval(self, values, x):
#values[0] = 0.0
#values[1] = 0.0
#if x[0] > 0.0 - DOLFIN_EPS:
#values[0] = self.t/self.end if self.t < self.end else 1.0
#def value_shape(self):
#return (2,)
def update(u, u0, v0, a0, beta, gamma, dt):
# Get vectors (references)
u_vec, u0_vec = u.vector(), u0.vector()
v0_vec, a0_vec = v0.vector(), a0.vector()
# Update acceleration and velocity
a_vec = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/(0.5*dt*dt) -(1.0-2.0*beta)*a0_vec )
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec
# Update (t(n) <-- t(n+1))
v0.vector()[:], a0.vector()[:] = v_vec, a_vec
u0.vector()[:] = u.vector()
-----------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------
Mesh and Material definition
-----------------------------------------------------------------------------------------------
Load mesh and define function space
mesh = UnitSquareMesh(32, 32)
import os
os.system ('dolfin-convert GMSH/plaque_1m.msh GMSH/plaque_1m.xml')
mesh = Mesh("GMSH/plaque_1m.xml")
Define function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)
Test and trial functions
u1, w = TrialFunction(V), TestFunction(V)
Proprietes Aluminium
E = 71*10**9
nu = 0.346
Mass density and viscous damping coefficient
rho = 2700
eta = 0
mu, lmbda = E/(2.0(1.0 + nu)), Enu/((1.0 + nu)(1.0 - 2.0nu))
cl = sqrt((lmbda + 2*mu)/rho)
ct = sqrt(mu/rho)
print('Longitudinal waves speed',cl)
print('Transversal waves speed',ct)
-----------------------------------------------------------------------------------------------
Time integration
-----------------------------------------------------------------------------------------------
Time stepping parameters
beta, gamma = 0.25, 0.5
dt = 1*10**(-6)
t, T = 0.0, 200*dt
Fields from previous time step (displacement, velocity, acceleration)
u0, v0, a0 = Function(V), Function(V), Function(V)
Exitation terms
fcb = 9000103 # 900 KHz
a1b = -(pi*fcb)2
tdb = 110**(-6)
h = Constant ((0,0))
Velocity and acceleration at t_(n+1)
v1 = (gamma/(betadt))(u1 - u0) - (gamma/beta - 1.0)v0 - dt(gamma/(2.0beta)- 1.0)a0
a1 = (1.0/(beta*dt**2))*(u1 - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0
-----------------------------------------------------------------------------------------------
Variational formulation
-----------------------------------------------------------------------------------------------
Stress tensor
def sigma(u):
return 2.0musym(grad(u)) + lmbdatr(sym(grad(u)))Identity(len(u))
Governing equation
F = rhodot(a1, w)dx + inner(sigma(u1), sym(grad(w)))dx - dot(h,w)ds
Extract bilinear and linear forms
a = lhs(F)
L = rhs(F)
Set up boundary condition at left end
zero = Constant((0.0, 0.0))
def left(x):
return x[0] < DOLFIN_EPS
bc = DirichletBC(V, zero, left)
A, b = assemble_system(a, L, bc)
-----------------------------------------------------------------------------------------------
Solver configuration
-----------------------------------------------------------------------------------------------
Set up PDE, advance in time and solve
u = Function(V)
solve(a == L, u, bc)
problem = LinearVariationalProblem(a, L, u,bc)
solver = LinearVariationalSolver(problem)
Save solution in VTK format
file = File("Results/displacement.pvd")
while t <= T:
A, b = assemble_system(a, L, bc)
if t<0.0001:
#Exitation = PointSource (V,Point(0,0),0.5 + a1b*((t-tdb)2)exp(a1b((t-tdb)2))
#Exitation.apply(b)
delta = PointSource(V, Point(0,0), sin(2000 * t)),
delta.apply(b)
#h.t = t
#solver.solve()
solve(A, u.vector(), b)
update(u, u0, v0, a0, beta, gamma, dt)
t+=dt
file << u
