diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 19c9b45f7..1f739bc5b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -10,6 +10,7 @@ Change Log
 - [CHANGED] Colebrook-White friction model now uses scipy.optimize.newton for solving the equation
 - [CHANGED] Fluid properties compressibility and viscosity can now process pressure and temperature as inputs, if given the "allow_2d" attribute
 - [CHANGED] Temperature change along pipes now based on exponential function and does not use qext_w anymore
+- [ADDED] Churchill friction factor for application in entire flow regime
 
 [0.12.0] - 2025-06-27
 -------------------------------
diff --git a/src/pandapipes/converter/stanet/preparing_steps.py b/src/pandapipes/converter/stanet/preparing_steps.py
index 77a00f29f..d22868db8 100644
--- a/src/pandapipes/converter/stanet/preparing_steps.py
+++ b/src/pandapipes/converter/stanet/preparing_steps.py
@@ -171,7 +171,7 @@ def get_net_params(net, stored_data):
     :return: net parameters
     :rtype:
     """
-    known_friction_models = {1: "swamee-jain", 3: "nikuradse", 5: "colebrook"}
+    known_friction_models = {1: "swamee-jain", 3: "nikuradse", 5: "colebrook", 7:"churchill"}
     compressibility_models = {0: "linear", 1: "AGA", 2: "GERG-88"}
     net_params = dict()
     net_data = stored_data["net_parameters"]
diff --git a/src/pandapipes/pf/derivative_calculation.py b/src/pandapipes/pf/derivative_calculation.py
index b2034b275..81bb3523a 100644
--- a/src/pandapipes/pf/derivative_calculation.py
+++ b/src/pandapipes/pf/derivative_calculation.py
@@ -199,6 +199,13 @@ def calc_lambda(m, eta, d, k, gas_mode, friction_model, lengths, options, area):
         # 1.325 instead of 0.25???
         lambda_swamee_jain = 0.25 / ((np.log10(k / (3.7 * d) + 5.74 / (re ** 0.9))) ** 2)
         return lambda_swamee_jain, re
+    
+    elif friction_model == "churchill":
+        paramA = (-2.457*np.log((7/re)**0.9 + 0.27*k/d))**16
+        paramB = (37530/re)**16
+        lambda_churchill = 8*((8/re)**12 + 1/(paramA+paramB)**1.5)**(1/12)
+        return lambda_churchill, re
+
     else:
         # lambda_tot = np.where(re > 2300, lambda_laminar + lambda_nikuradse, lambda_laminar)
         lambda_tot = lambda_laminar + lambda_nikuradse
@@ -256,6 +263,23 @@ def calc_der_lambda(m, eta, d, k, friction_model, lambda_pipe, area, re, lengths
         lambda_der[pos] = 0.5 * np.log(10) ** 2 / (np.log(param) ** 3) / param * 5.166 * (
                     (eta[pos] * area[pos]) / (d[pos])) ** 0.9 * np.abs(m[pos]) ** -1.9
         return lambda_der
+    
+    elif friction_model == "churchill":
+        param = (7*eta[pos]*area[pos]/(m[pos]*d[pos]))**0.9 + 0.27*k[pos]/d[pos]
+        partial_dparamdm = -0.9*7**0.9 * (d[pos]/(eta[pos] * area[pos]))**(-0.9) * m[pos]**(-1.9)
+
+        paramsAB = ((-2.457*np.log((7*eta[pos]*area[pos]/(m[pos] * d[pos]))**0.9 + 0.27*k[pos]/d[pos]))**16 + (37530*eta[pos]*area[pos]/(m[pos] * d[pos]))**16)
+        paramC = (8*eta[pos]*area[pos]/(m[pos]*d[pos]))**12 + paramsAB**(-1.5)
+        
+        partial_dAdm = 16*(2.457*np.log(param))**15 * 2.457*param**(-1) * partial_dparamdm
+        partial_dBdm = -16*37530**16*(eta[pos]*area[pos]/(m[pos]*d[pos]))**17
+
+        partial_dCdm = -12*(8*eta[pos]*area[pos]/(m[pos]*d[pos]))**11 * (8*eta[pos]*area[pos]/(m[pos]**2*d[pos])) - paramsAB**(-3)*1.5*paramsAB**0.5 * (partial_dAdm + partial_dBdm)
+
+        lambda_der[pos] = 2/3* paramC**(-11/12) * partial_dCdm
+
+        return lambda_der
+
     else:
         lambda_der[pos] = -(64 * eta[pos] * area[pos]) / (m[pos] ** 2 * d[pos])
         return lambda_der
@@ -309,3 +333,8 @@ def cw_derivative(lambda_cb, re_nz, k_nz, d_nz):
         converged = np.all(res.converged)
 
     return converged, lambda_res
+
+
+
+
+