Python Code for Mohr's Circle, - Stresse at a given angle theta, Principal stresses and their Planes, Maximum Shear Stress and its Plane

Hi,

You can use the Python to draw Mohr's Circle. In Structure Mechanics, taught to Civil Engineering and Mechanical

engineering students, one has to find the stresses at given oblique planes.

Principal stresses and maximum shear stresses are important to find, and also their 

corresponding planes. One can use the stress transformation equations or Mohr's circle method to 

find them. The following Python language code is prepared to be run on your IDE like Pycharm on your Computer.

One can also use online Jupyter's notebook to run the code. 

code: 

def mohcircle():
    import numpy as np
    import matplotlib.pyplot as plt
    import math

    σx = float(input('σx = '))
    σy = float(input('σy = '))
    τxy = float(input('τxy = '))
    u = input('stress unit = ')
    w = float(input("Angle( in degrees) of plane's axis from x axis(+ve counter clockwise), θ = "))
    θ = math.radians(w)
    R = np.sqrt(0.25 * (σx - σy) ** 2 + (τxy) ** 2)
    σavg = (σx + σy) / 2
    ψ = np.linspace(0, 2 * np.pi, 360)
    x = σavg + R * np.cos(ψ)
    y = R * (np.sin(ψ))
    φ1 = math.degrees(0.5 * math.atan(2 * τxy / (σx - σy)))
    φ2 = φ1 + 90
    σθ1 = σavg + R * np.cos(2 * np.radians(φ1) + 2 * θ)
    σθ2 = σavg + R * np.cos(2 * np.radians(φ1) + 2 * θ + np.pi)
    τθ = R * np.sin(2 * np.radians(φ1) + 2 * θ)
    print(f'''
       Radius, R = √(0.25*(σx-σy)^2 + τxy^2) 
               = √(0.25*({σx}-{σy})^2 + {τxy}^2)  ={R} {u}

       Average Stress,(acts at the Center of Mohr's Circle) 
               = σavg = (σx + σy)/2 = ({σx} + {σy})/2 = {σavg} {u}

       Principal Stresses
       σ1 = σavg + R = {σavg} + {R} = {σavg + R} {u}
       σ2 = σavg - R = {σavg} - {R} = {σavg - R} {u}
       Angle φ1 it makes with x-axis, 
       φ1 = 0.5(atan(2*τxy/(σx - σy)) = 0.5 * atan(2*{τxy}/({σx} - {σy})) = {φ1} degrees
       Angle σ2 makes with x-axis = φ2 = φ1 + 90 = {φ2} degrees

       Maximum Shear Stress = τmax = R = {R} {u}
       It occurs at, α = φ1 + 45 = {φ1 + 45} degrees

       Stresses at a plane with axis at θ anticlockwise from x axis, 
        σθ1 = σavg + R* Cos(2φ1 + 2θ) = {σavg} + {R}* Cos({2 * φ1 + 2 * θ})
           = {σθ1}, {u}
        σθ2 = σavg + R* Cos(2φ1 + 2θ + pi) = 
            {σθ2} {u}
        τθ = R*Sin(2*φ1 + 2*θ)  = {R * np.sin(2 * np.radians(φ1) + 2 * θ)} {u}

       ''')

    plt.plot(x, y)
    plt.plot([σavg - R - 10, σavg + R + 10], [0, 0], linestyle='--', color='black')
    plt.plot([σavg, σavg], [-R - 10, R + 10], linestyle='--', color='black')
    plt.plot([σx, σy], [τxy, -τxy], [σx, σx], [0, τxy], [σy, σy], [0, -τxy], linestyle='-', color='brown')
    plt.plot([σθ1, σθ2], [τθ, -τθ], [σθ1, σθ1], [0, τθ], [σθ2, σθ2], [0, -τθ], linestyle='--', color='red')
    plt.xlabel('σ (MPa)')
    plt.ylabel('τ (MPa)')
    plt.title("Mohr's Circle")
    plt.show()

mohcircle()
v = input('Exit? y/n ')
while v == "n":
    mohcircle()
    v = input('Exit? y/n ')
exit()

---- --------------

Thakns!

Solve Linear Equations in Three Variables - Using Python, Numpy and Sympy , With the Given Code

 Hi,

In the series of posts about Python for Civil Engineers, I have come with something from linear algebra.

While dealing with statics problems, such as finding unknown forces in space truss members, one has to deal with the system of linear equations in three variables. 

One can use the 'Sympy' and 'Numpy' libraries in Python can get a beautiful solution code for solving a system of linear equations in three variables. Computer programming is very much easy in Python, and you definitely need to try it. The code given below is separated into two parts: 1st part is the solution of the same problem using Sympy, and 2nd part gives solution code using Numpy library. 

The code is written specifically to solve the following three linear equations, but you can change the numbers accordingly.

x + 2y + 3z = 8

4x + 2y + 4z = 5

4x + 2y + 5z = 1

Simply change the coefficient for any other system. 

Code:

import numpy as np
import sympy as sp
from sympy import Eq, solve_linear_system, Matrix

#  Using sympy

x,y,z = sp.symbols('x y z')
row1 = [1, 2, 3, 8]
row2 = [4, 2, 4, 5]
row3 = [4, 2, 5, 1]
system = Matrix((row1, row2, row3))
print(system)
print(solve_linear_system(system,x,y,z))

print('___________________________________ ')

# Using numpy
row1 = np.array([1,2,3])
row2 = np.array([4,2,4])
row3 = np.array([4,2,5])
constants = np.array([8,5,1])
system = np.array((row1, row2, row3))
print(system)
xyz= np.linalg.solve(system, constants)
print(xyz)

Thanks!

Python Code to Find out the Dot Product of Vectors and Angle between them (Using end Point Co-ordinates too)

Hi, 

Sometimes you are given with the endpoints' co-ordinates of vectors, and sometimes you are directly given the vectors. For both conditions, the code given below will work. You have to input them as and when asked. End results will give the vectors, their magnitudes, their dot product and the angle between the two vectors. 

Python is a computer programming language that can be used to develip softwares, apps, and websites etc. The following python code provides a short solution to the Vector Algebra. 

Code: 

import numpy as np

import math

 a1 = input('For vector1, co-ordinates given? (y/n): ')

a2 = input('For vector2, co-ordinates given? (y/n): ')
u = input('unit: ')
if a1 == 'y':
    x11= float(input('enter co-ordinates of starting point, x11=  '))
    y11 = float(input('y11= '))
    z11 = float(input('z11 = '))
    V11 = np.array([x11, y11, z11])

    x21 = float( input('enter co-ordinates of end point, x21=  '))
    y21 = float(input('y21= '))
    z21 = float(input('z21 = '))
    V21 = np.array([x21, y21, z21])

    print(f'Starting point, V11 = {V11}')
    print(f'Ending point, V21 = {V21}')
    v1 = V21 - V11
    print(f'Vector1, V1 = V21 - V11 = {v1}')
    v1m = np.sqrt((x21-x11) ** 2 + (y21-y11) ** 2 + (z21-z11) ** 2)
    print(f" Magnitude of v1, v1m = {v1m} {u}")
else:
   print('enter the elements in vector1 ')
   x1 = float(input())
   y1 = float(input())
   z1 = float(input())
   v1 = [x1,y1,z1]
   print(f'v1 =({x1}i + {y1}j + {z1}k) {u}' '\n')
   v1m = np.sqrt(x1**2 + y1**2 + z1**2)
   print(f" Magnitude of v1, v1m = √(({x1})^2 + ({y1})^2 + ({z1})^2 = {v1m} {u}" '\n')

if a2 == 'y':
    x12= float( input('enter co-ordinates of starting point, x21 =  '))
    y12 = float(input('y12= '))
    z12 = float(input('z12 = '))
    V12 = np.array([x12, y12, z12])

    x22 = float( input('enter co-ordinates of end point, x22 =  '))
    y22 = float(input('y22= '))
    z22 = float(input('z22 = '))
    V22 = np.array([x22, y22, z22])
    print(f'Starting point, V12 = {V12}')
    print(f'Ending point, V22 = {V22}')
    v2 = V22 - V12
    print(f'Vector2, V2 = V22 - V12 = {v2}')
    v2m = np.sqrt((x22 - x12) ** 2 + (y22 - y12) ** 2 + (z22 - z12) ** 2)
    print(f" Magnitude of v2, v2m = {v2m} {u}")

else:
    print('enter the elements in the 2nd vector: ')
    x2 = float(input())
    y2 = float(input())
    z2  = float(input())
    v2 = [x2,y2,z2]
    print(f'v2 = ({x2}i + {y2}j + {z2}k) {u}' '\n')
    v2m = np.sqrt(x2**2 + y2**2 + z2**2)
    print(f" Magnitude of v2, v2m = √(({x2})^2 + ({y2})^2 + ({z2})^2 = {v2m} {u}" '\n')

d = float(np.dot(v1,v2))
e = d/(v1m*v2m)
a = math.acos(e)
k = math.degrees(a)

print(f'''Taking Dot product, v1.v2 = v1m * v2m * cosθ 

      => cosθ = v1.v2/(v1m*v2m) = {d}/({v1m}*{v2m}) = {e} 

      hence, Angle between the two vectors, θ = acos({e}) = {k}
''')

Python Code - Discharge in Rectangular Channel using Manning's Formula

 Hi,

Next in the series 'Python for Civil Engineers' is the code for Manning's formula for the rectangular channel.

One has to provide the needed inputs to get the discharge, in the given rectangular channel, as the output. The Manning's formula for discharge calculations can be used for various channels such as rectangular, trapezoidal etc. The Python code in this post is written for the rectangular channel, in particular. Inputs needed are Manning's n, slope s, bottom width b, and water depth h. Also, mention the units when asked. 

def discharge():
    n = float(input("n = Manning's value = "))
    S = float(input("S = Channel's slope = "))
    u = input("unit of width/height (m or ft)= ")
    b = float(input("b = bottom width = "))
    h = float(input("h = water depth = "))
    k = 1
    u1 = "Cumecs"
    if u.lower() == "ft":
        u1 = "cfs"
    if u.lower() == "ft":
        k = 1.49

    Q = (k/n) * (b * h) * (b * h / (2.0 * h + b)) ** (2.0 / 3.0) * S ** (0.5)
    Q = round(Q,3)
    print(f' Q = {Q} {u1}')
    return Q

print('''
    Manning's Formula for discharge in rectangular channel Q = (k * A * R^(2/3) * S^0.5) / n
        Q = Flow Rate, (in cfs, or m)
        k = 1.49 for discharge in cfs, and 1 for in cumecs
        A = Area of section = b*h , b = bottom width, h = water depth
        R = hydraulic Radius = Area/Wetted Perimeter = bh/(2h+b)
        n = Manning's Roughness coefficient
        ''')

print(f'Inputs needed to calculate the discharge Q')
discharge()

Thanks!

Vertical Equal Tangent Parabolic Curve - Solution Code in Python (Computer Program)

 Hi,

Now solving the vertical parabolic curve for the stations and elevations at the full stations is easier with the coding. A computer program that I have coded for you runs on Pycharm (Python). Pycharm is a software that is used for offline coding on your personal computer. If you have a good Internet connection, you can produce the code on online platforms such as Jupyter. 

 You have to put the value of the grades, length of the curve, station PVC and Elevation of PVC. It will calculate all other stations and elevations for you. 

Code: 

import math

g1 = input('Grade G1 (%) = ')
a = float(g1)/100
g2 = input('Grade G2 (%) = ')
b = float(g2)/100
Length = input('Length of the vertical curve = ')
L = float(Length)
u = input('unit = ')
A = float((b - a)/L)
print('Rate of change of grade = k = (G2-G1)/L =', A)
print('Let x be the distance from PVC to a point on the curve,'
      ' and Y be the respective elevation.')
c = float(input('Station(PVC) = Xo = ',))
d = float(input('Elevation at PVC = Yo = ',))
print(f"Equation of a parabolic equal tangent curve, "
      f"Y = Yo + G1x + kx^2/2 = {d} + {a}x + {A}x^2/2     ..eq1")
x = 100*((1+c//100) - c/100)
X = c + x,
print('1st Full Station is', X, 'at x =', x, u, 'from PVC')
Y = d + a*x + A*(x**2)/2
print('Putting the value of x in eq1, Elevation, Y= ', Y, u)
x = x + 100
while x < L:
    X = c + x
    Y = d + a * x + A * (x ** 2) / 2
    print('Similarly, next Full Station is', X, 'at x =', x, u, 'from PVC')
    print('Elevation, Y = ', Y, u)
    x = x + 100
X = c + L
print('End Station, EVC = PVC + L =', X, u, 'at x = L = ', L, u)
x = L
Y = d + a*x + A*(x**2)/2,
print('Put x in eq1, Elev EVC = ', Y, u)

Output Sample: 

Grade G1 (%) = 3
Grade G2 (%) = -2
Length of the vertical curve = 1000
unit = m
Rate of change of grade = k = (G2-G1)/L = -5e-05
Let x be the distance from PVC to a point on the curve, and Y be the respective elevation.
Station(PVC) = Xo = 340
Elevation at PVC = Yo = 100
Equation of a parabolic equal tangent curve, Y = Yo + G1x + kx^2/2 = 100.0 + 0.03x + -5e-05x^2/2     ..eq1
1st Full Station is (400.0,) at x = 60.00000000000001 m from PVC
Putting the value of x in eq1, Elevation, Y=  101.71 m
Similarly, next Full Station is 500.0 at x = 160.0 m from PVC
Elevation, Y =  104.16 m
Similarly, next Full Station is 600.0 at x = 260.0 m from PVC
Elevation, Y =  106.11 m
Similarly, next Full Station is 700.0 at x = 360.0 m from PVC
Elevation, Y =  107.56 m
Similarly, next Full Station is 800.0 at x = 460.0 m from PVC
Elevation, Y =  108.50999999999999 m
Similarly, next Full Station is 900.0 at x = 560.0 m from PVC
Elevation, Y =  108.96 m
Similarly, next Full Station is 1000.0 at x = 660.0 m from PVC
Elevation, Y =  108.91 m
Similarly, next Full Station is 1100.0 at x = 760.0 m from PVC
Elevation, Y =  108.36 m
Similarly, next Full Station is 1200.0 at x = 860.0 m from PVC
Elevation, Y =  107.31 m
Similarly, next Full Station is 1300.0 at x = 960.0 m from PVC
Elevation, Y =  105.76 m
End Station, EVC = PVC + L = 1340.0 m at x = L =  1000.0 m
Put x in eq1, Elev EVC =  (105.0,) m

Thanks!

Check Points of Buying a Flat/House

Hi,
Building a house is a mutual work of the builder, government and the buyers. A civil engineer plays its role in designing and constructing the house, and once he does his job, it is the buyer who has the last say if he wants to buy the house or not. If the house looks good, a buyer would love to buy the house, but there are other certain things that must be checked before buying the house from a builder. This blog post is about those important things that I have learned after buying my own flat but should have known them before.

 Check if the builder has got the project registered in the RERA website of the concerned state. For example, I bought a flat in Punjab, I checked the Punjab RERA website to see the project name under the registered projects list. It is very important that you do it because this will probably save you from any kind of possible cheating and fraud.

Things to check in the builder's RERA profile:

(A)
See if the project registration period has not expired. If it is expired, then you must avoid buying the flat.

(B)
See if there is a specified format of 'Agreement of Sale' in the uploaded documents section. You must insist to use this agreement format to finalise a deal with the builder. It has been seen that the builder asks you to sign the agreement as per his own terms and clauses which most of the times would be partial towards you. On the other hand, the format given by RERA will have the clauses which will give you all the rights that RERA has to offer to the buyers.
If you sign an agreement different from the one given by the RERA it is highly likely that builder would harass you in some way, like forfeiting your booking amount in case of cancellation of the deal, asking you for unfair maintenance charges, having partial interest rates for compensation in case of delay in possession of the flat etc.

(C)
If the project is shown as 'completed' in the RERA website then check for the completion and occupancy certificates. They would be uploaded in the website and hence you must check it. Remember that even if you get the registry papers of the flat you are still at risk. Until and unless you have the occupancy and completion certificates you are at the risk of losing your flat. You want to know why!? Please read about Campa Cola building in Mumbai.

(D)

 Don't give more than 10% of the money as the booking amount, this is the amount specified by RERA. If the builder asks you for more, please ask him to read the RERA rules.

(E)

 Ask the builder for the layout plan,  and as-built-plan and layout of the society.

(F)

 Check if the builder has all the clearance certificates and NOCs uploaded on the RERA account. You can see all these papers for free by signing up on the RERA website.

Once everything is okay, only then you must invest your money.

Thanks! 

All the basic Soil Water and Air Relations and Formulae, Three phase diagram II Soil Mechanics

Basic Soil and Water Relationships.
 Soil is considered a three-phase material, made of soil, water and air. The basic formulae used are written below. 

Sl No	Relationship in mass density	Relationship in unit weight
1	n = e/(1+e)	n = e/(1+e)
2	e = n/(1-n)	e = n/(1-n)
3	na = n ac	na = n ac
4	P= (G+Se)Pw/(1+e)	Y= (G+Se) Y w/(1+e)
5	Pd = GPw/(1+e)	Yd = G Yw/(1+e)
6	Psat = (G+e)Pw/(1+e)	Ysat = (G+e) Yw/(1+e)
7	P’ = (G-1)Pw/(1+e)	Y’ = (G-1) Y w/(1+e)
8	e = wG/s	e = wG/s
9	Pd = P/(1+w)	Yd = Y/(1+w)
10	Pd = (1-na)G Pw/(1+wG)	Yd = (1-na)G Y w/(1+wG)

Notations used: 
n = porosity, 
e = void ratio, 
na = air content,
P = Density,  Pd= dry density, 
Psat = saturated density, P' = submerged density
Y = Unit weight density, Yd = dry unit weight, 
Ysat = Saturated Unit weight, Y' = submerged unit weight
G = Specific Gravity
S = Saturation
w = water content

Thanks!