vibrate

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Stress >

vibrate

Previous pageReturn to chapter overviewNext page

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

 

}  

 

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 }  

   smoother = 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))  

 

boundary conditions

  'no load' : load(Vx)=0 load(Vy)=0

  'Y load'  : load(Vx)=0 load(Vy)=jiggle

  'freeze'  : value(Ux)=0 value(Uy)=0 value(Vx)=0 value(Vy)=0

 

initial values  

   Ux = 0             { start at rest }  

   Uy = 0  

   Vx = 0  

   Vy = 0  

 

equations             { define the displacement equations }  

   Ux:  Vx + smoother*div(grad(Ux)) = dt(Ux)  

   Uy:  Vy + smoother*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)

 

    use bc 'no load'   { free boundary on bottom, no normal stress }

    line to (L,-hW)

 

    use bc 'Y load'     { 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)

 

    use bc 'no load'   { free boundary on top, no normal stress }

    line to (0,hW)

 

    use bc 'freeze'     { freeze left end (both displacement and velocity) }

    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