I intend to design the foundation of an elevated water tank, consisting of a raft with the shape shown in the image below. To do so, I need to determine the raft diameter Dr by solving the following inequality:
NT∙Dr ≥ 8∙M (where NT and M are inputs in my code)
Based on my hand calculations, after taking into account all the necessary parameters to evaluate NT, the final equation to be solved is a cubic polynomial with the following expression:
5.694·Dr3 - 0.942·Dr2 +2540.505·Dr - 73612.480 = 0
I can easily solve this cubic polynominal equation by appling the following code which return the real solution Dr =17.360 that match my hand calculation.
import sys
from scipy.optimize import fsolve
def equation(Dr):
return 5.694 * Dr**3 - 0.942 * Dr**2 + 2540.505 * Dr - 73612.480
# initial guess for the solver
initial_guess = 0.0
root = fsolve(equation, initial_guess)
OUT = root[0]
However, I would like to automatically determinate the polynomial coefficients instead of computing them manually. This would allow me to reuse the code later to design another raft based on different inputs, without having to manually recalculate these coefficients, for that I appied the code below which return a wrong results as you can see:
Cubic polynominal equation:
6.676·Dr3 - 0.942·Dr2 +2540.505·Dr - 73612.480 = 0Real solution:
Dr = 16.744 m
import sys
import os
sys.path.append(r'C:\Users\nono\AppData\Local\Programs\Python\Python38\Lib\site-packages')
import numpy as np
import math
def vertical_load(H, a, b, c, di, de, gb, gs, P0, Dr):
# calculating the entire vertical force at the base
# H is embedment depth
# de: shaft outer diameter
# di: shaft inner diameter
# gb is the concrete density = 2.5T/m3
# gs is the soil density = 1.9T/m3
# P0 vertical load from the superstructure
# Dr raft diameter to find
h = H / 4
d1 = de + 2*a
d2 = d1 + 2*b
d3 = Dr - 2*c
# raft area calculation
s1 = de**2 - di**2
s2 = Dr**2
s3 = (d3**2 - d2**2) / 2
s4 = d2**2
s5 = (d2**2 - d1**2) / 2
s6 = d1**2
# raft weight calculation
Pr = gb * h * math.pi / 4 * (s1 + s2 + s3 + s4 + s5 + s6)
# soil area calculation
t1 = Dr**2 - d3**2
t2 = (d3**2 - d2**2) / 2
t3 = Dr**2 - d2**2
t4 = (d2**2 - d1**2) / 2
t5 = Dr**2 - de**2
# soil weight
Pt = gs * h * math.pi / 4 * (t1 + t2 + t3 + t4 + t5)
# Total vertical load
NT = P0 + Pr + Pt
return NT
def Equi_inequality(Dr, H, a, b, c, di, de, gb, gs, P0, M):
NT = vertical_load(H, a, b, c, di, de, gb, gs, P0, Dr)
return NT * Dr - 8 * M
def polynom(H, a, b, c, di, de, gb, gs, P0, M):
Dr_vals = np.array([0.0, 1.0, 2.0, 3.0])
F_vals = np.array([
Equi_inequality(Dr, H, a, b, c, di, de, gb, gs, P0, M)
for Dr in Dr_vals
])
A = np.array([
[0, 0, 0, 1],
[1, 1, 1, 1],
[8, 4, 2, 1],
[27, 9, 3, 1]
])
coeffs = np.linalg.solve(A, F_vals)
return coeffs
def solve_Dr(H, a, b, c, di, de, gb, gs, P0, M):
a3, a2, a1, a0 = polynom(
H, a, b, c, di, de, gb, gs, P0, M
)
roots = np.roots([a3, a2, a1, a0])
real_roots = [
r.real for r in roots
if abs(r.imag) < 1e-8 and r.real > 0
]
if not real_roots:
raise ValueError("no solution found")
return min(real_roots), (a3, a2, a1, a0)
# Inputs
H = 4.00
a = 1.05
b = 2.00
c = 1.00
di = 6.80
de = 7.60
gb = 2.5
gs = 1.9
P0 = 2492.52
M = 9201.56
Dr, coeffs = solve_Dr(H, a, b, c, di, de, gb, gs, P0, M)
print(f"{coeffs[0]:.3f}·Dr3 {coeffs[1]:+.3f}·Dr2 "
f"{coeffs[2]:+.3f}·Dr {coeffs[3]:+.3f}")
print(f"Dr = {Dr:.3f} m")
OUT = Dr, coeffs = solve_Dr(H, a, b, c, di, de, gb, gs, P0, M)
How to fix this code to return the right solution?
Remark: I'm using Cpython3 engine within Dynamo for Revit to execute the code.
Edit:
According to my hand calculations, and in order to show how I arrived at the final expression of the polynomial equation, I'd like to highlight the following key points:
NT∙Dr ≥ 8∙M:Equilibrum inequality to be checked whereMis an input andNTis the Totale vertical load defined asNT = P0 + Pr + Ptwhere:P0is an input value in the codePr: raft weight computed according to the geometry shown in the imagePt: soil weight computed according to the geometry shown in the image.
By taking into account all the input data and applying the equilibrium inequality, we obtain the cubic polynomial equation derived from my hand calculations,shown in this image: equilibrium inequality
1 Answer 1
In all likelihood, your code is right and your hand calculations are wrong. In the calculation of Pr, the order-2 term should have a non-zero coefficient due to s3 but it does not.
Your s3 values are negative which doesn't make sense for areas. This is because you use non-physical initial values for Dr in your polynomial that, while in theory more numerically stable, produce a less accurate result. If you pass in an estimate for this value that's closer to the real thing, your inequality error decreases by something like three orders of magnitude.
I will also offer some techniques that will improve the quality of the solver and make consistency checking more feasible.
Don't sys.path.append; that's a sign of a broken environment and should not be hacked in the script.
Don't import math; use Numpy instead.
Add type hints: most of your parameters are scalars, and Dr may be either a scalar or a Numpy array. Add unit descriptions to all of your parameters.
Treat Dr as being vectorised, and do not write a loop comprehension. All four values of Dr will be processed at the same time, and so stepping through will have improved visibility for the intermediate variables like Pr, etc.
If you move most of the system parameters to an immutable class, then it will be easier to pass them around (rather than as a long argument list).
Don't write A manually. That's just a broadcast outer-power. Likewise, do not write Dr_vals manually; that's just a range.
Use a mask rather than a loop comprehension for real_roots.
Do not write # Inputs etc. in the global namespace. Those symbols are visible to all of your functions and so create shadowing problems. Put them in a function.
I implore you to use a real debugger and step through, checking every one of your quantities against your hand calculations (that you still have not shown).
All together, this could look like
import typing
import numpy as np
# This is an approximation; typing in Numpy is rough
type TNum = float | np.ndarray
class System(typing.NamedTuple):
H: float # embedment depth (m)
a: float # ? (m)
b: float # ? (m)
c: float # ? (m)
de: float # shaft outer diameter (m)
di: float # shaft inner diameter (m)
gb: float # concrete density (T/m3)
gs: float # soil density (T/m3)
P0: float # vertical load from the superstructure
M: float # ?
# Dr is the raft diameter to find
def vertical_load[T: TNum](self, Dr: T) -> T:
"""calculate the entire vertical force at the base"""
h = self.H/4
d1 = self.de + 2*self.a
d2 = d1 + 2*self.b
d3 = Dr - 2*self.c
# raft areas
s1 = self.de**2 - self.di**2
s2 = Dr**2
s3 = (d3**2 - d2**2)/2
s4 = d2**2
s5 = (d2**2 - d1**2)/2
s6 = d1**2
# raft weight
mu = self.gb * h * np.pi/4
Pr_old = mu*(s1 + s2 + s3 + s4 + s5 + s6)
# p2 = mu*( 1 + 0 ) # artificially remove the second term
# should be:
p2 = mu*( 1 + 0.5*1 )
p1 = mu*( 0.5*-4*self.c )
p0 = mu*(s1 + 0.5*(4*self.c**2 - d2**2) + s4 + s5 + s6)
Pr = p2*Dr**2 + p1*Dr + p0
if not np.isclose(p2, 1.963, atol=0, rtol=1e-3):
print(f'Warning: p2 {p2} != 1.963')
if not np.isclose(p1, -3.927, atol=0, rtol=1e-3):
print(f'Warning: p1 {p1} != -3.927')
if not np.isclose(p0, 487.448, atol=0, rtol=1e-3):
print(f'Warning: p0 {p0} != 487.448')
if not np.allclose(Pr, Pr_old, atol=0, rtol=1e-3):
print(f'Warning: Pr diverge, {Pr}, {Pr_old}')
# soil areas
t1 = Dr**2 - d3**2
t2 = (d3**2 - d2**2)/2
t3 = Dr**2 - d2**2
t4 = (d2**2 - d1**2)/2
t5 = Dr**2 - self.de**2
# soil weight
Pt = self.gs * h * np.pi/4 * (t1 + t2 + t3 + t4 + t5)
# Total vertical load
NT = self.P0 + Pr + Pt
return NT
def equi_inequality[T: TNum](self, Dr: T) -> T:
NT = self.vertical_load(Dr)
return NT*Dr - 8*self.M
def polynom(self) -> np.ndarray:
Dr_min = self.de + 2*(self.a + self.b + self.c)
Dr = np.linspace(start=Dr_min, stop=1.1*Dr_min, num=4)
A = Dr[:, np.newaxis]**np.arange(3, -1, -1)
F = self.equi_inequality(Dr)
return np.linalg.solve(A, F)
def solve_Dr(self) -> tuple[
float, # lowest real root of Dr
np.ndarray, # polynomial coefficients
]:
a = self.polynom()
roots = np.roots(a)
mask = (abs(roots.imag) < 1e-8) & (roots.real > 0)
if not mask.any():
raise ValueError('no solution found')
return roots[mask].real.min(), a
def main() -> None:
system = System(
H=4., a=1.05, b=2., c=1.,
di=6.8, de=7.6, gb=2.5, gs=1.9,
P0=2492.52, M=9201.56,
)
Dr, (m, n, o, p) = system.solve_Dr()
print(f'{m:.3f}·Dr3 {n:+.3f}·Dr2 {o:+.3f}·Dr {p:+.3f}')
print(f'Dr = {Dr:.3f} m')
print(f'NT = {system.vertical_load(Dr):.3f}')
print(f'ineq = {system.equi_inequality(Dr):.1e}')
if __name__ == '__main__':
main()
The correct output should be
Warning: p2 2.945243112740431 != 1.963
6.676·Dr3 -0.942·Dr2 +2540.505·Dr -73612.480
Dr = 16.744 m
Warning: p2 2.945243112740431 != 1.963
NT = 4396.369
Warning: p2 2.945243112740431 != 1.963
ineq = 1.5e-11
Using the artificially zeroed term, the incorrect but consistent output is
Warning: Pr diverge, [ 886.99923331 953.16466525 1024.37409873 1100.62753375], [1116.8165534 1218.28999611 1327.32944113 1443.93488845]
5.694·Dr3 -0.942·Dr2 +2540.505·Dr -73612.480
Dr = 17.360 m
Warning: Pr diverge, 1011.0370702193181, 1306.9188505015122
NT = 4240.258
Warning: Pr diverge, 1011.0370702193181, 1306.9188505015122
ineq = -2.9e-11
After your other modifications, and taking your new diagram into account:
3D diagram
you should note the following changes:
TNummust be aTypeVarof aUnion, not aUnionembedment depthwas misspelled- Don't write
None; keepnp.newaxis - With the radius change, your
R_minwas incorrect - Don't remove
main; return a dictionary from it - Use
pprintto print the dictionary - Using Python 3.8 is a bad idea; you need to upgrade. I show the modifications necessary to make Python 3.8 compatible, but this should only be temporary.
import pprint
import typing
import numpy as np
# This is an approximation; typing in Numpy is rough
TNum = typing.TypeVar('TNum', bound=typing.Union[float, np.ndarray])
class System(typing.NamedTuple):
H: float # embedment depth (m)
a: float # extension distance needed for radius r1 (m)
b: float # extension distance needed for radius r2 (m)
c: float # extension distance needed for radius r3 (m)
re: float # shaft outer radius (m)
ri: float # shaft inner radius (m)
gb: float # concrete density (T/m3)
gs: float # soil density (T/m3)
P0: float # vertical load from the superstructure (T)
M: float # overturning bending moment (T.m)
# R is the raft radius to find
def vertical_load(self, R: TNum) -> TNum:
"""Calculate the total vertical force at the base"""
h = self.H/4
r1 = self.re + self.a
r2 = r1 + self.b
r3 = R - self.c
# raft areas
s1 = R**2
s2 = (r3**2 + r2**2 + r3*r2)/3
s3 = (r2**2 + r1**2 + r2*r1)/3
s4 = self.re**2 - self.ri**2
# raft weight
mu = self.gb * h * np.pi
Pr = mu * (s1 + s2 + s3 + s4)
# soil areas
t1 = R**2 - r3**2
t2 = r3**2 - s2
t3 = R**2 - r2**2
t4 = r2**2 - s3
t5 = R**2 - self.re**2
# soil weight
Pt = self.gs * h * np.pi * (t1 + t2 + t3 + t4 + t5)
# Total vertical load
NT = self.P0 + Pr + Pt
return NT
def equi_inequality(self, R: TNum) -> TNum:
NT = self.vertical_load(R)
return NT*R - 4*self.M
def polynom(self) -> np.ndarray:
R_min = self.re + self.a + self.b + self.c
R = np.linspace(start=R_min, stop=1.1*R_min, num=4)
A = R[:, np.newaxis]**np.arange(3, -1, -1)
F = self.equi_inequality(R)
return np.linalg.solve(A, F)
def solve_R(self) -> typing.Tuple[
float, # lowest real root of R
np.ndarray, # polynomial coefficients
]:
coeffs = self.polynom()
roots = np.roots(coeffs)
mask = (np.abs(roots.imag) < 1e-8) & (roots.real > 0)
if not mask.any():
raise ValueError('No physical solution found')
return roots[mask].real.min(), coeffs
def main() -> typing.Dict[str, typing.Any]:
system = System(
H=4., a=1.05, b=2., c=1.,
ri=6.8/2, re=7.6/2, gb=2.5, gs=1.9,
P0=2492.52, M=9201.56,
)
R, coeffs = system.solve_R()
m, n, o, p = coeffs
out = {
"Polynomial": f"{m:.3f}·R3 {n:+.3f}·R2 {o:+.3f}·R {p:+.3f}",
"R (m)": R,
"NT": float(system.vertical_load(R)),
"Inequality": float(system.equi_inequality(R)),
"Coefficients": coeffs.tolist(),
}
pprint.pprint(out)
return out
if __name__ == '__main__':
OUT = main()
{'Coefficients': [26.389378290387,
3.0473448683118107,
2519.8895552440827,
-36806.24000012433],
'Inequality': -7.275957614183426e-12,
'NT': 4395.6427346302935,
'Polynomial': '26.389·R3 +3.047·R2 +2519.890·R -36806.240',
'R (m)': 8.373346566596176}
Comments
Explore related questions
See similar questions with these tags.
x^3is always dominant thencbrt(-a0/a3) ~ 23.5is a better starting guess. numpy has done its job the solved root is correct given the equation.rand heighthformulaV = pi*r^2*h/3for the concrete volumes and volume of a cylinder diameter Dr - volume of concrete for the soil volume for each of the layers in the diagram.P_r = density * pi * d^2 * h / 4with theds being diameters. Most water towers in my neighbourhood are round (most volume for least material cost).