1026 lines
38 KiB
Python
1026 lines
38 KiB
Python
|
"""
|
||
|
This module can be used to solve probelsm related to 2D parabolic arches
|
||
|
"""
|
||
|
from sympy.core.sympify import sympify
|
||
|
from sympy.core.symbol import Symbol,symbols
|
||
|
from sympy import diff, sqrt, cos , sin, atan, rad, Min
|
||
|
from sympy.core.relational import Eq
|
||
|
from sympy.solvers.solvers import solve
|
||
|
from sympy.functions import Piecewise
|
||
|
from sympy.plotting import plot
|
||
|
from sympy import limit
|
||
|
from sympy.utilities.decorator import doctest_depends_on
|
||
|
from sympy.external.importtools import import_module
|
||
|
|
||
|
numpy = import_module('numpy', import_kwargs={'fromlist':['arange']})
|
||
|
|
||
|
class Arch:
|
||
|
"""
|
||
|
This class is used to solve problems related to a three hinged arch(determinate) structure.\n
|
||
|
An arch is a curved vertical structure spanning an open space underneath it.\n
|
||
|
Arches can be used to reduce the bending moments in long-span structures.\n
|
||
|
|
||
|
Arches are used in structural engineering(over windows, door and even bridges)\n
|
||
|
because they can support a very large mass placed on top of them.
|
||
|
|
||
|
Example
|
||
|
========
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.get_shape_eqn
|
||
|
5 - (x - 5)**2/5
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,1),crown_x=6)
|
||
|
>>> a.get_shape_eqn
|
||
|
9/5 - (x - 6)**2/20
|
||
|
"""
|
||
|
def __init__(self,left_support,right_support,**kwargs):
|
||
|
self._shape_eqn = None
|
||
|
self._left_support = (sympify(left_support[0]),sympify(left_support[1]))
|
||
|
self._right_support = (sympify(right_support[0]),sympify(right_support[1]))
|
||
|
self._crown_x = None
|
||
|
self._crown_y = None
|
||
|
if 'crown_x' in kwargs:
|
||
|
self._crown_x = sympify(kwargs['crown_x'])
|
||
|
if 'crown_y' in kwargs:
|
||
|
self._crown_y = sympify(kwargs['crown_y'])
|
||
|
self._shape_eqn = self.get_shape_eqn
|
||
|
self._conc_loads = {}
|
||
|
self._distributed_loads = {}
|
||
|
self._loads = {'concentrated': self._conc_loads, 'distributed':self._distributed_loads}
|
||
|
self._loads_applied = {}
|
||
|
self._supports = {'left':'hinge', 'right':'hinge'}
|
||
|
self._member = None
|
||
|
self._member_force = None
|
||
|
self._reaction_force = {Symbol('R_A_x'):0, Symbol('R_A_y'):0, Symbol('R_B_x'):0, Symbol('R_B_y'):0}
|
||
|
self._points_disc_x = set()
|
||
|
self._points_disc_y = set()
|
||
|
self._moment_x = {}
|
||
|
self._moment_y = {}
|
||
|
self._load_x = {}
|
||
|
self._load_y = {}
|
||
|
self._moment_x_func = Piecewise((0,True))
|
||
|
self._moment_y_func = Piecewise((0,True))
|
||
|
self._load_x_func = Piecewise((0,True))
|
||
|
self._load_y_func = Piecewise((0,True))
|
||
|
self._bending_moment = None
|
||
|
self._shear_force = None
|
||
|
self._axial_force = None
|
||
|
# self._crown = (sympify(crown[0]),sympify(crown[1]))
|
||
|
|
||
|
@property
|
||
|
def get_shape_eqn(self):
|
||
|
"returns the equation of the shape of arch developed"
|
||
|
if self._shape_eqn:
|
||
|
return self._shape_eqn
|
||
|
|
||
|
x,y,c = symbols('x y c')
|
||
|
a = Symbol('a',positive=False)
|
||
|
if self._crown_x and self._crown_y:
|
||
|
x0 = self._crown_x
|
||
|
y0 = self._crown_y
|
||
|
parabola_eqn = a*(x-x0)**2 + y0 - y
|
||
|
eq1 = parabola_eqn.subs({x:self._left_support[0], y:self._left_support[1]})
|
||
|
solution = solve((eq1),(a))
|
||
|
parabola_eqn = solution[0]*(x-x0)**2 + y0
|
||
|
if(parabola_eqn.subs({x:self._right_support[0]}) != self._right_support[1]):
|
||
|
raise ValueError("provided coordinates of crown and supports are not consistent with parabolic arch")
|
||
|
|
||
|
elif self._crown_x:
|
||
|
x0 = self._crown_x
|
||
|
parabola_eqn = a*(x-x0)**2 + c - y
|
||
|
eq1 = parabola_eqn.subs({x:self._left_support[0], y:self._left_support[1]})
|
||
|
eq2 = parabola_eqn.subs({x:self._right_support[0], y:self._right_support[1]})
|
||
|
solution = solve((eq1,eq2),(a,c))
|
||
|
if len(solution) <2 or solution[a] == 0:
|
||
|
raise ValueError("parabolic arch cannot be constructed with the provided coordinates, try providing crown_y")
|
||
|
parabola_eqn = solution[a]*(x-x0)**2+ solution[c]
|
||
|
self._crown_y = solution[c]
|
||
|
|
||
|
else:
|
||
|
raise KeyError("please provide crown_x to construct arch")
|
||
|
|
||
|
return parabola_eqn
|
||
|
|
||
|
@property
|
||
|
def get_loads(self):
|
||
|
"""
|
||
|
return the position of the applied load and angle (for concentrated loads)
|
||
|
"""
|
||
|
return self._loads
|
||
|
|
||
|
@property
|
||
|
def supports(self):
|
||
|
"""
|
||
|
Returns the type of support
|
||
|
"""
|
||
|
return self._supports
|
||
|
|
||
|
@property
|
||
|
def left_support(self):
|
||
|
"""
|
||
|
Returns the position of the left support.
|
||
|
"""
|
||
|
return self._left_support
|
||
|
|
||
|
@property
|
||
|
def right_support(self):
|
||
|
"""
|
||
|
Returns the position of the right support.
|
||
|
"""
|
||
|
return self._right_support
|
||
|
|
||
|
@property
|
||
|
def reaction_force(self):
|
||
|
"""
|
||
|
return the reaction forces generated
|
||
|
"""
|
||
|
return self._reaction_force
|
||
|
|
||
|
def apply_load(self,order,label,start,mag,end=None,angle=None):
|
||
|
"""
|
||
|
This method adds load to the Arch.
|
||
|
|
||
|
Parameters
|
||
|
==========
|
||
|
|
||
|
order : Integer
|
||
|
Order of the applied load.
|
||
|
|
||
|
- For point/concentrated loads, order = -1
|
||
|
- For distributed load, order = 0
|
||
|
|
||
|
label : String or Symbol
|
||
|
The label of the load
|
||
|
- should not use 'A' or 'B' as it is used for supports.
|
||
|
|
||
|
start : Float
|
||
|
|
||
|
- For concentrated/point loads, start is the x coordinate
|
||
|
- For distributed loads, start is the starting position of distributed load
|
||
|
|
||
|
mag : Sympifyable
|
||
|
Magnitude of the applied load. Must be positive
|
||
|
|
||
|
end : Float
|
||
|
Required for distributed loads
|
||
|
|
||
|
- For concentrated/point load , end is None(may not be given)
|
||
|
- For distributed loads, end is the end position of distributed load
|
||
|
|
||
|
angle: Sympifyable
|
||
|
The angle in degrees, the load vector makes with the horizontal
|
||
|
in the counter-clockwise direction.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
For applying distributed load
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.apply_load(0,'C',start=3,end=5,mag=-10)
|
||
|
|
||
|
For applying point/concentrated_loads
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.apply_load(-1,'C',start=2,mag=15,angle=45)
|
||
|
|
||
|
"""
|
||
|
y = Symbol('y')
|
||
|
x = Symbol('x')
|
||
|
x0 = Symbol('x0')
|
||
|
# y0 = Symbol('y0')
|
||
|
order= sympify(order)
|
||
|
mag = sympify(mag)
|
||
|
angle = sympify(angle)
|
||
|
|
||
|
if label in self._loads_applied:
|
||
|
raise ValueError("load with the given label already exists")
|
||
|
|
||
|
if label in ['A','B']:
|
||
|
raise ValueError("cannot use the given label, reserved for supports")
|
||
|
|
||
|
if order == 0:
|
||
|
if end is None or end<start:
|
||
|
raise KeyError("provide end greater than start")
|
||
|
|
||
|
self._distributed_loads[label] = {'start':start, 'end':end, 'f_y': mag}
|
||
|
self._points_disc_y.add(start)
|
||
|
|
||
|
if start in self._moment_y:
|
||
|
self._moment_y[start] -= mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
|
||
|
self._load_y[start] += mag*(Min(end,x)-start)
|
||
|
else:
|
||
|
self._moment_y[start] = -mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
|
||
|
self._load_y[start] = mag*(Min(end,x)-start)
|
||
|
|
||
|
self._loads_applied[label] = 'distributed'
|
||
|
|
||
|
if order == -1:
|
||
|
|
||
|
if angle is None:
|
||
|
raise TypeError("please provide direction of force")
|
||
|
height = self._shape_eqn.subs({'x':start})
|
||
|
|
||
|
self._conc_loads[label] = {'x':start, 'y':height, 'f_x':mag*cos(rad(angle)), 'f_y': mag*sin(rad(angle)), 'mag':mag, 'angle':angle}
|
||
|
self._points_disc_x.add(start)
|
||
|
self._points_disc_y.add(start)
|
||
|
|
||
|
if start in self._moment_x:
|
||
|
self._moment_x[start] += self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
|
||
|
self._load_x[start] += self._conc_loads[label]['f_x']
|
||
|
else:
|
||
|
self._moment_x[start] = self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
|
||
|
self._load_x[start] = self._conc_loads[label]['f_x']
|
||
|
|
||
|
if start in self._moment_y:
|
||
|
self._moment_y[start] -= self._conc_loads[label]['f_y']*(x0-start)
|
||
|
self._load_y[start] += self._conc_loads[label]['f_y']
|
||
|
else:
|
||
|
self._moment_y[start] = -self._conc_loads[label]['f_y']*(x0-start)
|
||
|
self._load_y[start] = self._conc_loads[label]['f_y']
|
||
|
|
||
|
self._loads_applied[label] = 'concentrated'
|
||
|
|
||
|
|
||
|
def remove_load(self,label):
|
||
|
"""
|
||
|
This methods removes the load applied to the arch
|
||
|
|
||
|
Parameters
|
||
|
==========
|
||
|
|
||
|
label : String or Symbol
|
||
|
The label of the applied load
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.apply_load(0,'C',start=3,end=5,mag=-10)
|
||
|
>>> a.remove_load('C')
|
||
|
removed load C: {'start': 3, 'end': 5, 'f_y': -10}
|
||
|
"""
|
||
|
y = Symbol('y')
|
||
|
x = Symbol('x')
|
||
|
x0 = Symbol('x0')
|
||
|
|
||
|
if label in self._distributed_loads :
|
||
|
|
||
|
self._loads_applied.pop(label)
|
||
|
start = self._distributed_loads[label]['start']
|
||
|
end = self._distributed_loads[label]['end']
|
||
|
mag = self._distributed_loads[label]['f_y']
|
||
|
self._points_disc_y.remove(start)
|
||
|
self._load_y[start] -= mag*(Min(x,end)-start)
|
||
|
self._moment_y[start] += mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
|
||
|
val = self._distributed_loads.pop(label)
|
||
|
print(f"removed load {label}: {val}")
|
||
|
|
||
|
elif label in self._conc_loads :
|
||
|
|
||
|
self._loads_applied.pop(label)
|
||
|
start = self._conc_loads[label]['x']
|
||
|
self._points_disc_x.remove(start)
|
||
|
self._points_disc_y.remove(start)
|
||
|
self._moment_y[start] += self._conc_loads[label]['f_y']*(x0-start)
|
||
|
self._moment_x[start] -= self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
|
||
|
self._load_x[start] -= self._conc_loads[label]['f_x']
|
||
|
self._load_y[start] -= self._conc_loads[label]['f_y']
|
||
|
val = self._conc_loads.pop(label)
|
||
|
print(f"removed load {label}: {val}")
|
||
|
|
||
|
else :
|
||
|
raise ValueError("label not found")
|
||
|
|
||
|
def change_support_position(self, left_support=None, right_support=None):
|
||
|
"""
|
||
|
Change position of supports.
|
||
|
If not provided , defaults to the old value.
|
||
|
Parameters
|
||
|
==========
|
||
|
|
||
|
left_support: tuple (x, y)
|
||
|
x: float
|
||
|
x-coordinate value of the left_support
|
||
|
|
||
|
y: float
|
||
|
y-coordinate value of the left_support
|
||
|
|
||
|
right_support: tuple (x, y)
|
||
|
x: float
|
||
|
x-coordinate value of the right_support
|
||
|
|
||
|
y: float
|
||
|
y-coordinate value of the right_support
|
||
|
"""
|
||
|
if left_support is not None:
|
||
|
self._left_support = (left_support[0],left_support[1])
|
||
|
|
||
|
if right_support is not None:
|
||
|
self._right_support = (right_support[0],right_support[1])
|
||
|
|
||
|
self._shape_eqn = None
|
||
|
self._shape_eqn = self.get_shape_eqn
|
||
|
|
||
|
def change_crown_position(self,crown_x=None,crown_y=None):
|
||
|
"""
|
||
|
Change the position of the crown/hinge of the arch
|
||
|
|
||
|
Parameters
|
||
|
==========
|
||
|
|
||
|
crown_x: Float
|
||
|
The x coordinate of the position of the hinge
|
||
|
- if not provided, defaults to old value
|
||
|
|
||
|
crown_y: Float
|
||
|
The y coordinate of the position of the hinge
|
||
|
- if not provided defaults to None
|
||
|
"""
|
||
|
self._crown_x = crown_x
|
||
|
self._crown_y = crown_y
|
||
|
self._shape_eqn = None
|
||
|
self._shape_eqn = self.get_shape_eqn
|
||
|
|
||
|
def change_support_type(self,left_support=None,right_support=None):
|
||
|
"""
|
||
|
Add the type for support at each end.
|
||
|
Can use roller or hinge support at each end.
|
||
|
|
||
|
Parameters
|
||
|
==========
|
||
|
|
||
|
left_support, right_support : string
|
||
|
Type of support at respective end
|
||
|
|
||
|
- For roller support , left_support/right_support = "roller"
|
||
|
- For hinged support, left_support/right_support = "hinge"
|
||
|
- defaults to hinge if value not provided
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
For applying roller support at right end
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.change_support_type(right_support="roller")
|
||
|
|
||
|
"""
|
||
|
support_types = ['roller','hinge']
|
||
|
if left_support:
|
||
|
if left_support not in support_types:
|
||
|
raise ValueError("supports must only be roller or hinge")
|
||
|
|
||
|
self._supports['left'] = left_support
|
||
|
|
||
|
if right_support:
|
||
|
if right_support not in support_types:
|
||
|
raise ValueError("supports must only be roller or hinge")
|
||
|
|
||
|
self._supports['right'] = right_support
|
||
|
|
||
|
def add_member(self,y):
|
||
|
"""
|
||
|
This method adds a member/rod at a particular height y.
|
||
|
A rod is used for stability of the structure in case of a roller support.
|
||
|
"""
|
||
|
if y>self._crown_y or y<min(self._left_support[1], self._right_support[1]):
|
||
|
raise ValueError(f"position of support must be between y={min(self._left_support[1], self._right_support[1])} and y={self._crown_y}")
|
||
|
x = Symbol('x')
|
||
|
a = diff(self._shape_eqn,x).subs(x,self._crown_x+1)/2
|
||
|
x_diff = sqrt((y - self._crown_y)/a)
|
||
|
x1 = self._crown_x + x_diff
|
||
|
x2 = self._crown_x - x_diff
|
||
|
self._member = (x1,x2,y)
|
||
|
|
||
|
def shear_force_at(self, pos = None, **kwargs):
|
||
|
"""
|
||
|
return the shear at some x-coordinates
|
||
|
if no x value provided, returns the formula
|
||
|
"""
|
||
|
if pos is None:
|
||
|
return self._shear_force
|
||
|
else:
|
||
|
x = Symbol('x')
|
||
|
if 'dir' in kwargs:
|
||
|
dir = kwargs['dir']
|
||
|
return limit(self._shear_force,x,pos,dir=dir)
|
||
|
return self._shear_force.subs(x,pos)
|
||
|
|
||
|
def bending_moment_at(self, pos = None, **kwargs):
|
||
|
"""
|
||
|
return the bending moment at some x-coordinates
|
||
|
if no x value provided, returns the formula
|
||
|
"""
|
||
|
if pos is None:
|
||
|
return self._bending_moment
|
||
|
else:
|
||
|
x0 = Symbol('x0')
|
||
|
if 'dir' in kwargs:
|
||
|
dir = kwargs['dir']
|
||
|
return limit(self._bending_moment,x0,pos,dir=dir)
|
||
|
return self._bending_moment.subs(x0,pos)
|
||
|
|
||
|
|
||
|
def axial_force_at(self,pos = None, **kwargs):
|
||
|
"""
|
||
|
return the axial/normal force generated at some x-coordinate
|
||
|
if no x value provided, returns the formula
|
||
|
"""
|
||
|
if pos is None:
|
||
|
return self._axial_force
|
||
|
else:
|
||
|
x = Symbol('x')
|
||
|
if 'dir' in kwargs:
|
||
|
dir = kwargs['dir']
|
||
|
return limit(self._axial_force,x,pos,dir=dir)
|
||
|
return self._axial_force.subs(x,pos)
|
||
|
|
||
|
def solve(self):
|
||
|
"""
|
||
|
This method solves for the reaction forces generated at the supports,\n
|
||
|
and bending moment and generated in the arch and tension produced in the member if used.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
|
||
|
>>> a.apply_load(0,'C',start=3,end=5,mag=-10)
|
||
|
>>> a.solve()
|
||
|
>>> a.reaction_force
|
||
|
{R_A_x: 8, R_A_y: 12, R_B_x: -8, R_B_y: 8}
|
||
|
|
||
|
>>> from sympy import Symbol
|
||
|
>>> t = Symbol('t')
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(16,0),crown_x=8,crown_y=5)
|
||
|
>>> a.apply_load(0,'C',start=3,end=5,mag=t)
|
||
|
>>> a.solve()
|
||
|
>>> a.reaction_force
|
||
|
{R_A_x: -4*t/5, R_A_y: -3*t/2, R_B_x: 4*t/5, R_B_y: -t/2}
|
||
|
|
||
|
>>> a.bending_moment_at(4)
|
||
|
-5*t/2
|
||
|
"""
|
||
|
y = Symbol('y')
|
||
|
x = Symbol('x')
|
||
|
x0 = Symbol('x0')
|
||
|
|
||
|
discontinuity_points_x = sorted(self._points_disc_x)
|
||
|
discontinuity_points_y = sorted(self._points_disc_y)
|
||
|
|
||
|
self._moment_x_func = Piecewise((0,True))
|
||
|
self._moment_y_func = Piecewise((0,True))
|
||
|
|
||
|
self._load_x_func = Piecewise((0,True))
|
||
|
self._load_y_func = Piecewise((0,True))
|
||
|
|
||
|
accumulated_x_moment = 0
|
||
|
accumulated_y_moment = 0
|
||
|
|
||
|
accumulated_x_load = 0
|
||
|
accumulated_y_load = 0
|
||
|
|
||
|
for point in discontinuity_points_x:
|
||
|
cond = (x >= point)
|
||
|
accumulated_x_load += self._load_x[point]
|
||
|
accumulated_x_moment += self._moment_x[point]
|
||
|
self._load_x_func = Piecewise((accumulated_x_load,cond),(self._load_x_func,True))
|
||
|
self._moment_x_func = Piecewise((accumulated_x_moment,cond),(self._moment_x_func,True))
|
||
|
|
||
|
for point in discontinuity_points_y:
|
||
|
cond = (x >= point)
|
||
|
accumulated_y_moment += self._moment_y[point]
|
||
|
accumulated_y_load += self._load_y[point]
|
||
|
self._load_y_func = Piecewise((accumulated_y_load,cond),(self._load_y_func,True))
|
||
|
self._moment_y_func = Piecewise((accumulated_y_moment,cond),(self._moment_y_func,True))
|
||
|
|
||
|
moment_A = self._moment_y_func.subs(x,self._right_support[0]).subs(x0,self._left_support[0]) +\
|
||
|
self._moment_x_func.subs(x,self._right_support[0]).subs(y,self._left_support[1])
|
||
|
|
||
|
moment_hinge_left = self._moment_y_func.subs(x,self._crown_x).subs(x0,self._crown_x) +\
|
||
|
self._moment_x_func.subs(x,self._crown_x).subs(y,self._crown_y)
|
||
|
|
||
|
moment_hinge_right = self._moment_y_func.subs(x,self._right_support[0]).subs(x0,self._crown_x)- \
|
||
|
self._moment_y_func.subs(x,self._crown_x).subs(x0,self._crown_x) +\
|
||
|
self._moment_x_func.subs(x,self._right_support[0]).subs(y,self._crown_y) -\
|
||
|
self._moment_x_func.subs(x,self._crown_x).subs(y,self._crown_y)
|
||
|
|
||
|
net_x = self._load_x_func.subs(x,self._right_support[0])
|
||
|
net_y = self._load_y_func.subs(x,self._right_support[0])
|
||
|
|
||
|
if (self._supports['left']=='roller' or self._supports['right']=='roller') and not self._member:
|
||
|
print("member must be added if any of the supports is roller")
|
||
|
return
|
||
|
|
||
|
R_A_x, R_A_y, R_B_x, R_B_y, T = symbols('R_A_x R_A_y R_B_x R_B_y T')
|
||
|
|
||
|
if self._supports['left'] == 'roller' and self._supports['right'] == 'roller':
|
||
|
|
||
|
if self._member[2]>=max(self._left_support[1],self._right_support[1]):
|
||
|
|
||
|
if net_x!=0:
|
||
|
raise ValueError("net force in x direction not possible under the specified conditions")
|
||
|
|
||
|
else:
|
||
|
eq1 = Eq(R_A_x ,0)
|
||
|
eq2 = Eq(R_B_x, 0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
|
||
|
eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)
|
||
|
|
||
|
eq5 = Eq(moment_hinge_right + R_B_y*(self._right_support[0]-self._crown_x) +\
|
||
|
T*(self._member[2]-self._crown_y),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._left_support[1]:
|
||
|
eq1 = Eq(R_A_x ,0)
|
||
|
eq2 = Eq(R_B_x, 0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
T*(self._member[2]-self._left_support[1])+moment_A,0)
|
||
|
eq5 = Eq(T+net_x,0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._right_support[1]:
|
||
|
eq1 = Eq(R_A_x ,0)
|
||
|
eq2 = Eq(R_B_x, 0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])+\
|
||
|
T*(self._member[2]-self._left_support[1])+moment_A,0)
|
||
|
eq5 = Eq(T-net_x,0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._supports['left'] == 'roller':
|
||
|
if self._member[2]>=max(self._left_support[1], self._right_support[1]):
|
||
|
eq1 = Eq(R_A_x ,0)
|
||
|
eq2 = Eq(R_B_x+net_x,0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)
|
||
|
eq5 = Eq(moment_hinge_left + R_A_y*(self._left_support[0]-self._crown_x) -\
|
||
|
T*(self._member[2]-self._crown_y),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._left_support[1]:
|
||
|
eq1 = Eq(R_A_x ,0)
|
||
|
eq2 = Eq(R_B_x+ T +net_x,0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
R_B_x*(self._right_support[1]-self._left_support[1])-\
|
||
|
T*(self._member[2]-self._left_support[0])+moment_A,0)
|
||
|
eq5 = Eq(moment_hinge_left + R_A_y*(self._left_support[0]-self._crown_x)-\
|
||
|
T*(self._member[2]-self._crown_y),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._right_support[0]:
|
||
|
eq1 = Eq(R_A_x,0)
|
||
|
eq2 = Eq(R_B_x- T +net_x,0)
|
||
|
eq3 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq4 = Eq(moment_hinge_left+R_A_y*(self._left_support[0]-self._crown_x),0)
|
||
|
eq5 = Eq(moment_A+R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
R_B_x*(self._right_support[1]-self._left_support[1])+\
|
||
|
T*(self._member[2]-self._left_support[1]),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._supports['right'] == 'roller':
|
||
|
if self._member[2]>=max(self._left_support[1], self._right_support[1]):
|
||
|
eq1 = Eq(R_B_x,0)
|
||
|
eq2 = Eq(R_A_x+net_x,0)
|
||
|
eq3 = Eq(R_A_y+R_B_y+net_y,0)
|
||
|
eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x)+\
|
||
|
T*(self._member[2]-self._crown_y),0)
|
||
|
eq5 = Eq(moment_A+R_B_y*(self._right_support[0]-self._left_support[0]),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._left_support[1]:
|
||
|
eq1 = Eq(R_B_x,0)
|
||
|
eq2 = Eq(R_A_x+T+net_x,0)
|
||
|
eq3 = Eq(R_A_y+R_B_y+net_y,0)
|
||
|
eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x),0)
|
||
|
eq5 = Eq(moment_A-T*(self._member[2]-self._left_support[1])+\
|
||
|
R_B_y*(self._right_support[0]-self._left_support[0]),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
|
||
|
elif self._member[2]>=self._right_support[1]:
|
||
|
eq1 = Eq(R_B_x,0)
|
||
|
eq2 = Eq(R_A_x-T+net_x,0)
|
||
|
eq3 = Eq(R_A_y+R_B_y+net_y,0)
|
||
|
eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x)+\
|
||
|
T*(self._member[2]-self._crown_y),0)
|
||
|
eq5 = Eq(moment_A+T*(self._member[2]-self._left_support[1])+\
|
||
|
R_B_y*(self._right_support[0]-self._left_support[0]))
|
||
|
solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
|
||
|
else:
|
||
|
eq1 = Eq(R_A_x + R_B_x + net_x,0)
|
||
|
eq2 = Eq(R_A_y + R_B_y + net_y,0)
|
||
|
eq3 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
|
||
|
R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)
|
||
|
eq4 = Eq(moment_hinge_right + R_B_y*(self._right_support[0]-self._crown_x) -\
|
||
|
R_B_x*(self._right_support[1]-self._crown_y),0)
|
||
|
solution = solve((eq1,eq2,eq3,eq4),(R_A_x,R_A_y,R_B_x,R_B_y))
|
||
|
|
||
|
for symb in self._reaction_force:
|
||
|
self._reaction_force[symb] = solution[symb]
|
||
|
|
||
|
self._bending_moment = - (self._moment_x_func.subs(x,x0) + self._moment_y_func.subs(x,x0) -\
|
||
|
solution[R_A_y]*(x0-self._left_support[0]) +\
|
||
|
solution[R_A_x]*(self._shape_eqn.subs({x:x0})-self._left_support[1]))
|
||
|
|
||
|
angle = atan(diff(self._shape_eqn,x))
|
||
|
|
||
|
fx = (self._load_x_func+solution[R_A_x])
|
||
|
fy = (self._load_y_func+solution[R_A_y])
|
||
|
|
||
|
axial_force = fx*cos(angle) + fy*sin(angle)
|
||
|
shear_force = -fx*sin(angle) + fy*cos(angle)
|
||
|
|
||
|
self._axial_force = axial_force
|
||
|
self._shear_force = shear_force
|
||
|
|
||
|
@doctest_depends_on(modules=('numpy',))
|
||
|
def draw(self):
|
||
|
"""
|
||
|
This method returns a plot object containing the diagram of the specified arch along with the supports
|
||
|
and forces applied to the structure.
|
||
|
|
||
|
Examples
|
||
|
========
|
||
|
|
||
|
>>> from sympy import Symbol
|
||
|
>>> t = Symbol('t')
|
||
|
>>> from sympy.physics.continuum_mechanics.arch import Arch
|
||
|
>>> a = Arch((0,0),(40,0),crown_x=20,crown_y=12)
|
||
|
>>> a.apply_load(-1,'C',8,150,angle=270)
|
||
|
>>> a.apply_load(0,'D',start=20,end=40,mag=-4)
|
||
|
>>> a.apply_load(-1,'E',10,t,angle=300)
|
||
|
>>> p = a.draw()
|
||
|
>>> p # doctest: +ELLIPSIS
|
||
|
Plot object containing:
|
||
|
[0]: cartesian line: 11.325 - 3*(x - 20)**2/100 for x over (0.0, 40.0)
|
||
|
[1]: cartesian line: 12 - 3*(x - 20)**2/100 for x over (0.0, 40.0)
|
||
|
...
|
||
|
>>> p.show()
|
||
|
|
||
|
"""
|
||
|
x = Symbol('x')
|
||
|
markers = []
|
||
|
annotations = self._draw_loads()
|
||
|
rectangles = []
|
||
|
supports = self._draw_supports()
|
||
|
markers+=supports
|
||
|
|
||
|
xmax = self._right_support[0]
|
||
|
xmin = self._left_support[0]
|
||
|
ymin = min(self._left_support[1],self._right_support[1])
|
||
|
ymax = self._crown_y
|
||
|
|
||
|
lim = max(xmax*1.1-xmin*0.8+1, ymax*1.1-ymin*0.8+1)
|
||
|
|
||
|
rectangles = self._draw_rectangles()
|
||
|
|
||
|
filler = self._draw_filler()
|
||
|
rectangles+=filler
|
||
|
|
||
|
if self._member is not None:
|
||
|
if(self._member[2]>=self._right_support[1]):
|
||
|
markers.append(
|
||
|
{
|
||
|
'args':[[self._member[1]+0.005*lim],[self._member[2]]],
|
||
|
'marker':'o',
|
||
|
'markersize': 4,
|
||
|
'color': 'white',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
if(self._member[2]>=self._left_support[1]):
|
||
|
markers.append(
|
||
|
{
|
||
|
'args':[[self._member[0]-0.005*lim],[self._member[2]]],
|
||
|
'marker':'o',
|
||
|
'markersize': 4,
|
||
|
'color': 'white',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
|
||
|
|
||
|
markers.append({
|
||
|
'args':[[self._crown_x],[self._crown_y-0.005*lim]],
|
||
|
'marker':'o',
|
||
|
'markersize': 5,
|
||
|
'color':'white',
|
||
|
'markerfacecolor':'none',
|
||
|
})
|
||
|
|
||
|
if lim==xmax*1.1-xmin*0.8+1:
|
||
|
|
||
|
sing_plot = plot(self._shape_eqn-0.015*lim,
|
||
|
self._shape_eqn,
|
||
|
(x, self._left_support[0], self._right_support[0]),
|
||
|
markers=markers,
|
||
|
show=False,
|
||
|
annotations=annotations,
|
||
|
rectangles = rectangles,
|
||
|
xlim=(xmin-0.05*lim, xmax*1.1),
|
||
|
ylim=(xmin-0.05*lim, xmax*1.1),
|
||
|
axis=False,
|
||
|
line_color='brown')
|
||
|
|
||
|
else:
|
||
|
sing_plot = plot(self._shape_eqn-0.015*lim,
|
||
|
self._shape_eqn,
|
||
|
(x, self._left_support[0], self._right_support[0]),
|
||
|
markers=markers,
|
||
|
show=False,
|
||
|
annotations=annotations,
|
||
|
rectangles = rectangles,
|
||
|
xlim=(ymin-0.05*lim, ymax*1.1),
|
||
|
ylim=(ymin-0.05*lim, ymax*1.1),
|
||
|
axis=False,
|
||
|
line_color='brown')
|
||
|
|
||
|
return sing_plot
|
||
|
|
||
|
|
||
|
def _draw_supports(self):
|
||
|
support_markers = []
|
||
|
|
||
|
xmax = self._right_support[0]
|
||
|
xmin = self._left_support[0]
|
||
|
ymin = min(self._left_support[1],self._right_support[1])
|
||
|
ymax = self._crown_y
|
||
|
|
||
|
if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
|
||
|
max_diff = 1.1*xmax-0.8*xmin
|
||
|
else:
|
||
|
max_diff = 1.1*ymax-0.8*ymin
|
||
|
|
||
|
if self._supports['left']=='roller':
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._left_support[0]],
|
||
|
[self._left_support[1]-0.02*max_diff]
|
||
|
],
|
||
|
'marker':'o',
|
||
|
'markersize':11,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
else:
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._left_support[0]],
|
||
|
[self._left_support[1]-0.007*max_diff]
|
||
|
],
|
||
|
'marker':6,
|
||
|
'markersize':15,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
if self._supports['right']=='roller':
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._right_support[0]],
|
||
|
[self._right_support[1]-0.02*max_diff]
|
||
|
],
|
||
|
'marker':'o',
|
||
|
'markersize':11,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
else:
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._right_support[0]],
|
||
|
[self._right_support[1]-0.007*max_diff]
|
||
|
],
|
||
|
'marker':6,
|
||
|
'markersize':15,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._right_support[0]],
|
||
|
[self._right_support[1]-0.036*max_diff]
|
||
|
],
|
||
|
'marker':'_',
|
||
|
'markersize':15,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
support_markers.append(
|
||
|
{
|
||
|
'args':[
|
||
|
[self._left_support[0]],
|
||
|
[self._left_support[1]-0.036*max_diff]
|
||
|
],
|
||
|
'marker':'_',
|
||
|
'markersize':15,
|
||
|
'color':'black',
|
||
|
'markerfacecolor':'none'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
return support_markers
|
||
|
|
||
|
def _draw_rectangles(self):
|
||
|
member = []
|
||
|
|
||
|
xmax = self._right_support[0]
|
||
|
xmin = self._left_support[0]
|
||
|
ymin = min(self._left_support[1],self._right_support[1])
|
||
|
ymax = self._crown_y
|
||
|
|
||
|
if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
|
||
|
max_diff = 1.1*xmax-0.8*xmin
|
||
|
else:
|
||
|
max_diff = 1.1*ymax-0.8*ymin
|
||
|
|
||
|
if self._member is not None:
|
||
|
if self._member[2]>= max(self._left_support[1],self._right_support[1]):
|
||
|
member.append(
|
||
|
{
|
||
|
'xy':(self._member[0],self._member[2]-0.005*max_diff),
|
||
|
'width':self._member[1]-self._member[0],
|
||
|
'height': 0.01*max_diff,
|
||
|
'angle': 0,
|
||
|
'color':'brown',
|
||
|
}
|
||
|
)
|
||
|
|
||
|
elif self._member[2]>=self._left_support[1]:
|
||
|
member.append(
|
||
|
{
|
||
|
'xy':(self._member[0],self._member[2]-0.005*max_diff),
|
||
|
'width':self._right_support[0]-self._member[0],
|
||
|
'height': 0.01*max_diff,
|
||
|
'angle': 0,
|
||
|
'color':'brown',
|
||
|
}
|
||
|
)
|
||
|
|
||
|
else:
|
||
|
member.append(
|
||
|
{
|
||
|
'xy':(self._member[1],self._member[2]-0.005*max_diff),
|
||
|
'width':abs(self._left_support[0]-self._member[1]),
|
||
|
'height': 0.01*max_diff,
|
||
|
'angle': 180,
|
||
|
'color':'brown',
|
||
|
}
|
||
|
)
|
||
|
|
||
|
if self._distributed_loads:
|
||
|
for loads in self._distributed_loads:
|
||
|
|
||
|
start = self._distributed_loads[loads]['start']
|
||
|
end = self._distributed_loads[loads]['end']
|
||
|
|
||
|
member.append(
|
||
|
{
|
||
|
'xy':(start,self._crown_y+max_diff*0.15),
|
||
|
'width': (end-start),
|
||
|
'height': max_diff*0.01,
|
||
|
'color': 'orange'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
|
||
|
return member
|
||
|
|
||
|
def _draw_loads(self):
|
||
|
load_annotations = []
|
||
|
|
||
|
xmax = self._right_support[0]
|
||
|
xmin = self._left_support[0]
|
||
|
ymin = min(self._left_support[1],self._right_support[1])
|
||
|
ymax = self._crown_y
|
||
|
|
||
|
if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
|
||
|
max_diff = 1.1*xmax-0.8*xmin
|
||
|
else:
|
||
|
max_diff = 1.1*ymax-0.8*ymin
|
||
|
|
||
|
for load in self._conc_loads:
|
||
|
x = self._conc_loads[load]['x']
|
||
|
y = self._conc_loads[load]['y']
|
||
|
angle = self._conc_loads[load]['angle']
|
||
|
mag = self._conc_loads[load]['mag']
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':'',
|
||
|
'xy':(
|
||
|
x+cos(rad(angle))*max_diff*0.08,
|
||
|
y+sin(rad(angle))*max_diff*0.08
|
||
|
),
|
||
|
'xytext':(x,y),
|
||
|
'fontsize':10,
|
||
|
'fontweight': 'bold',
|
||
|
'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'blue','edgecolor':'blue'}
|
||
|
}
|
||
|
)
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':f'{load}: {mag} N',
|
||
|
'fontsize':10,
|
||
|
'fontweight': 'bold',
|
||
|
'xy': (x+cos(rad(angle))*max_diff*0.12,y+sin(rad(angle))*max_diff*0.12)
|
||
|
}
|
||
|
)
|
||
|
|
||
|
for load in self._distributed_loads:
|
||
|
start = self._distributed_loads[load]['start']
|
||
|
end = self._distributed_loads[load]['end']
|
||
|
mag = self._distributed_loads[load]['f_y']
|
||
|
x_points = numpy.arange(start,end,(end-start)/(max_diff*0.25))
|
||
|
x_points = numpy.append(x_points,end)
|
||
|
for point in x_points:
|
||
|
if(mag<0):
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':'',
|
||
|
'xy':(point,self._crown_y+max_diff*0.05),
|
||
|
'xytext': (point,self._crown_y+max_diff*0.15),
|
||
|
'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'orange','edgecolor':'orange'}
|
||
|
}
|
||
|
)
|
||
|
else:
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':'',
|
||
|
'xy':(point,self._crown_y+max_diff*0.2),
|
||
|
'xytext': (point,self._crown_y+max_diff*0.15),
|
||
|
'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'orange','edgecolor':'orange'}
|
||
|
}
|
||
|
)
|
||
|
if(mag<0):
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':f'{load}: {abs(mag)} N/m',
|
||
|
'fontsize':10,
|
||
|
'fontweight': 'bold',
|
||
|
'xy':((start+end)/2,self._crown_y+max_diff*0.175)
|
||
|
}
|
||
|
)
|
||
|
else:
|
||
|
load_annotations.append(
|
||
|
{
|
||
|
'text':f'{load}: {abs(mag)} N/m',
|
||
|
'fontsize':10,
|
||
|
'fontweight': 'bold',
|
||
|
'xy':((start+end)/2,self._crown_y+max_diff*0.125)
|
||
|
}
|
||
|
)
|
||
|
return load_annotations
|
||
|
|
||
|
def _draw_filler(self):
|
||
|
x = Symbol('x')
|
||
|
filler = []
|
||
|
xmax = self._right_support[0]
|
||
|
xmin = self._left_support[0]
|
||
|
ymin = min(self._left_support[1],self._right_support[1])
|
||
|
ymax = self._crown_y
|
||
|
|
||
|
if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
|
||
|
max_diff = 1.1*xmax-0.8*xmin
|
||
|
else:
|
||
|
max_diff = 1.1*ymax-0.8*ymin
|
||
|
|
||
|
x_points = numpy.arange(self._left_support[0],self._right_support[0],(self._right_support[0]-self._left_support[0])/(max_diff*max_diff))
|
||
|
|
||
|
for point in x_points:
|
||
|
filler.append(
|
||
|
{
|
||
|
'xy':(point,self._shape_eqn.subs(x,point)-max_diff*0.015),
|
||
|
'width': (self._right_support[0]-self._left_support[0])/(max_diff*max_diff),
|
||
|
'height': max_diff*0.015,
|
||
|
'color': 'brown'
|
||
|
}
|
||
|
)
|
||
|
|
||
|
return filler
|