Sample Problems > Applications > Stress

{ 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

vibrate

vibrate