Continuum Mechanics: Create a Rich 2D Beam Solving System

About me

Name: Keshri Kumar Rushyam

University: Indian Institute of Technology, Kharagpur

Email: [email protected]

GitHub username: rushyam

Time-zone : UTC + 5:30

Personal Information

My name is Keshri Kumar Rushyam, a second-year undergraduate student at IIT Kharagpur, India. I am pursuing a degree in Ocean Engineering Naval Architecture. I am very much interested in learning Mechanics and Mathematics since my high school days. After joining the university, I got exposed towards the field of Computer Science. I have had taken courses like Mechanics, Mechanics of Solid, Ship Structural Analysis & Design, Probability and Stochastic processes, Advanced Engineering Mathematics, Numerical Methods, Programming and Data Structure, Fluid Mechanics, Marine Hydrodynamics.

Programming Details

I use Ubuntu 16.04 LTS on my workstation and Sublime Text (with some really cool extensions) as my primary text editor. Its multi-edit feature really assists me to search or edit anything quite easily. I am quite familiar with the concept of version control system and have been using Git and GitHub for quite some time and have enough knowledge to carry out the project successfully. I've been coding in C and C++ for over two years and in Python for under a year and am proficient in all three of them now.

I am very much interested in robotics and work for the Autonomous Underwater Vehicle group of my college. I have worked on and implemented Image enhancement Algorithms to enhance the underwater images.

My mathematics and physics background is sufficient for the project, with relevant courses.

My favorite feature of SymPy is pprint() it is very impressive and it makes mathematical expression easy to visualize.

>>> a, b = symbols('a, b')
>>> p = Polygon((0, 0), (a, 0), (a, b), (0,b))
>>> pprint(p.second_moment_of_area())
⎛   3   3     ⎞
⎜a⋅b   a ⋅b   ⎟
⎜────, ────, 0⎟
⎝ 12    12    ⎠

Contributions

I was introduced to Sympy by Shivam Vats in December 2017 and have been contributing ever since. Here is a list of all the merged and unmerged pull requests in chronological order.

Sympy

• (Merged) Implemented the Method for the second moment and product moment of area of a 2d polygon. #13939. Despite being my first pull request, it is my most important contribution to Sympy. I was introduced to the Sympy workflow while working on this PR and hence it probably took a much longer time than it should have.

• (Merged) Implemented the Method for the second moment and product moment of area of an ellipse. #14190

• (Unmerged) introductions of centroid() of ellipse #14439.for fixing the consistency issue.

• (Open) Implemented Shape class that can construct closed 2D Figure having methods area, centroid, the second moment of area and product moment of area #14434

• (Open) continuum_mechanics: statically indeterminate beam problem can be solved #14467

Project

Singularity functions are a popular tool for solving beam bending stress and deflection problems in mechanical design. In GSoC 2016 Singularity function was added to sympy and was used for solving fundamental 2D Beam problems. Current Beam module is unable to solve for statically indeterminate problems(the static equilibrium equations are insufficient for determining the internal forces and reactions on that structure). It solves for the Beams whose moment of area and elastic modulus is constant and does not have any cross-section related property(like neutral axis, normal stress, normal strain, etc).

Motivation

In beam bending problems in mechanics(statically determinate or indeterminate) the use of manual calculations is advocated which is very error-prone. Sometimes people use small simple functions and many assumptions to promote manual calculations with ease but by this way, the problems become an over-simplified model which cannot adequately portray physical systems.

The best way to tackle this issue would be by implementing a full-fledged computer algebra system supporting all the mathematical functionality of Singularity Functions. In this way, productivity would get significantly increased for solving problems that are more realistic and challenging, by giving more focus to the physical concepts rather than delegating with the horrendous mathematics involved.

Implementation details and Theory

I have proposed a solution in PR #14467 for solving statically indeterminate beam problems.

Current Beam class takes length, elastic modulus and moment of area as parameters and solves the problem whose elastic modulus and second moment of area are constant.

>>> b = Beam(L, E, I)

For solving the beam problems(statically determinate or indeterminate) whose E (elastic modulus) and I (second moment of area) varies along the length of the beam as shown in the diagram.

One possible solution can be:

We can express elastic modulus and second moment of area in terms of Singularity function. For the given diagram

E = E1*SingularityFunction(x, 0, 0) - E1*SingularityFunction(x, l, 0) + E2*SingularityFunction(x, l, 0) - E2*SingularityFunction(x, 2*l, 0)

I = I1*SingularityFunction(x, 0, 0) - I1*SingularityFunction(x, l, 0) + I2*SingularityFunction(x, l, 0) -I2*SingularityFunction(x, 2*l, 0)

In the present package modification in elastic_modulus setter method, second_moment setter method, Deflection method and Slope method will be done for solving variable elastic modulus and variable second moment of area problems.

>>> l, E1, E2, I1, I2 = symbols(l, E1, E2, I1, I2)
>>> b = Beam(l)
>>> b.elastic_modulus((E1, initial_1, order), (E2, initial_2, order))
>>> b.second_moment((I1, initial_1, order), (I2, initial_2, order))

• • E1 = elastic modulus of member AB
  • E2 = elastic modulus of member BC
  • ... so on.
• • I1 = second moment of area of member AB
  • I2 = second moment of area member BC
  • ... so on.
• • initial_1 = X coordinate of A
  • initial_2 = X coordinate of B
  • ... so on.
• • For constant , order=0
  • For linearly varying, order=1
  • For parabolically varying, order=2
  • ... so on.

Plotting

I will use Plotting module for plotting of Shear force, bending moment, slope, deflection of the beam. shear_force_plot, bending_moment_plot, slope_plot and deflection_plot Method will be implemented. Different color code will be used for shear force, bending moment, slope and deflection.

We will plot for 0<=x<=l. l is the length of the beam. deflection, slope method in Beam class does not give zero value for x > l (length of the beam). It will be improved.

>>> b = Beam(10, E, I)
>>> b.apply_load(R1, 0, -1)
>>> b.apply_load(M1, 0, -2)
>>> b.apply_load(10, 5, -1)
>>> b.bc_deflection = [(0, 0)]
>>> b.bc_slope = [(0, 0)]
>>> b.solve_for_reaction_loads(R1,M1)
>>> b.slope().subs(x, 20)
125/(E*I)                     #it should be zero

Arbitrary Cross Sections

In Geometry package Polygon, Ellipse and Circle have methods for calculating centroid, second moment of area and product moment of area. In pr #14434 I have implemented Shape class that can construct closed shape by boolean operations of parabola Ellipse and Circle with the Horizontal and the Vertical lines, and it has methods for finding centroid, second moment of area and product moment of area.

We can use Geometry package for the cross_section class.

>>> s1 = Section((geometry element), E1)
>>> s2 = Section((geometry element), E2)
>>> C1 = Cross_section((s1, s2), (1, 1))
>>> C2 = Cross_section((s1, s2), (1, -1))
>>> C3 = Cross_section(s1)

Geometry element can be Polygon, Ellipse, circle or closed constructed shape. E is the elastic modulus of that section.

section class will create an object which will have geometric property(geometry element) and Material property(elastic modulus (E)). Different sections will be added or subtracted(1 for addition and -1 for subtraction), and the cross-section of the beam will be generated using cross_section class. This implementation will help for solving composite beam and hollow beam problems.

• Example:

>>> p1 = Polygon((0, 0), (4, 0), (4, 3), (0, 3))
>>> p2 = Polygon((1, 1), (3, 1), (3, 2), (1, 2))
>>> p3 = Polygon((1, 3), (3, 3), (3, 4), (1, 4))
>>> s1 = Section(p1, E1)
>>> s2 = Section(p2, E1)
>>> s3 = Section(p3, E2)
>>> C1 = Cross_section((s1, s2), (1, -1))
>>> C2 = Cross_section((s1, s3), (1, 1))
>>> C3 = Cross_section((s1, s2, s3), (1, -1, 1))

If we are joining section A and section B (A and B are gemetric elements): if A and B touching each other then we can join and If they are overlapping or completely separated then it will return an error.

If we are trimming section A from section B, then it will be checked whether section A is completely inside the section B or not and whether the elastic modulus of section A is same as that of section B or not. If both are true then we can trim else it will raise an error

Strain at any point P is the sum of strain due to normal force(N) and due to the bending moment in y and z-direction. K_y and K_z are the radius of curvatures along the y-axis and z-axis respectively.

stress at any point P is equal to elastic modulus times strain.

# basic formula 1

workflow of the class will be like this:

Superposition of different sections will be used for finding the area, centroid and second moment of area of the cross-section.

• cross_section class will have methods:

  • area

    >>> C1.area
    area_of_the_shape
    

  • centroid

    >>> C1.centroid
    (x, y)
    

  • second_moment_of_area

    It will return second moment and product moment of area about p and if p is None then it will return about the centroid.

    >>> p = Point(x1, y1)
    >>> C1.second_moment_of_area(p)
    I_xx, I_yy, I_xy
    

  • apply_moment

    >>> N = CoordSys3D('N')
    >>> M = M_y*N.j + M_z*N.k
    >>> C1.apply_moment(M)
    

    Parameter: Moment vector

  • radius_of_curvature It will return radius of curvature vector.

    >>> C1.radius_of_curvature
    k_y*N.j + k_z*N.k
    

  • neutral_axis

    It will return the equation of Neutral axis or line. The line on which stress and strain are zero.

    >>> C1.neutral_axis
    a*y + b*z +c = 0
    

  • normal_stress

    It will return the normal stress at any point p on the cross-section.

    >>> p = Point2D(x, y)
    >>> C1.normal_stress(p)
    value
    

  • normal_strain

    It will return the normal strain at any point p on the cross-section.

    >>> p = Point2D(x, y)
    >>> C1.normal_strain(p)
    value
    

  • max_stress

    It will return the value of maximum stress.

    >>> C1.max_stress
    value
    

  • max_strain

    It will return the value of maximum strain.

    >>> C1.max_strain
    value
    

  • location_of_max_stress

    It will either return a 2D line segment or a 2D point depending upon the cross-section and moment.

    >>> C1.location_of_max_stress
    Line2D(Point2D(x1, y1), Point2D(x2, y2))
    >>> C4.loaction_of_max_stress
    Point2D(x1, y1)
    

Sphere

3-Dimensional Geometry The current geometry module restricts the functionality of it's entities to 2-dimensional space.I would like to extend the 3-dimensional shapes and will start with the sphere

>>> r = symbols('r', positive=True) 
>>> s = Sphere(Point3D(0, 0, 0), r)

It will have method like equation, tangent_plane at any point p, normal_plane at any point p, is_tangent will check whether a plane is a tangent plane to sphere or not by using the property perpendicular distance of plane from centre is equal to radius, intersection of line and plane with sphere. The intersection of the sphere with the line will return Point or Points and with the Plane will return point or circle.

Sphere is the basic shape that we always use in rotational mechanics. So I will implement methods for finding mass, the moment of inertia about x-axis, y-axis, and z-axis for the variable density of sphere with respect to distance from the center of the sphere.

Density can be represented as Singularity function. For a hollow sphere with the inner radius and a and outer radius b and having uniform density rho. r is the distance from the center of the sphere.

d = rho*SingularityFunction(r, a, 0) - rho*SingularityFunction(r, b, 0) 

mass = integrate(d4Pi*r**2,(r, a, b)

>>> s.mass(d)
((rho*4*Pi*b**3)/3) - ((rho*4*Pi*a**3)/3)

moment of inertia is the second moment of mass with respect to distance from an axis

Moment of inertia of inertia about any axis can be found using parallel axis theorem. 1st moment of inertia about centroid will be found using integration and then using parallel axis theorem we can find the moment of inertia about any axis.

>>> p = Point3D(0, 0, 0)
>>> s.moment_of_inertia(d, p)
(-8*pi*a**5*rho/15 + 8*pi*b**5*rho/15, -8*pi*a**5*rho/15 + 8*pi*b**5*rho/15, -8*pi*a**5*rho/15 + 8*pi*b**5*rho/15)

Timeline (Tentative)

Any Plans/Commitment (During GSoC):-

• I am having End Semester Examination from 19 April to 27 April. But this would not be an issue for having Community Bonding. My summer vacations start by 28 April. I would be able to devote 40 - 50 hours a week during the project since I have no other significant commitments other than GSoC.

I will be free for the whole summer break. I can give a definite assurance that this vacation will not affect the progress of my project. As this project will form the core of all my working and learning throughout the Summer.

• My college reopens on July 17, but for the first month, I would still be able to devote around 40 hours, since there would be no tests/exams. I can continue with same promptness.

Community Bonding Period - Week 1 (22 April - 29 May)

Goal: Community Bonding

• I will have end term examination till last week of April, after that, I will be productive and can start right Away. My principal focus will be to discuss more the project with mentors and also to figure out how to make more improvement in the 2d beam class.

• This will be a great time to get to know everybody and the fellow students.

• I will try to get all the PR merged related to my project.

Week 2-3

Goal: Improvement in Beam class will be done for solving variable elastic modulus and variable second moment of inertia will be done.

*Modification in elastic_modulus setter method, second_moment setter method, Deflection method and Slope method will be done

• I would also check if any other improvements are required else I can proceed towards the second phase of my project.

• I will be writing the tests along with the development. I will also prefer to document the newly added functionalities.

Week 3-Week 4

Goal: Plotting

• shear_force_plot, bending_moment_plot, slope_plot and deflection_plot Method will be implemented.

• I will write tests as well the documentation of the newly implemented function.

Week 5 - Week 6

Goal: Start with Implementing Cross-Section Class.

• This week I would start with implementing Beam Cross-section Class.

• Following that, I would start with creating the initial phase of section class. I will then implement the cross-section class`. The condition check of joining and trimming of more than one sections (as mentioned above) will be done

• I will be writing tests along with all the function.

Week 7 - Week 9

Goal: Complete the Implementation part of Cross-section class.

• I would continue the work of the previous weeks.
• I will work on implementing area, centroid, second_moment_o_area, apply_moment, neutral_axis, normal_stress, normal_strain and location_of_max_stress.
• I will write tests for these implemented methods. I will fix the bugs and issues.
• By the end of this time, I will have successfully implemented a Cross-section class

Week 10 - Week 11

Goal: Buffer Period and Wrap Up

This week would be a compensation for everything that might delay the completion. This will mostly involve :

• Preparing documentation of examples and tutorials for implemented class. A user will use software more if its documentation is having plenty of examples and tutorials depicting the uses of every function.
• Finish pending tasks for this week are the buffer.
• I will be cleaning up the code and refactoring and fix bugs.

Week 12 - Week 13

• If I would finish all the work mentioned above in time then I will start implementing sphere class.
• I will write tests and documentation for these implemented methods.

Post GSoC :

• I will be fixing bugs and reviewing pull requests in Sympy and would look after the work done by me.
• I will implement classes like cone, cuboid in geometry package it will have methods like center of mass, volume, the moment of inertia, etc.

Reference

cross-section

Application Sampad Kumar Saha: Singularity Functions

problems of 2d beam