﻿ Sample Problems > Applications > Stress > harmonic

# harmonic

Navigation:  Sample Problems > Applications > Stress >

# harmonic

{ 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 }

{ Displacements }

Ur

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 }

line to (L,-hW)

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