<< Click to Display Table of Contents >> 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 }  

    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 }