{ HARMONIC.PDE
******************************************************
This example shows the use of FlexPDE in harmonic analysis of
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(U) + mu*del2(dt(U))
Fy1 = Fy0 - rho*dtt(V) + mu*del2(dt(V))
where rho is the material mass density, mu is the viscosity, and U and V
are the material displacements in the x and y directions.
If we assume that the displacement is harmonic in time (all transients
have died out), then we can assert
U(t) = U0*exp(-i*omega*t)
V(t) = V0*exp(-i*omega*t)
Here U0(x,y) and V0(x,y) are the complex amplitude distributions, and
omega is the angular velocity of the oscillation.
Substituting this assumption into the stress equations and dividing out
the common exponential factors, we get (implying U0 by U and V0 by V)
dx(Sx) + dy(Txy) + Fx0 + rho*omega^2*U - i*omega*mu*del2(U) = 0
dx(Txy) + dy(Sy) + Fy0 + rho*omega^2*V - i*omega*mu*del2(V) = 0
All the terms in this equation are now complex. Separating into real
and imaginary parts gives
U = Ur + i*Ui
Sx = Srx + i*Six
Sy = Sry + i*Siy
etc...
Expressed in terms of the constitutive relations of the material,
Srx = [C11*dx(Ur) + C12*dy(Vr)]
Sry = [C12*dx(Ur) + C22*dy(Vr)]
Trxy = C33*[dy(Ur) + dx(Vr)]
etc...
The final result is a set of four equations in Ur,Vr,Ui and Vi.
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.
We run the problem in three stages, first with no viscosity, then with increasing
viscosities to show the mixing of real and imaginary components.
***************************************** }
title "Harmonic Stress analysis"
variables { Recall that the declared variable range, if too large,
will affect the interpretation of error, and thus the
timestep and solution accuracy }
Ur { Displacements }
Ui
Vr
Vi
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 = staged(0,1e3,1e4) { Estimated viscosity Kg/M/sec }
cvel = sqrt(E/rho) { sound velocity, M/sec }
tau = L/cvel { transit time }
tone = 0.25/tau { approximate resonant frequency }
omega = 2*pi*tone { driving angular velocity }
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
Um = sqrt(Ur^2+Ui^2) { X-displacement amplitude }
Vm = sqrt(Vr^2+Vi^2) { X-displacement amplitude }
Srx = [C11*dx(Ur) + C12*dy(Vr)] { Real Stresses }
Sry = [C12*dx(Ur) + C22*dy(Vr)]
Trxy = C33*[dy(Ur) + dx(Vr)]
Six = [C11*dx(Ui) + C12*dy(Vi)] { Imaginary Stresses }
Siy = [C12*dx(Ui) + C22*dy(Vi)]
Tixy = C33*[dy(Ui) + dx(Vi)]
Sxm = sqrt(Srx^2+Six^2)
Sym = sqrt(Sry^2+Siy^2)
Txym = sqrt(Trxy^2+Tixy^2)
equations { define the displacement equations }
Ur: dx(Srx) + dy(Trxy) + rho*omega^2*Ur + omega*mu*del2(Ui) = 0
Ui: dx(Six) + dy(Tixy) + rho*omega^2*Ui - omega*mu*del2(Ur) = 0
Vr: dx(Trxy) + dy(Sry) + rho*omega^2*Vr + omega*mu*del2(Vi) = 0
Vi: dx(Tixy) + dy(Siy) + rho*omega^2*Vi - omega*mu*del2(Vr) = 0
boundaries
region 1
start (0,-hW)
load(Ur)=0 { free boundary on bottom, no normal stress }
load(Ui)=0
load(Vr)=0
load(Vi)=0
line to (L,-hW)
load(Vr) = force { Apply oscillatory vertical load on end. }
line to (L,hW)
load(Vr)=0 { free boundary on top, no normal stress }
line to (0,hW)
value(Ur) = 0 { clamp the left end }
value(Ui) = 0
value(Vr) = 0
value(Vi) = 0
line to close
monitors
elevation(Vr,Vi) from(0,0) to (L,0)
report(omega) report(mu)
plots
grid(x+mag*Ur,y+mag*Vr) as "Real displacement" { show final deformed grid }
report(omega) report(mu)
grid(x+mag*Ui,y+mag*Vi) as "Imag displacement"
report(omega) report(mu)
elevation(Vr,Vi) from(0,0) to (L,0)
report(omega) report(mu)
contour(Ur) as "Real X-displacement(M)"
report(omega) report(mu)
contour(Vr) as "Real Y-displacement(M)"
report(omega) report(mu)
contour(Ui) as "Imag X-displacement(M)"
report(omega) report(mu)
contour(Vi) as "Imag Y-displacement(M)"
report(omega) report(mu)
contour(Sxm) as "X-Stress amplitude"
report(omega) report(mu)
contour(Sym) as "Y-Stress amplitude"
report(omega) report(mu)
contour(Txym) as "Shear Stress amplitude"
report(omega) report(mu)
end