Analysing mechanical stresses in Haskell.

StressAnalysis.hs

-- Define a module with the same name as the file.
module StressAnalysis where

-- A Beam is defined by its geometry. Think of this as a class declaration.
data Beam = Circle {
    diameter :: Float -- Member named 'diameter' of type 'Float'
 } deriving (Show)

-- Define functions for common properties of beams. The functions are
-- impleented down the bottom of the file. These functions take a Beam and
-- convert it into a Float.
beamI :: Beam -> Float
beamJ :: Beam -> Float
-- These two functions have two 'arguments' - that is, they take a Beam, and
-- return a function that takes a Float and returns a Float. In other words, if
-- you give them a Beam and a Float, they give you a Float.
beamT :: Beam -> Float -> Float
beamQ :: Beam -> Float -> Float

-- A Section through the beam with forces acting on it.
data Section = Section {
    name :: String,   -- Identifier for this Section
    shearY  :: Float, -- Shear forces acting on the beam
    shearZ  :: Float,
    momentX :: Float, -- Bending moments acting on the beam
    momentY :: Float,
    momentZ :: Float
 } deriving (Eq, Ord)

-- Custom show routine for Sections.
instance Show Section where
    show = name

-- Constant Section used like a constructor.
newSection = Section {
    name = "",
    shearY = 0,
    shearZ = 0,
    momentX = 0,
    momentY = 0,
    momentZ = 0
 }

-- An Element is a small finite area on a Section.
data Element = Element {
    distance :: Float, -- From centre of shaft
    angle    :: Float  -- CCW from horizontal
 } deriving (Eq, Ord)

instance Show Element where
    show e = "distance: " ++ (show $ distance e) ++ ", angle: " ++ (show $ angle e)

newElement = Element {
    distance = 0,
    angle = 0
 }

-- Calculate principal stress of a condition on a beam using Mohr's circle.
principalStress :: Beam -> Section -> Element -> Float
principalStress b s e =
    -- Mohr's circle
    1/2 * sigma + sqrt (1/4 * sigma^2 + tau^2)
    where
        sigma = bendingStress b s e
        tau   = shearStress b s e

-- Calculate the bending stress on a beam given a section and an element on that
-- section. Calculates sum of both moments given in the Section.
bendingStress :: Beam -> Section -> Element -> Float
bendingStress b s e =
    -- sigma = Mc / I
    ((momentY s) * zDist + (momentZ s) * yDist) / (beamI b)
    where
        -- Distance of this element from the neutral axes
        yDist = (distance e) * cos (angle e)
        zDist = (distance e) * sin (angle e)

-- Calculate shear stress on a beam given a section and element. Takes torque
-- into account.
shearStress :: Beam -> Section -> Element -> Float
shearStress b s e =
    -- Sum X and Y components as vectors
    sqrt (tauY^2 + tauZ^2)
    where
        -- tau = VQ / It
        tauY
            | tY == 0   = 0
            | otherwise = (shearY s) * (beamQ b hY) / (beamI b) / tY
                        + torsion * sin (angle e)
        tauZ
            | tZ == 0   = 0
            | otherwise = (shearZ s) * (beamQ b hZ) / (beamI b) / tZ
                        + torsion * cos (angle e)
        -- Thicknesses
        tY = (beamT b hY)
        tZ = (beamT b hZ)
        -- Distances from neutral axes
        hY = (distance e) * cos (angle e)
        hZ = (distance e) * sin (angle e)
        -- tau = Tr / J
        torsion = (momentX s) * (distance e) / (beamJ b)

-- These functions operate on Circles with diameter d. To account for beams of
-- different geometries, add more overloads with different pattern matches.
beamI (Circle d) = 1/4 * (d/2)^4
beamJ (Circle d) = 1/2 * (d/2)^4

beamT (Circle d) h
    -- Guard against trying to calculate thickness outside the cross section
    | (abs h) >= r = 0
    | otherwise    = 2 * sqrt (r^2 - h^2)
    where
        -- Define r that is used above
        r = d/2

beamQ (Circle d) h
    -- Guard against cases where h is outside the cross-section
    | (abs h) >= r = 0
    | otherwise    = area * centroid
    -- Now define what we mean when we used 'area' and 'centroid'
    where
        area     = r^2 / 2 * (theta - sin theta)
        centroid = (4 * r * (sin (theta / 2))^3) / (3 * theta - 3 * sin theta)
        -- Also define theta and r since we used them in the last two definitions
        theta    = 2 * acos (h / r)
        r        = d / 2