I performed some testing on your code. It appears that the solve function is having an issue with the sqrt function having as its argument the f symbol variable. I modified your code slightly and brought it to its minimum requirements for testing purposes. Entering a value that is independent of the f variable inside the sqrt function (inside the log function) appears to alleviate the previously observed issues. The temporary variable for now is temp.

I would advise you to recheck the formula to determine if it is set up correctly or if the f symbol variable is needed in that particular location of the formula. If it is needed there, then you may have to re-write the formula into smaller expressions, for example, to get the same expected results.

from sympy import symbols, solve, log, sqrt
# import sympy
from math import pi
f_T = 0.027
T = 60 # [F]
g = 32.2 # [ft/s^2]
H = 44 # [ft]
H_2 = 0 # [ft]
Q = 26 # [ft^3/s]
L = 1700 # [ft]
K_i = 0.5
K_l = 0.2
K_e = 1
D_i = 0.5 # [ft]
f_i = 0.02
print('D = ' + str(D_i) + ' [ft]')
epsilon = 0.0005
nu = 1.21*10**(-5) # [ft^2/s]
y = epsilon / D_i
for i in range(0, 5):
V = symbols('V')
D = symbols('D')
f = symbols('f')
A = (pi / 4) * D_i ** 2 # [ft^2] value of Area is obtained
V = solve(A * V - Q, V) # [ft/s] with Area velocity can be obtained
Re = V[0] * D_i / nu # with that velocity Reynolds is calculated
temp = float(input('Enter f: ')) # f is obtained from Moody Diagram
D = solve(H_2 + ((V[0] ** 2) / (2 * g)) * ((1700 / D) * f_i + (4 * K_l + K_i + K_e)) - H, D)
D_i = D[0] # [ft]
f = solve((-2 * log((epsilon / D_i) / 3.7 + 2.51 / (Re * sqrt(temp)), 10)) - 1 / sqrt(f), f)
print(f)

The equation has the form of those that can be solved using Lambert’s function. So, solve is computing a change of variable that turns the equation into the one defining Lambert’s function.

"""
Given an expression assumed to be in the form
``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
where X = g(x) and x = g^-1(X), return the Lambert solution,
``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.
"""

In your case (a,b,c,d,f)=(-2, 458717838856613, -8499189896411140, log(10), 42*log(10))

It is computing the part exp(c*d/a/b)*exp(-f/a) as exp(c*d-b*f)/a/b)

In the process, it falls into evaluating

13882669564194443*log(10)

and for this it evaluates

10**13882669564194443

The “never stops” is the time that it is evaluating this.

Try replacing some/all of the variables that are parameters of the equation, with symbols. That should make it more manageable. Unfortunately, even then there are going to be 9287 solutions.
I stopped it when it was computing all 9287 roots of a polynomial of the form x^{9287}-C=0, for certain constant C. It is taking a lot of time trying some heuristics. Maybe it will finish.