{ VIBRATE.PDE
******************************************************
This example shows the use of FlexPDE in transient Stress problems.
The equations of Stress/Strain in a material medium can be given as
dx(Sx) + dy(Txy) + Fx = 0
dx(Txy) + dy(Sy) + Fy = 0
where Sx and Sy are the stresses in the x- and y- directions,
Txy is the shear stress, and Fx and Fy are the body forces in the
x- and y- directions.
In a time-dependent problem, the material acceleration and viscous force
act as body forces, and are included in a new body force term
Fx1 = Fx0 - rho*dtt(Ux) + mu*del2(dt(Ux))
Fy1 = Fy0 - rho*dtt(Uy) + mu*del2(dt(Uy))
where rho is the material mass density, mu is the viscosity, and Ux and Uy
are the material displacements.
The second time derivative in the acceleration term cannot be modelled
directly in FlexPDE, but the problem can still be solved.
Define Vx and Vy as the velocities in the x and y directions; then
Vx = dt(Ux)
and Vy = dt(Uy)
The body forces are then
Fx1 = Fx0 - rho*dt(Vx) + mu*del2(Vx)
Fy1 = Fy0 - rho*dt(Vy) + mu*del2(Vy)
This results in a set of four equations in Ux,Uy,Vx and Vy.
Notice that the stress-balance equation is the Velocity equation, and it
is to this equation that boundary loads must be applied.
In the problem considered here, we have an aluminum bar one meter long and
5 cm thick suspended on the left, and driven on the right by an oscillatory
load. The load frequency is chosen to be near the resonant frequency of
the bar.
The problem models four cycles of the driving load, and runs in less about a half an hour
on a Pentium computer.
***************************************** }
title "Transient Stress analysis"
select
deltat=1.0e-7 { Start out at careful timestep, it will grow. }
ngrid=21 { Grid a little more densely than default }
regrid = off { Cell splitting causes instantaneous changes in the
effective material properties. These changes act
like small earthquakes in the material, and propagate
high-frequency noise. To avoid these effects, we
supress grid refinement. }
variables { Recall that the declared variable range, if too large,
will affect the interpretation of error, and thus the
timestep and solution accuracy }
Ux (threshold=1e-7) { Displacements }
Uy (threshold=1e-7)
Vx (threshold=1e-5) { Velocities }
Vy (threshold=1e-5)
definitions
L = 1 { the bar length, in Meters }
hL = L/2
W = 0.05 { the bar thickness, in Meters }
hW = W/2
eps = 0.01*L
nu = 0.3 { Poisson's Ratio }
E = 6.7e+10 { Young's Modulus for Aluminum (N/M^2) }
{ plane strain coefficients }
E1 = E/((1+nu)*(1-2*nu))
C11 = E1*(1-nu)
C12 = E1*nu
C22 = E1*(1-nu)
C33 = E1*(1-2*nu)/2
rho = 2700 { Kg/M^3 }
mu = 1e3 { Estimated viscosity Kg/M/sec }
smooth = 1 { artificial diffusion to smooth results (M^2/sec) }
cvel = sqrt(E/rho) { sound velocity, M/sec }
tau = L/cvel { transit time }
tone = 0.25/tau { approximate resonant frequency }
freq = 1.1*tone { driving frequency }
period = 1/freq
amplitude=1e-8 { a guess for plot scaling }
mag=1/amplitude
force = 25 { loading force in Newtons (~1 pound force) }
{ distribute the force uniformly over the driven end: }
fdist = force/W
{ the driving force is sinusoidal in time: }
jiggle = force*sin(2*pi*freq*t)
Sx = [C11*dx(Ux) + C12*dy(Uy)] { Stresses }
Sy = [C12*dx(Ux) + C22*dy(Uy)]
Txy = C33*[dy(Ux) + dx(Uy)]
initial values
Ux = 0 { start at rest }
Uy = 0
Vx = 0
Vy = 0
equations { define the displacement equations }
Ux: Vx + smooth*div(grad(Ux)) = dt(Ux)
Uy: Vy + smooth*div(grad(Uy)) = dt(Uy)
Vx: dx(Sx) + dy(Txy) + mu*div(grad(Vx)) = rho*dt(Vx)
Vy: dx(Txy) + dy(Sy) + mu*div(grad(Vy)) = rho*dt(Vy)
boundaries
region 1
start (0,-hW)
load(Vx)=0 { free boundary on bottom, no normal stress }
load(Vy)=0
line to (L,-hW)
load(Vx) = 0
load(Vy) = jiggle { Apply oscillatory vertical load on end.
Note that this driving force must be applied to the
equation which contains the stress divergence. }
line to (L,hW)
load(Vx)=0 { free boundary on top, no normal stress }
load(Vy)=0
line to (0,hW)
value(Ux) = 0 { freeze left end (both displacement and velocity) }
value(Uy) = 0
value(Vx) = 0
value(Vy) = 0
line to close
feature { a "Gridding Feature" to force grid refinement near
the mount }
start (hW/2,-hW) line to (hW/2,hW)
start (L-hW/2,-hW) line to (L-hW/2,hW)
time 0 to 4*period
monitors
for cycle=5
elevation(Uy) from(0,0) to (L,0) range=(-amplitude,amplitude)
plots
for t= period/2 by period/2 to endtime
grid(x+mag*Ux,y+mag*Uy) as "deformation" { show final deformed grid }
vector(Ux,Uy) as "displacement" { show displacement field }
vector(Vx,Vy) as "velocity" { show velocity field }
contour(Ux) as "X-displacement(M)"
contour(Uy) as "Y-displacement(M)"
contour(Vx) as "X-velocity(M/s)"
contour(Vy) as "Y-velocity(M/s)"
contour(Sx) as "X-Stress"
contour(Sy) as "Y-Stress"
contour(Txy) as "Shear Stress"
histories
history(Ux) at (L,0) (0.8*L,0) (hL,0) as "Horizontal Displacement(M)"
history(Vx) at (L,0) (0.8*L,0) (hL,0) as "Horizontal Velocity(M/s)"
history(Sx) at (eps,hW-eps) (eps,-hW+eps) (L-eps,hW-eps) (L-eps,-hW+eps)
as "Horizontal Stress"
history(Uy) at (L,0) (0.8*L,0) (hL,0) as "Vertical Displacement(M)"
history(Vy) at (L,0) (0.8*L,0) (hL,0) as "Vertical Velocity(M/s)"
history(Sy) at (eps,hW-eps) (eps,-hW+eps) (L-eps,hW-eps) (L-eps,-hW+eps)
as "Vertical Stress"
history(Txy) at (eps,hW-eps) (eps,-hW+eps) (L-eps,hW-eps) (L-eps,-hW+eps)
as "Shear Stress"
end