Source code for sympy.physics.continuum_mechanics.beam

"""
This module can be used to solve 2D beam bending problems with
singularity functions in mechanics.

"""

from __future__ import print_function, division

from sympy.core import S, Symbol, diff
from sympy.solvers import linsolve
from sympy.printing import sstr
from sympy.functions import SingularityFunction
from sympy.core import sympify
from sympy.integrals import integrate
from sympy.series import limit


[docs]class Beam(object): """ A Beam is a structural element that is capable of withstanding load primarily by resisting against bending. Beams are characterized by their cross sectional profile(Second moment of area), their length and their material. .. note:: While solving a beam bending problem, a user should choose its own sign convention and should stick to it. The results will automatically follow the chosen sign convention. Examples ======== There is a beam of length 4 meters. A constant distributed load of 6 N/m is applied from half of the beam till the end. There are two simple supports below the beam, one at the starting point and another at the ending point of the beam. The deflection of the beam at the end is resticted. Using the sign convention of downwards forces being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols, Piecewise >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(4, E, I) >>> b.apply_load(R1, 0, -1) >>> b.apply_load(6, 2, 0) >>> b.apply_load(R2, 4, -1) >>> b.bc_deflection = [(0, 0), (4, 0)] >>> b.boundary_conditions {'deflection': [(0, 0), (4, 0)], 'slope': []} >>> b.load R1*SingularityFunction(x, 0, -1) + R2*SingularityFunction(x, 4, -1) + 6*SingularityFunction(x, 2, 0) >>> b.solve_for_reaction_loads(R1, R2) >>> b.load -3*SingularityFunction(x, 0, -1) + 6*SingularityFunction(x, 2, 0) - 9*SingularityFunction(x, 4, -1) >>> b.shear_force() -3*SingularityFunction(x, 0, 0) + 6*SingularityFunction(x, 2, 1) - 9*SingularityFunction(x, 4, 0) >>> b.bending_moment() -3*SingularityFunction(x, 0, 1) + 3*SingularityFunction(x, 2, 2) - 9*SingularityFunction(x, 4, 1) >>> b.slope() (-3*SingularityFunction(x, 0, 2)/2 + SingularityFunction(x, 2, 3) - 9*SingularityFunction(x, 4, 2)/2 + 7)/(E*I) >>> b.deflection() (7*x - SingularityFunction(x, 0, 3)/2 + SingularityFunction(x, 2, 4)/4 - 3*SingularityFunction(x, 4, 3)/2)/(E*I) >>> b.deflection().rewrite(Piecewise) (7*x - Piecewise((x**3, x > 0), (0, True))/2 - 3*Piecewise(((x - 4)**3, x - 4 > 0), (0, True))/2 + Piecewise(((x - 2)**4, x - 2 > 0), (0, True))/4)/(E*I) """ def __init__(self, length, elastic_modulus, second_moment, variable=Symbol('x')): """Initializes the class. Parameters ========== length : Sympifyable A Symbol or value representing the Beam's length. elastic_modulus : Sympifyable A SymPy expression representing the Beam's Modulus of Elasticity. It is a measure of the stiffness of the Beam material. second_moment : Sympifyable A SymPy expression representing the Beam's Second moment of area. It is a geometrical property of an area which reflects how its points are distributed with respect to its neutral axis. variable : Symbol, optional A Symbol object that will be used as the variable along the beam while representing the load, shear, moment, slope and deflection curve. By default, it is set to ``Symbol('x')``. """ self.length = length self.elastic_modulus = elastic_modulus self.second_moment = second_moment self.variable = variable self._boundary_conditions = {'deflection': [], 'slope': []} self._load = 0 self._reaction_loads = {} def __str__(self): str_sol = 'Beam({}, {}, {})'.format(sstr(self._length), sstr(self._elastic_modulus), sstr(self._second_moment)) return str_sol @property def reaction_loads(self): """ Returns the reaction forces in a dictionary.""" return self._reaction_loads @property def length(self): """Length of the Beam.""" return self._length @length.setter def length(self, l): self._length = sympify(l) @property def variable(self): """ A symbol that can be used as a variable along the length of the beam while representing load distribution, shear force curve, bending moment, slope curve and the deflection curve. By default, it is set to ``Symbol('x')``, but this property is mutable. Examples ======== >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> x, y, z = symbols('x, y, z') >>> b = Beam(4, E, I) >>> b.variable x >>> b.variable = y >>> b.variable y >>> b = Beam(4, E, I, z) >>> b.variable z """ return self._variable @variable.setter def variable(self, v): if isinstance(v, Symbol): self._variable = v else: raise TypeError("""The variable should be a Symbol object.""") @property def elastic_modulus(self): """Young's Modulus of the Beam. """ return self._elastic_modulus @elastic_modulus.setter def elastic_modulus(self, e): self._elastic_modulus = sympify(e) @property def second_moment(self): """Second moment of area of the Beam. """ return self._second_moment @second_moment.setter def second_moment(self, i): self._second_moment = sympify(i) @property def boundary_conditions(self): """ Returns a dictionary of boundary conditions applied on the beam. The dictionary has three kewwords namely moment, slope and deflection. The value of each keyword is a list of tuple, where each tuple contains loaction and value of a boundary condition in the format (location, value). Examples ======== There is a beam of length 4 meters. The bending moment at 0 should be 4 and at 4 it should be 0. The slope of the beam should be 1 at 0. The deflection should be 2 at 0. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> b = Beam(4, E, I) >>> b.bc_deflection = [(0, 2)] >>> b.bc_slope = [(0, 1)] >>> b.boundary_conditions {'deflection': [(0, 2)], 'slope': [(0, 1)]} Here the deflection of the beam should be ``2`` at ``0``. Similarly, the slope of the beam should be ``1`` at ``0``. """ return self._boundary_conditions @property def bc_slope(self): return self._boundary_conditions['slope'] @bc_slope.setter def bc_slope(self, s_bcs): self._boundary_conditions['slope'] = s_bcs @property def bc_deflection(self): return self._boundary_conditions['deflection'] @bc_deflection.setter def bc_deflection(self, d_bcs): self._boundary_conditions['deflection'] = d_bcs
[docs] def apply_load(self, value, start, order, end=None): """ This method adds up the loads given to a particular beam object. Parameters ========== value : Sympifyable The magnitude of an applied load. start : Sympifyable The starting point of the applied load. For point moments and point forces this is the location of application. order : Integer The order of the applied load. - For moments, order= -2 - For point loads, order=-1 - For constant distributed load, order=0 - For ramp loads, order=1 - For parabolic ramp loads, order=2 - ... so on. end : Sympifyable, optional An optional argument that can be used if the load has an end point within the length of the beam. Examples ======== There is a beam of length 4 meters. A moment of magnitude 3 Nm is applied in the clockwise direction at the starting point of the beam. A pointload of magnitude 4 N is applied from the top of the beam at 2 meters from the starting point and a parabolic ramp load of magnitude 2 N/m is applied below the beam starting from 2 meters to 3 meters away from the starting point of the beam. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> b = Beam(4, E, I) >>> b.apply_load(-3, 0, -2) >>> b.apply_load(4, 2, -1) >>> b.apply_load(-2, 2, 2, end = 3) >>> b.load -3*SingularityFunction(x, 0, -2) + 4*SingularityFunction(x, 2, -1) - 2*SingularityFunction(x, 2, 2) + 2*SingularityFunction(x, 3, 0) + 2*SingularityFunction(x, 3, 2) """ x = self.variable value = sympify(value) start = sympify(start) order = sympify(order) self._load += value*SingularityFunction(x, start, order) if end: if order == 0: self._load -= value*SingularityFunction(x, end, order) elif order.is_positive: self._load -= value*SingularityFunction(x, end, order) + value*SingularityFunction(x, end, 0) else: raise ValueError("""Order of the load should be positive.""")
@property def load(self): """ Returns a Singularity Function expression which represents the load distribution curve of the Beam object. Examples ======== There is a beam of length 4 meters. A moment of magnitude 3 Nm is applied in the clockwise direction at the starting point of the beam. A pointload of magnitude 4 N is applied from the top of the beam at 2 meters from the starting point and a parabolic ramp load of magnitude 2 N/m is applied below the beam starting from 3 meters away from the starting point of the beam. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> b = Beam(4, E, I) >>> b.apply_load(-3, 0, -2) >>> b.apply_load(4, 2, -1) >>> b.apply_load(-2, 3, 2) >>> b.load -3*SingularityFunction(x, 0, -2) + 4*SingularityFunction(x, 2, -1) - 2*SingularityFunction(x, 3, 2) """ return self._load
[docs] def solve_for_reaction_loads(self, *reactions): """ Solves for the reaction forces. Examples ======== There is a beam of length 30 meters. A moment of magnitude 120 Nm is applied in the clockwise direction at the end of the beam. A pointload of magnitude 8 N is applied from the top of the beam at the starting point. There are two simple supports below the beam. One at the end and another one at a distance of 10 meters from the start. The deflection is restricted at both the supports. Using the sign convention of upward forces and clockwise moment being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols, linsolve, limit >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(30, E, I) >>> b.apply_load(-8, 0, -1) >>> b.apply_load(R1, 10, -1) # Reaction force at x = 10 >>> b.apply_load(R2, 30, -1) # Reaction force at x = 30 >>> b.apply_load(120, 30, -2) >>> b.bc_deflection = [(10, 0), (30, 0)] >>> b.load R1*SingularityFunction(x, 10, -1) + R2*SingularityFunction(x, 30, -1) - 8*SingularityFunction(x, 0, -1) + 120*SingularityFunction(x, 30, -2) >>> b.solve_for_reaction_loads(R1, R2) >>> b.reaction_loads {R1: 6, R2: 2} >>> b.load -8*SingularityFunction(x, 0, -1) + 6*SingularityFunction(x, 10, -1) + 120*SingularityFunction(x, 30, -2) + 2*SingularityFunction(x, 30, -1) """ x = self.variable l = self.length shear_curve = limit(self.shear_force(), x, l) moment_curve = limit(self.bending_moment(), x, l) reaction_values = linsolve([shear_curve, moment_curve], reactions).args self._reaction_loads = dict(zip(reactions, reaction_values[0])) self._load = self._load.subs(self._reaction_loads)
[docs] def shear_force(self): """ Returns a Singularity Function expression which represents the shear force curve of the Beam object. Examples ======== There is a beam of length 30 meters. A moment of magnitude 120 Nm is applied in the clockwise direction at the end of the beam. A pointload of magnitude 8 N is applied from the top of the beam at the starting point. There are two simple supports below the beam. One at the end and another one at a distance of 10 meters from the start. The deflection is restricted at both the supports. Using the sign convention of upward forces and clockwise moment being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(30, E, I) >>> b.apply_load(-8, 0, -1) >>> b.apply_load(R1, 10, -1) >>> b.apply_load(R2, 30, -1) >>> b.apply_load(120, 30, -2) >>> b.bc_deflection = [(10, 0), (30, 0)] >>> b.solve_for_reaction_loads(R1, R2) >>> b.shear_force() -8*SingularityFunction(x, 0, 0) + 6*SingularityFunction(x, 10, 0) + 120*SingularityFunction(x, 30, -1) + 2*SingularityFunction(x, 30, 0) """ x = self.variable return integrate(self.load, x)
[docs] def bending_moment(self): """ Returns a Singularity Function expression which represents the bending moment curve of the Beam object. Examples ======== There is a beam of length 30 meters. A moment of magnitude 120 Nm is applied in the clockwise direction at the end of the beam. A pointload of magnitude 8 N is applied from the top of the beam at the starting point. There are two simple supports below the beam. One at the end and another one at a distance of 10 meters from the start. The deflection is restricted at both the supports. Using the sign convention of upward forces and clockwise moment being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(30, E, I) >>> b.apply_load(-8, 0, -1) >>> b.apply_load(R1, 10, -1) >>> b.apply_load(R2, 30, -1) >>> b.apply_load(120, 30, -2) >>> b.bc_deflection = [(10, 0), (30, 0)] >>> b.solve_for_reaction_loads(R1, R2) >>> b.bending_moment() -8*SingularityFunction(x, 0, 1) + 6*SingularityFunction(x, 10, 1) + 120*SingularityFunction(x, 30, 0) + 2*SingularityFunction(x, 30, 1) """ x = self.variable return integrate(self.shear_force(), x)
[docs] def slope(self): """ Returns a Singularity Function expression which represents the slope the elastic curve of the Beam object. Examples ======== There is a beam of length 30 meters. A moment of magnitude 120 Nm is applied in the clockwise direction at the end of the beam. A pointload of magnitude 8 N is applied from the top of the beam at the starting point. There are two simple supports below the beam. One at the end and another one at a distance of 10 meters from the start. The deflection is restricted at both the supports. Using the sign convention of upward forces and clockwise moment being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(30, E, I) >>> b.apply_load(-8, 0, -1) >>> b.apply_load(R1, 10, -1) >>> b.apply_load(R2, 30, -1) >>> b.apply_load(120, 30, -2) >>> b.bc_deflection = [(10, 0), (30, 0)] >>> b.solve_for_reaction_loads(R1, R2) >>> b.slope() (-4*SingularityFunction(x, 0, 2) + 3*SingularityFunction(x, 10, 2) + 120*SingularityFunction(x, 30, 1) + SingularityFunction(x, 30, 2) + 4000/3)/(E*I) """ x = self.variable E = self.elastic_modulus I = self.second_moment if not self._boundary_conditions['slope']: return diff(self.deflection(), x) C3 = Symbol('C3') slope_curve = integrate(self.bending_moment(), x) + C3 bc_eqs = [] for position, value in self._boundary_conditions['slope']: eqs = slope_curve.subs(x, position) - value bc_eqs.append(eqs) constants = list(linsolve(bc_eqs, C3)) slope_curve = slope_curve.subs({C3: constants[0][0]}) return S(1)/(E*I)*slope_curve
[docs] def deflection(self): """ Returns a Singularity Function expression which represents the elastic curve or deflection of the Beam object. Examples ======== There is a beam of length 30 meters. A moment of magnitude 120 Nm is applied in the clockwise direction at the end of the beam. A pointload of magnitude 8 N is applied from the top of the beam at the starting point. There are two simple supports below the beam. One at the end and another one at a distance of 10 meters from the start. The deflection is restricted at both the supports. Using the sign convention of upward forces and clockwise moment being positive. >>> from sympy.physics.continuum_mechanics.beam import Beam >>> from sympy import symbols >>> E, I = symbols('E, I') >>> R1, R2 = symbols('R1, R2') >>> b = Beam(30, E, I) >>> b.apply_load(-8, 0, -1) >>> b.apply_load(R1, 10, -1) >>> b.apply_load(R2, 30, -1) >>> b.apply_load(120, 30, -2) >>> b.bc_deflection = [(10, 0), (30, 0)] >>> b.solve_for_reaction_loads(R1, R2) >>> b.deflection() (4000*x/3 - 4*SingularityFunction(x, 0, 3)/3 + SingularityFunction(x, 10, 3) + 60*SingularityFunction(x, 30, 2) + SingularityFunction(x, 30, 3)/3 - 12000)/(E*I) """ x = self.variable E = self.elastic_modulus I = self.second_moment if not self._boundary_conditions['deflection'] and not self._boundary_conditions['slope']: return S(1)/(E*I)*integrate(integrate(self.bending_moment(), x), x) elif not self._boundary_conditions['deflection']: return integrate(self.slope(), x) elif not self._boundary_conditions['slope'] and self._boundary_conditions['deflection']: C3 = Symbol('C3') C4 = Symbol('C4') slope_curve = integrate(self.bending_moment(), x) + C3 deflection_curve = integrate(slope_curve, x) + C4 bc_eqs = [] for position, value in self._boundary_conditions['deflection']: eqs = deflection_curve.subs(x, position) - value bc_eqs.append(eqs) constants = list(linsolve(bc_eqs, (C3, C4))) deflection_curve = deflection_curve.subs({C3: constants[0][0], C4: constants[0][1]}) return S(1)/(E*I)*deflection_curve C4 = Symbol('C4') deflection_curve = integrate((E*I)*self.slope(), x) + C4 bc_eqs = [] for position, value in self._boundary_conditions['deflection']: eqs = deflection_curve.subs(x, position) - value bc_eqs.append(eqs) constants = list(linsolve(bc_eqs, C4)) deflection_curve = deflection_curve.subs({C4: constants[0][0]}) return S(1)/(E*I)*deflection_curve