Skip to content
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

Introduced a class Column which provides functionality for column buckling calculations #17122

Open
wants to merge 6 commits into
base: master
from

Conversation

Projects
None yet
4 participants
@ishanaj
Copy link
Contributor

commented Jun 28, 2019

As discussed in issue #17072 , I have introduced a non-mutable class Column capable of doing calculations related to column buckling.
This class takes in all the required input during its object creation (i.e in the __init__), and further calculates the moment, slope, deflection equations accordingly.

TODO's:

  • method to determine critical load
  • method to determine the unknown reactions, with the help of boundary conditions (In another PR)
  • tests
  • documentation

Release Notes

  • physics.continuum_mechanics
    • A new class Column has been introduced in the continuum mechanics module which is capable of performing calculations related to column buckling
@sympy-bot

This comment has been minimized.

Copy link

commented Jun 28, 2019

Hi, I am the SymPy bot (v147). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • physics.continuum_mechanics
    • A new class Column has been introduced in the continuum mechanics module which is capable of performing calculations related to column buckling (#17122 by @ishanaj)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.5.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

As discussed in issue #17072 , I have introduced a non-mutable class `Column` capable of doing calculations related to column buckling.
This class takes in all the required input during its object creation (i.e in the `__init__`), and further calculates the moment, slope, deflection equations accordingly.

TODO's:

- [x] method to determine critical load
- [x] ~~method to determine the unknown reactions, with the help of boundary conditions~~ (In another PR)
- [x] tests
- [x] documentation

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
- physics.continuum_mechanics
    - A new class `Column` has been introduced in the continuum mechanics module which is capable of performing calculations related to column buckling
<!-- END RELEASE NOTES -->

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 28, 2019

Example:

>>> l, E, I  = symbols('l, E, I')
>>> P = Symbol('P')
>>> c = Column(l, E, I, P, top="pinned", bottom="fixed")
>>> c.solve_slope_deflection()
>>> c.moment()
-F*(l - x) + P*y(x)
>>> c.slope()
-F/P - sqrt(-P/(E*I))*(-F*l/(2*P) - F/(2*P*sqrt(-P/(E*I))))*exp(-x*sqrt(-P/(E*I))) + sqrt(-P/(E*I))*(-F*l/(2*P) + F/(2*P*sqrt(-P/(E*I))))*exp(x*sqrt(-P/(E*I)))
>>> c.deflection()
F*l/P - F*x/P + (-F*l/(2*P) - F/(2*P*sqrt(-P/(E*I))))*exp(-x*sqrt(-P/(E*I))) + (-F*l/(2*P) + F/(2*P*sqrt(-P/(E*I))))*exp(x*sqrt(-P/(E*I)))

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jun 28, 2019

After the differential equation is solved the solution to y(x) comes out to be:

C1*exp(-x*sqrt(-P/(E*I))) + C2*exp(x*sqrt(-P/(E*I))) + F*l/P - F*x/P

We can convert this equation in the form of B1*cos(z) + B2*sin(z) + F*l/P - F*x/P, where z is x*sqrt(P/(E*I)) and this is the general form that we use, but I am not able to convert it into it into this form.
The importance of converting it into this form is that it is easier to solve it for critical load using the boundary conditions.
Otherwise, with the exp form, after the C1, C2 are substituted, it becomes a lenghty eq, and after substituting the boundary condition (at x = l, y = 0) it becomes an eq that even solve or solveset is not able to solve it for load.
I have tried rewrite(), it converts it into hyperbolic cosh and sinh which again doesn't help.
Is there a way to convert it into a simpler form?
Or any other way to determine critical load?

ping @jashan498 @moorepants @oscarbenjamin

@jashan498

This comment has been minimized.

Copy link
Contributor

commented Jun 28, 2019

I'm not sure but if P is positive then the terms can be expressed as C1*e^(I*x) + C2*e^(-I*x) which should get converted into cos and sin.

>>> expr = C1*exp(-x*sqrt(-P/(E*I2))) + C2*exp(x*sqrt(-P/(E*I2))) + F*l/P - F*x/P
>>> expr2 = expr.subs(sqrt(-P/(E*I2)), I*sqrt(P/(E*I2)))
>>> pprint(expr2)
             ______               ______            
            ╱  P                 ╱  P               
    -ⅈ⋅x⋅  ╱  ────        ⅈ⋅x⋅  ╱  ────             
         ╲╱   E⋅I₂            ╲╱   E⋅I₂    F⋅l   F⋅x
C₁⋅ℯ                + C₂⋅ℯ               + ─── - ───
                                            P     P 

>>> pprint(expr2.rewrite(cos).expand())
          ⎛      ______⎞         ⎛      ______⎞           ⎛      ______⎞      
          ⎜     ╱  P   ⎟         ⎜     ╱  P   ⎟           ⎜     ╱  P   ⎟      
- ⅈ⋅C₁⋅sin⎜x⋅  ╱  ──── ⎟ + C₁⋅cos⎜x⋅  ╱  ──── ⎟ + ⅈ⋅C₂⋅sin⎜x⋅  ╱  ──── ⎟ + C₂⋅
          ⎝  ╲╱   E⋅I₂ ⎠         ⎝  ╲╱   E⋅I₂ ⎠           ⎝  ╲╱   E⋅I₂ ⎠      

   ⎛      ______⎞            
   ⎜     ╱  P   ⎟   F⋅l   F⋅x
cos⎜x⋅  ╱  ──── ⎟ + ─── - ───
   ⎝  ╲╱   E⋅I₂ ⎠    P     P 

Now as you said boundary condition is x = l, y = 0 then at x = l both real and img part should be zero giving you 2 equations with 2 variables.
But again I'm not quite sure if this is the right way to do it.

@codecov

This comment has been minimized.

Copy link

commented Jun 28, 2019

Codecov Report

Merging #17122 into master will increase coverage by 0.017%.
The diff coverage is 83.478%.

@@              Coverage Diff              @@
##            master    #17122       +/-   ##
=============================================
+ Coverage   74.495%   74.513%   +0.017%     
=============================================
  Files          623       624        +1     
  Lines       161220    161661      +441     
  Branches     37847     37926       +79     
=============================================
+ Hits        120102    120459      +357     
- Misses       35801     35851       +50     
- Partials      5317      5351       +34
@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 28, 2019

You need to declare the symbols as real/positive to get the sin/cos output from dsolve:

In [1]: x = Symbol('x')                                                                                                           

In [2]: f = Function('f')                                                                                                         

In [3]: a = Symbol('a')                                                                                                           

In [4]: dsolve(f(x).diff(x, 2) + a**2 * f(x))                                                                                     
Out[4]: 
           -ⅈ⋅a⋅x       ⅈ⋅a⋅x
f(x) = C₁⋅ℯ       + C₂⋅ℯ     

In [5]: a = Symbol('a', positive=True)                                                                                            

In [6]: dsolve(f(x).diff(x, 2) + a**2 * f(x))                                                                                     
Out[6]: f(x) = C₁⋅sin(a⋅x) + C₂⋅cos(a⋅x)

When working with complex numbers the exponential form is preferred. Both forms are valid for real or complex numbers though. Ideally dsolve would have an option to control which form you want.

@oscarbenjamin

This comment has been minimized.

Copy link
Contributor

commented Jun 28, 2019

The ics argument to dsolve can be used to get the value of the constants.

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jul 1, 2019

>>> x, y = symbols('x, y', positive=True)
>>> a, b = symbols('a, b', positive=True)
>>> solve(a*sin(x)-b*cos(x), x)
[-2*atan((a - sqrt(a**2 + b**2))/b)]

Shouldn't we simply be getting atan(b/a) for this?
Is there a way I could change its form?

ishanaj added some commits Jul 3, 2019

solve_slope_deflection now handles cases of trivial solutions where
deflection equation becomes zero. Also added some documentation
@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jul 10, 2019

Will be adding tests and some examples in the documentation next

ishanaj added some commits Jul 11, 2019

added tests for Column class. Also added an XFAIL for the problem in
determining the critical_load for pinned-fixed condition

@ishanaj ishanaj changed the title [WIP] Introduced a class Column which provides functionality for column buckling calculationss Introduced a class Column which provides functionality for column buckling calculations Jul 13, 2019

@ishanaj

This comment has been minimized.

Copy link
Contributor Author

commented Jul 13, 2019

The work on this PR has been completed. Its ready for review!

78000*y(x)
>>> c.solve_slope_deflection()
>>> c.deflection()
C1*sin(20*sqrt(195)*x/(sqrt(E)*sqrt(I)))

This comment has been minimized.

Copy link
@jashan498

jashan498 Jul 13, 2019

Contributor

Shouldn't the value of C1 get evaluated using boundary_conditions here?

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jul 14, 2019

Author Contributor

This is from the case where C1 and C2 comes out to be zero and hence even the deflection comes out to be zero. I have made a special case for it in the code.
Here C2 comes out to be zero and C1 remains unsolved.
I haven't found any resource of how to calculate the leftover value of C1, but I think maybe if the user provides some more boundary conditions we can solve it.
We could leave this for solve_for_unknowns() to do.

>>> c.boundary_conditions
{'deflection': [(0, 0), (3, 0)], 'slope': [(0, 0)]}
>>> c.moment()
78000*y(x)

This comment has been minimized.

Copy link
@jashan498

jashan498 Jul 13, 2019

Contributor

What represents y(x) here? Shouldn't the Function be well defined here?

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jul 14, 2019

Author Contributor

y(x) is the deflection (in the y-direction) which is a function of x.
Should I make it defl(x) instead?

This comment has been minimized.

Copy link
@jashan498

jashan498 Jul 18, 2019

Contributor

Nah, I think y(x) should be fine here. But is there a purpose of representing the moment in terms of deflection rather than in terms of x?
Why not just substitute the value of y(x) here with C1*sin(20*sqrt(195)*x/(sqrt(E)*sqrt(I)))?

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jul 21, 2019

Author Contributor

But the deflection is determined only after the user calls solve_slope_deflection().
So I guess for doing this we might need to create a variable to check whether deflection has been determined and then accordingly return the moment.
I am not sure whether to do this.

# differnetial equation of Column buckling
eq = E*I*y(x).diff(x, 2) + self._moment

defl_sol = dsolve(Eq(eq, 0), y(x)).args[1]

This comment has been minimized.

Copy link
@oscarbenjamin

oscarbenjamin Jul 13, 2019

Contributor

You should check that the solution has been returned in explicit form (that args[0] is just the function on the lhs of the equation).

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jul 14, 2019

Author Contributor

Does this mean to check whether args[0] comes out to be y(x)?

defl_sol = dsolve(Eq(eq, 0), y(x)).args[1]
slope_sol = defl_sol.diff(x)

constants = list(linsolve([defl_sol.subs(x, 0), slope_sol.subs(x, 0)], C1, C2).args[0])

This comment has been minimized.

Copy link
@oscarbenjamin

oscarbenjamin Jul 13, 2019

Contributor

You should use the ics argument to dsolve to find the integration constants rather than explicitly solving for them here.

This comment has been minimized.

Copy link
@ishanaj

ishanaj Jul 14, 2019

Author Contributor

Sure, will do it in the next commit

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.