New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update beam.py to allow solving reaction moments #14681

Closed
wants to merge 6 commits into
base: master
from

Conversation

Projects
None yet
4 participants
@cym1
Contributor

cym1 commented May 1, 2018

The previous version doesn't allow applying reaction moments.

Example, a beam has length l and fixed supports at both ends. A 1000N load is applied at the center.

E = Symbol('E')
I = Symbol('I')
l = Symbol('l', positive=True)
R1, R2 = symbols('R1, R2')
M1, M2 = symbols('M1, M2')

b = Beam(l, E, I)
b.bc_deflection = [(0, 0),(l, 0)]
b.bc_slope = [(0, 0),(l, 0)]

b.apply_load(R1, 0, -1)
b.apply_load(M1, 0, -2)
b.apply_load(R2, l, -1)
b.apply_load(M2, l, -2)
b.apply_load(-1000, 0.5*l, -1)

Method solve_for_reaction_loads will not solve M1 and M2. A new method solve() is proposed here to find M1, M2.

A new method extrema() also proposed to find the maximum deflection points

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Update beam.py
The previous version doesn't allow applying reaction moments.

Example, a beam has length l and fixed supports at both ends. A 1000N load is applied at the center.

E = Symbol('E')
I = Symbol('I')
l = Symbol('l', positive=True)
R1, R2 = symbols('R1, R2')
M1, M2 = symbols('M1, M2')

b = Beam(l, E, I)
b.bc_deflection = [(0, 0),(l, 0)]
b.bc_slope = [(0, 0),(l, 0)]

b.apply_load(R1, 0, -1)
b.apply_load(M1, 0, -2)
b.apply_load(R2, l, -1)
b.apply_load(M2, l, -2)
b.apply_load(-1000, 0.5*l, -1)

Method solve_for_reaction_loads will not solve M1 and M2. A new method solve() is proposed here to find M1, M2.

A new method extrema() also proposed to find the maximum deflection points
@moorepants

This comment has been minimized.

Member

moorepants commented May 1, 2018

These look like great improvements. It will need some tuning to get merged. First is that you need to follow the deprecation policy for API changes (i.e. removing the method). Secondly, please improve the docstrings for the new methods, ideally with parameters, returns, examples etc.

@moorepants

This comment has been minimized.

Member

moorepants commented May 1, 2018

@moorepants

This comment has been minimized.

Member

moorepants commented May 1, 2018

@jashan498 It would be good to have you review this.

self._reaction_loads = {}
self.C3 = Symbol('C3')
self.C4 = Symbol('C4')

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

is there any particular reason to declare C3 and C4 inside the constructor? I think declaring them inside the function itself should work.

This comment has been minimized.

@cym1

cym1 May 2, 2018

Contributor

If number of fixed supports >= 2, number of bc_slope applied will be also >=2, when solving C3 in slope(), the system of equations will be in the form of:
[C3, C3 + an expression with symbols R1&M1]

The solution is C3 =0 or an expression of R1&M1, linsolve will return an empty set.

In this case, reactions R1 and M1 must be solved before calling slope(). In the process of solving R1 and M1, C3 and C4 are also solved, so why don't just remember the value of C3 and C4 such that re-calculation is not needed in slope() and deflection()?

This comment has been minimized.

@jashan498

jashan498 May 2, 2018

Contributor

Yeah it does make sense, but I think there would be some test failures due to this. Let's see if there are some during travis build up.

-8*SingularityFunction(x, 0, -1) + 6*SingularityFunction(x, 10, -1)
+ 120*SingularityFunction(x, 30, -2) + 2*SingularityFunction(x, 30, -1)
"""
def solve(self):

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

I think we can go with the same name. In that way, we don't need to deprecate the API.

else:
for position, value in self._boundary_conditions['slope']:
eqs = sympify(slope_curve.subs(x, position) - value)

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

Shouldn't there be simplify instead of sympify?

This comment has been minimized.

@moorepants

moorepants May 1, 2018

Member

Best to avoid simplify calls in library code, as they can be time consuming. It is better to let the user call simplify if they want to.

bc_eqs.append(eqs)
constants = list(linsolve(bc_eqs, C4))
deflection_curve = deflection_curve.subs({C4: constants[0][0]})

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

After removing the above the chunk of code, how would we find the values of integration constants(C3 and C4)?

deflection_curve = deflection_curve.subs({C4: constants[0][0]})
slope_curve = integrate(self.bending_moment(), x) + self.C3
deflection_curve = integrate(slope_curve, x) + self.C4

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

C3 and C4 can be evaluated by using self._boundary_conditions['slope'] and self._boundary_conditions['deflection'] respectively.

def extrema(self):
"""
Return a list of extrema in form of tuple of x and deflection

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

Also, add an example for docstring.

return S(1)/(E*I)*deflection_curve
def extrema(self):

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

Can we name it something like max_delfection to make it clear whose extrema we are talking about?

else:
deflection_curve = integrate(slope_curve, x) + C4
for position, value in self._boundary_conditions['deflection']:

This comment has been minimized.

@jashan498

jashan498 May 1, 2018

Contributor

I think there is no need to keep this loop inside if-else block. If not self._boundary_conditions['deflection'] is true then the loop would never be used.
Same goes for above loop.

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 1, 2018

@cym1 this issue is not specific to reaction moments as you mentioned in the description.
As far I know this applies to all the beam systems having more than 2 reaction forces(as the previous code used only 2 equations to solve for reactions). I proposed something similar here and i think this is what your code is doing.

Also you need to add required test cases to show that this PR does its work here in this file.

@cym1 cym1 closed this May 2, 2018

@@ -432,27 +462,17 @@ def slope(self):
>>> 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.solve()

This comment has been minimized.

@jashan498

jashan498 May 2, 2018

Contributor

Now as we are keeping the same name, you can change it back to solve_for_reaction_loads

This comment has been minimized.

@moorepants

moorepants May 2, 2018

Member

Yes, I think we shouldn't change the name.

def max_deflection(self):
"""
Return the maximum deflection point in form of tuple of x and deflection

This comment has been minimized.

@jashan498

jashan498 May 2, 2018

Contributor

You should also include a working example in the method's docstring(see other methods in the class).

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 2, 2018

Also, add few test cases to show that your PR can solve the problem it proposes.

@cym1 cym1 reopened this May 2, 2018

cym1 added some commits May 2, 2018

assert p == q
p = b5.max_deflection()
q = (l/2, -F*l**3/(192*E*I))

This comment has been minimized.

@cym1

cym1 May 2, 2018

Contributor

The answer can be checked here

cym1 added some commits May 2, 2018

@Abdullahjavednesar

This comment has been minimized.

Contributor

Abdullahjavednesar commented May 5, 2018

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 5, 2018

except for the failing test, this PR looks fine to me

@Abdullahjavednesar

This comment has been minimized.

Contributor

Abdullahjavednesar commented May 6, 2018

test_beam.py seems to fail.

@cym1

This comment has been minimized.

Contributor

cym1 commented May 8, 2018

The codes of solvers.py changed in sympy-1.1.2.dev.
It causes error as there are interval solutions (x beyond the boundaries) in solving slope() = 0. My idea may be replacing the value of slope() at x beyond the boundaries to S.NaN right before passing it to solve.

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 8, 2018

I think it just meant solution to be coming as Interval instance? Does using solveset instead of solve helps?

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 8, 2018

For example

>>> solveset(SingularityFunction(x,2,1).rewrite(Piecewise), x,domain=S.Reals)
Interval(-oo, 2)

but using solve for the same gives NotImplementedError:

>>> solve(SingularityFunction(x,2,1).rewrite(Piecewise), x,domain=S.Reals)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jashan/sympy/sympy/solvers/solvers.py", line 1155, in solve
    solution = _solve(f[0], *symbols, **flags)
  File "/home/jashan/sympy/sympy/solvers/solvers.py", line 1458, in _solve
    'solve cannot represent interval solutions')
NotImplementedError: solve cannot represent interval solutions
@cym1

This comment has been minimized.

Contributor

cym1 commented May 9, 2018

solveset doesn't work for multivariate expressions, e.g. SingularityFunction(x,b,1)

It requires a new multivariate as_set() in sympy\logic\boolalg.py

@jashan498

This comment has been minimized.

Contributor

jashan498 commented May 29, 2018

Sorry for the late response, was busy with exams.
Btw see if this diff is of any help. Worked for me in #14753

+ moment_curve = Piecewise((float("nan"), self.variable<=0),
+                (self.bending_moment(), self.variable<self.length),
+                (float("nan"), True))

- extreme_points = solve(self.slope().rewrite(Piecewise), self.variable, domain=S.Reals)
+ points = solve(moment_curve.rewrite(Piecewise), self.variable, domain=S.Reals)
@@ -171,3 +171,26 @@ def test_Beam():
raises(ValueError, lambda: b4.apply_load(-3, 0, -1, end=3))
with raises(TypeError):
b4.variable = 1
M1, M2 = symbols('M1, M2')
F = Symbol('F')

This comment has been minimized.

@jashan498

jashan498 May 29, 2018

Contributor

I think we should create separate test functions for both the methods(something like test_statically_indeterminate() and test_max_deflection()).
All tests inside test_beam function look confusing.

@cym1

This comment has been minimized.

Contributor

cym1 commented May 29, 2018

This works when the boundaries/fixtures are applied at x= 0 or x= length.

But someone could apply fixed support at other points, e.g. if someone applies fixed support at x = length/10, then there will still be an interval solution between 0 and length/10.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment