GSoC 2018 Application Subhash Saurabh: Improving 2D Beams and Implementing 3D beams - sympy/sympy GitHub Wiki
Name : Subhash Saurabh
University : Indian Institute of Technology, Guwahati
Major: Civil Engineering
Email: [email protected]
Github: Subhash-Saurabh
Timezone: IST (UTC +5:30)
My name is Subhash Saurabh a third year undergraduate at IIT Guwahati,India. I was very much interested in learning Mathematics since my high schools. After joining the university, I got exposed towards the field of Computer Science and I was very much interested in it. Even though my major is in Civil Engineering I have therefore taken many courses related to Computer Science and Mathematics such as Probability and Statistics , Real and Complex Analysis, Data Structures and Algorithms and Machine Learning.I have also taken courses such as Solid Mechanics, Structural Engineering (I & II) that will help me in implementing the project.
I am currently using Ubuntu 14.04 LTS and Atom as my primary text editor as it is open source, there are many plugins available and most important it has autocomplete feature which helps to write code faster.
I have 2+ year of experience with python and I have used it for Data Science(pandas, numpy, sklearn) , web scrapping (BeautifulSoup), and for many other programming related work. One of the projects that I did related to programming is Genome Assembling (designing an algorithm to find the genome from its read). I have a good level of understanding of git and I think it will be sufficient to complete this project.
One of the most interesting feature of sympy is it's ability to solve indefinite integrals. For example
>>> Integral(x**4*sin(3*x),x).doit()
-x**4*cos(3*x)/3 + 4*x**3*sin(3*x)/9 + 4*x**2*cos(3*x)/9 - 8*x*sin(3*x)/27 - 8*cos(3*x)/81
I got introduced to sympy in Nov 2017. Since then I have made several pull request to the sympy repository. Here is the list of my pull request.
- Added simplification of mul having args as units, prefixes and other arguments
- Corrected latex printing of units containing underscores
- Faster implementation of sample function in drv_types
- Made changes to resolve errors with mul and div in prefixes
- Fixed issue of wrong prefixes under multiplication
- Made changes to ensure no precision loss of Float on computing Jordan_form
- Removed backslash-escaping of square brackets
- Added abbreviated printing for units
- Resolved error during substitution of Quantities
- Fixed a typo in mechanics documentation
- Added angular quantities printing in latex and pretty
- User defined values in mul_symbol in latex
- Implemented presentation mathML printer
- Fixed issue of simplify error in sinc
- Fixed Documentation of gamma function
- Fix missing parentheses in LaTex and pprint output
- Added support of continuum_mechanics with units
- Added ascii printing for symbols
- Fix Problem with substitution and paranthesis handling (Sympy's Turn)
In GSOC 2016 Continuum mechanics package was added to sympy. In this project I plan to improve this package. In PR#14467 solution for indeterminate beam has been proposed.Therefore it is not included in the proposal. However it will be made sure that beam can be solved for all indeterminate cases before GSOC starts.
Methods like plot_shear_force
, plot_bending_moment
, plot_deflection
, plot_slope
will be implemented which will return the plots of shear force, bending moment, deflection and slope respectively.
Implementation
For plot_shear_force
- First compute shear_force using the method
shear_force
. - Using offset value of Singularity Function compute equation of shear_force in piecewise form.
- Using plot class of plotting module plot the shear_force.
Other plotting method will also be implemented in the same way.
Currently the methods shear_force
, bending_moment
, slope
and deflection
returns only the equation. Modification will be done to these function such that they can also take position of beam as argument and return the respective values of equation at that point of beam.
The expected behaviour would be
>>> b.shear_force()
-3*SingularityFunction(x, 0, 0) + 6*SingularityFunction(x, 2, 1) - 9*SingularityFunction(x, 4, 0)
>>> b.shear_force(2)
-3
Similar behavior is expected for bending_moment, slope and deflection.
Theory
An influence line is a diagram whose ordinates, which are plotted as a function of distance along the span, give the value of an internal force at a particular point in a structure as a unit load move across the structure. It represents the variation of either the reaction, shear or moment at a specific point in a member as a unit concentrated force moves over the member.
Influence lines can be plotted by deriving the general equations based on singularity function.
- First the given moving load is applied at any arbitrary position on the beam (say x0).
- Then shear force and bending moment can be calculated in terms of x0 and singularity function.
- Then we substitute value of x in singularity function by the point where ild has to be calculated.
For example
Implementation
- A method
apply_moving_load
will be implemented that will apply a moving load on the beam. This method will take two argument magnitude of force applied and order (if order is -2,-1) or magnitude of force ,length on which force is applied and order (if order is 0,1,2). - A method
ild
will be implemented which takes the argument as position of beam and return the ILD for shear force and bending moment in form of a list. - Methods like
max_shear_force
andmax_bending_moment
will be implemented that will return the maximum shear force and maximum bending moment at a point in case of moving load. - Another method
plot_ild
will be implemented which takes the argument as position of beam and plot the ILD for shear force and bending moment.
Example
Let us consider a beam of length 10 meter. A unit moving point load is applied on the beam. Find the ild at 6 meter of the beam.
The expected behaviour would be
>>> b = Beam(10,'E','I')
>>> b.apply_load(R1,0,-1)
>>> b.apply_load(R2,0,-1)
>>> b.apply_moving_load(-1,-1)
>>> b.bc_deflection = [(0,0),(10,0)]
>>> b.ild(6)
[-1*SingularityFunction(6, x, 0) - 0.1*x + 1, -1*SingularityFunction(6, x, 1) - 0.6*x + 6 ]
>>> b.max_shear_force()
-0.6
CrossSection
class will be implemented that can compute various parameters that are needed for 3D analysis of beam.
-
The shear area is a calculated reduction of a cross-sectional area. By using this value, we can consider the shear deformation when determining the internal forces. The shear areas calculated will only be used to determine the internal forces.
-
The static moment of area simply measures the distribution of a beam sections's area relative to an axis. It is calculated by taking the summation of all areas, multiplied by its distance from a particular axis.
-
Bending centre is the point on the cross-section through which if normal force is applied then only a uniform axial displacement is induced over the cross-section.
-
The principal axes of the cross section are defined as the rectangular axes for which the product of inertia vanishes(Iyz = 0). Moment about principle axes can be related to moment about original axis as
Calculating Moment about principle axes taking bending center as origin is important for 3D analysis of beams.
- Method
area
will be implemented that can compute the area of cross-section. - Method
shear_area
will be implemented that will compute the shear area of the cross-section(Ay, Az). - Method
static_moment
will be implemented that will compute the static moment of the cross-sectional area(Sy, Sz). - Method
moment_of_inertia
will be implemented that will calculate the three moment of inertia (Iyy, Iyz, Izz). - Method
bending_center
will be implemented that will compute the bending center of the beam (yB, zB). - Method
moment_about_principle_axis
will be implemented that will compute the moment_of_inertia taking bending center as origin and principle axes orientation.
Expected output
>>> p1,p2,p3,p4 = [(0, 0), (2, 0), (0, 1), (2, 1)]
>>> c = CrossSection(Polygon(p1,p2,p3,p4))
>>> c.area()
2
>>> c.shear_area()
[1.66666666666667, 1.66666666666667]
>>> c.static_moment()
[1, 2]
>>> c.moment_of_inertia()
[2.66666666666667, 1, 0.66666666666667]
>>> c.bending_center()
[1, 0.5]
>>> c.moment_about_principle_axis()
[0.66666666666667, 0.166666666666667]
Also a method set_second_moment
will be implemented in beam class that takes a list as argument and sets the second moment of inertia of the 2D beam.
The expected behaviour will be
>>> p1,p2,p3,p4 = [(0, 0), (3, 0), (0, 4), (3, 4)]
>>> b = Beam(4,E,I)
>>> b.set_second_moment(CrossSection(Polygon(p1,p2,p3,p4)))
>>> b.second_moment
16
The 3D beam class can be instantiated by passing five arguments length of beam, Cross-section of beam (A cross-section object), Elasticity of Beam, Shear modulus of beam and a fifth optional argument specifying whether the Beam is to be modeled by Timoshenko or Bernoulli-Euler theory. The default value will be set to 'Bernoulli-Euler'.
>>> E, G = symbols('E, G')
>>> p1, p2, p3, p4 = [(0, 0), (2, 0), (0, 1), (2, 1)]
>>> c = CrossSection(Polygon(p1, p2, p3, p4))
>>> b = 3DBeam(10, c, E, G, 'Timoshenko')
Theory
When we take (x, y, z) - coordinate system as principal axes coordinate system with origin at the bending center the Benoulli-Euler and Timoshenko beam reduce into following subsystem.
- Timoshenko beams
Axial deformation is governed by
whereas bending deformation in the y-direction is defined by the coupled equations
where the double index yy on the bending moment of inertia has been replaces by a single index y in order to indicate that the principal-axes coordinates are utilised.
Similarly, the flexural deformations in the z-directions are determined by
- Bernoulli-Euler Beam Bernoulli-Euler beam is governed by the following equations
- Method
boundary_condition
will be implemented that will allow the user to specify the boundary conditions. - Method
solve_beam
will be implemented that will solve the beam and provide the unknown slope ,deflection, force and moment at any point.
- We will first shift the axes to principle axes and origin to bending center.
- Using the new axes our equations will reduce to three subsystems.
- Each subsystem will then be solved using the boundary_conditions.
- Method
force
will be implemented which will take the position on beam as argument and return the five dimensional vector(Qx, Qy, Qz, My, Mz). - Method
displacement
will be implemented which will take the position on beam as argument and return the five dimensional vector(ex, ey, ez, slopey, slopez).
My end-sem exams will get over on 5th May and my college will reopen on 24 July.Therefore I will be able to contribute 40-45 hours weekly to my project. Even after my college opens work load will be fairly less for initial weeks and therefore I will be able to complete my project.
- During this period I will get myself familiar with the plotting module, Integration module and other modules that that will be required in my project . I will make sure that all unmerged PR (such as proposed solution of indeterminate beams) relating to my project gets merged with addition/changes if required. During this period I will also discuss the project with the mentors.
- Implement the plotting of shear_force, bending_moment, slope, deflection.
- PR will be submitted at the end of week 1 that will contain the mentioned functionality, tests and documentation.
- It will be made sure that plotting works for both determinate and indeterminate beams.
-
shear_force
,bending_moment
,deflection
,slope
methods will be modified as mentioned in the plan. - Appropriate test and documentation will be added.
- In this week
apply_moving_load
method will be implemented that will be necessary for ILD. - In this week itself work will also start on implementing
ild
method.
- In these two week
ild
,max_shear_force
andmax_bending_moment
method will be implemented. - The
ild
method will be able to compute ild for all loading conditions(order = -2,-1,0,1,2). - It will also work for fixed loading + moving load condition.
- Proper test and documentation will also be added.
- In this week
plot_ild
method will be implemented that can plot the ild. - All relevant test and documentation will be added.
- In these two weeks
CrossSection
class will be implemented. - All methods of of CrossSection class
area
,shear_area
,static_moment
,moment_of_inertia
,bending_center
,rotation_principle_axis
,moment_about_principle_axis
will be implemented. - Also
set_bending_moment
method will be implemented. - All tests and documentation will also be added.
- In this weeks
3DBeam
Class will be implemented. - Method
boundary_condition
will also be implemented. - All the tests and documentation will be added.
- All the remaining methods
solve_beam()
,force()
,displacement()
will be implemented. - All the tests and documentation will also be added.
- Any remaining work will be completed.
Final Submission
During Gsoc only beams without torsion will be implemented. After GSOC I plan to improve it to include torsion in beams.Also I plan to create a sympy subpackage dealing with truss, frame and other structures. This project is related to my academics so I will be regularly contributing to it.
- GSoC 2016 Application Sampad Kumar Saha: Singularity Functions : For design of proposal
- Influence line : ILD
- http://www.uap-bd.edu/ce/anam/Anam_files/Structural%20Engineering%20I.pdf : ILD
- List of second moments of area : Moment of inertia
- http://old.ceric.net/old_data/trend/data/new/FEA%20KOREA/ShearArea.pdf : Shear Area
- http://homes.civil.aau.dk/jc/FemteSemester/Beams3D.pdf : 3D beams and Cross-Section