<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Stress >


Previous pageReturn to chapter overviewNext page



 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"



  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)



   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 }



   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)



  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



  for cycle=5

    elevation(Uy) from(0,0) to (L,0) range=(-amplitude,amplitude)



  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"



  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"