beamsolve

The beamsolve proc solves a (possibly multi-span) Euler beam equation:

subject to initial and boundary conditions.  The solution  is the

transverse deflection of the beam at point  at time , subject to the load
density (i.e., load per unit length) given by . The coefficient  

is the beam's mass density (mass per unit length),  is the Young's modulus of

the beam's material, and  is the beam's cross-sectional moment of inertia

about the neutral axis.  The figure below illustrates a 3-span beam (drawn in green)
supported on four supports, and loaded by a variable density load (drawn in gray)
which may vary in time.  The objective is to determine the deformed shape of the
beam as a function of time.

The number of spans, their lengths, and the nature of the supports may be specified

as arguments to beamsolve.

In this worksheet we assume that  are constants for simplicity. Since only
the product of the coefficients  and  enters the calculations, we lump the two

together into a single variable which we indicate with the two-letter symbol .

Commonly,  is referred to as the beam's rigidity.

The PDE needs to be supplied with boundary conditions, two at each end, each

condition prescribing a value (possibly time-dependent) for one of  
(that's 36 possible combinations!) where I have used subscripts to indicate

derivatives.  Thus, for a single-span beam of length , the following is an admissible

set of boundary conditions:
.   (Oops, coorection, that last
condition was meant to be u_xxx(L,t) = sin t.)
Additionally, the PDE needs to be supplied with initial conditions which express

the initial displacement and the initial velocity:

The PDE is solved through the Method of Lines.  Thus, each span is subdivided into

subintervals and the PDE's spatial derivatives are approximated through finite differences.

The time derivatives, however, are not discretized.  This reduces the PDE into a set of

ODEs which are solved with Maple's dsolve().  

Calling sequence:

        beamsolve(, options)

Parameters:

        :  List of span lengths, in order from left to right, as in .

         The number of subintervals in the shortest span (for the finite difference approximation)

Notes:

Options:

        All options are of the form option_name=value, and have built-in default values.

        Only options that are other than the defaults need be specified.

        rho: the beam's (constant) mass density per unit length (default: rho = 1)

        EI: the beam's (constant) rigidity (default: EI = 1)

        T: solve the PDE over the time interval  (default: T = 1)

        F: an expression in  and  that describes the applied force   (default: F = 0)
        IC: the list of the initial displacement and velocity,  as
                expressions in  (default: IC = [0,0])

        BC_left: a list consisting of a pair of boundary conditions at the left end of
                the overall (possibly multi-span beam.  These can be any two of
                . The right-hand sides of these equations

                can be any expression in .  The left-hand sides should be entered literally as indicated.

                If a right-hand side is omitted, it is taken to be zero.   (default: BC_left =  which

                corresponds to a simple support).

        BC_right: like BC_left, but for the right end of the overall beam (default: BC_right =

The returned module:

        A call to beamsolve returns a module which presents three methods.  The methods are:

        plot (t, refine=val, options)

                plots the solution  at time .  If the discretization in the  direction

                is too coarse and the graph looks non-smooth, the refine option

                (default: refine=1) may be specified to smooth out the graph by introducing

                val number of intermediate points within each discretized subinterval.

                All other options are assumed to be plot options and are passed to plots:-display.

        plot3d (m=val, options)

                plots the surface .  The optional  specification requests

                a grid consisting of  subintervals in the time direction (default:

                Note that this grid is for plotting purposes only; the solution is computed

                as a continuous (not discrete) function of time. All other options are assumed

                to be plot3d options and are passed to plots:-display.

        animate (frames=val, refine=val, options)

                produces an animation of the beam's motion.  The frames option (default = 50)

                specifies the number of animation frames.  The refine option is passed to plot
                (see the description above. All other options are assumed to be plot options and
                are passed to plots:-display.

Note:

        In specifying the boundary conditions, the following reminder can be helpful.  If the beam

        is considered to be horizontal, then  is the vertical displacement,  is the slope,  

        is the bending moment, and  is the transverse shear force.

A single-span simply-supported beam with initial velocity

The function  is an exact solution of a simply supported beam with

  The solution is periodic in time with period .

The initial condition   does not lead to a separable form, and

therefore the motion is more complex.

A single-span cantilever beam

A cantilever beam with initial condition  where  is the
first eigenmode of its free vibration (calculated in another spreadsheet).  The motion is
periodic in time, with period

If the initial condition is not an eigenmode, then the solution is rather chaotic.

A single-span cantilever beam with a weight hanging from its free end

A single-span cantilever beam with oscillating support

A dual-span simply-supported beam with moving load

Load moves across a dual-span beam.

The beam continues oscillating after the load leaves.

A triple-span simply-supported beam with moving load

Load moves across a triple-span beam.

The beam continues oscillating after the load leaves.

A triple-span beam, moving load falling off the cantilever end

In this demo the load move across a multi-span beam with a cantilever section at the right.

As it skips past the cantilever end, the beam snaps back violently.

rand := proc(r::{posint, numeric .. numeric}, $)
   1   if nargs = 0 then
   2     RandomTools:-MersenneTwister:-GenerateInteger()
       else
   3     try
   4       if type(r,('numeric') .. ('numeric')) then
   5         if rhs(r) < lhs(r) then
   6           error "invalid range, end points are out of order"
             elif not type(r,('integer') .. ('integer')) then
   7           return () -> RandomTools:-Generate(float(('range') = r,('method') = ('uniform')))
             end if
           end if;
   8       return RandomTools:-MersenneTwister:-NewGenerator(('range') = r)
         catch :
   9       error
         end try
       end if
end proc

RandomTools:-MersenneTwister:-GenerateInteger := proc({range::{integer, integer .. integer} := 1000000000000}, $)
local bits, divisor, offset;
   1   if nargs = 0 then
   2     return MTKernelInterface(4)
       end if;
   3   if type(range,'integer') then
   4     divisor := range;
   5     offset := 0
       else
   6     if rhs(range) < lhs(range) then
   7       error "illegal range, end points out of order"
         end if;
   8     if lhs(range) = rhs(range) then
   9       return lhs(range)
         end if;
  10     offset := lhs(range);
  11     divisor := rhs(range)-offset+1
       end if;
  12   if divisor = 0 then
  13     return offset
       end if;
  14   bits := MTKernelInterface(5,divisor);
  15   MTKernelInterface(6,divisor,bits)+offset
end proc

Error, (in showstat) procedure name expected

Example: Solving a Double Inequality

Example: Solving an Equation Involving an Absolute Value