RHEOLOGY

Overview

The RHEOLOGY function fits experimental data to established rheological models that describe the flow behavior of non-Newtonian fluids and glass-forming materials. Rheology is the study of how matter flows and deforms, particularly materials that exhibit complex viscosity behavior dependent on shear rate, stress, or temperature.

This implementation supports six fundamental models used across materials science, chemical engineering, and glass research:

Fluid Viscosity Models: - Bingham Viscoplastic: Describes yield-stress fluids that behave as rigid solids below a threshold stress, then flow linearly above it. Common in drilling muds, toothpaste, and mayonnaise. The model is: \tau = \tau_0 + \eta \dot{\gamma}. - Carreau Shear-Thinning: Models fluids with viscosity that decreases with increasing shear rate, transitioning between zero-shear and infinite-shear Newtonian plateaus. Used for polymer solutions and blood. - Cross Pseudoplastic: Similar to Carreau but with a different mathematical form for describing shear-thinning behavior in pseudoplastic fluids. - Herschel-Bulkley Yield Power Law: A generalization of the Bingham plastic that combines yield stress with power-law behavior: \tau = \tau_0 + K\dot{\gamma}^n. Widely used for food products, cosmetics, and drilling fluids.

Glass Viscosity Models: - MYEGA: The Mauro-Yue-Ellison-Gupta-Allan equation provides improved accuracy for glass viscosity-temperature relationships, particularly near the glass transition. - Vogel-Fulcher-Tammann (VFT): The classical equation for describing the super-Arrhenius temperature dependence of glass viscosity: \eta = 10^{A + B/(T - T_0)}. See VFT equation for details.

The function uses scipy.optimize.curve_fit with non-linear least squares optimization to determine optimal model parameters. It returns the fitted parameter values along with standard errors derived from the covariance matrix, enabling assessment of fit quality and parameter uncertainty.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=RHEOLOGY(xdata, ydata, rheology_model)
  • xdata (list[list], required): The xdata value
  • ydata (list[list], required): The ydata value
  • rheology_model (str, required): The rheology_model value

Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error string.

Examples

Example 1: Demo case 1

Inputs:

rheology_model xdata ydata
bingham_viscoplastic_yield_stress 0.1 0.8933205880269406
1.325 4.110814531343567
2.55 7.666783682520016
3.775 11.244045759090662
5 14.194223121451643

Excel formula:

=RHEOLOGY("bingham_viscoplastic_yield_stress", {0.1;1.325;2.55;3.775;5}, {0.8933205880269406;4.110814531343567;7.666783682520016;11.244045759090662;14.194223121451643})

Expected output:

y0 A
0.5994 2.754
0.1549 0.05024

Example 2: Demo case 2

Inputs:

rheology_model xdata ydata
carreau_shear_thinning_viscosity 0.1 2.75
1.3250000000000002 2.75
2.5500000000000003 2.75
3.7750000000000004 2.75
5 2.75

Excel formula:

=RHEOLOGY("carreau_shear_thinning_viscosity", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

A1 A2 t a n
2.75 2.75 1 1 0.5

Example 3: Demo case 3

Inputs:

rheology_model xdata ydata
cross_pseudoplastic_viscosity 0.1 2.75
1.3250000000000002 2.75
2.5500000000000003 2.75
3.7750000000000004 2.75
5 2.75

Excel formula:

=RHEOLOGY("cross_pseudoplastic_viscosity", {0.1;1.3250000000000002;2.5500000000000003;3.7750000000000004;5}, {2.75;2.75;2.75;2.75;2.75})

Expected output:

A1 A2 t m
2.75 2.75 1 1
0 0 0 0

Example 4: Demo case 4

Inputs:

rheology_model xdata ydata
herschel_bulkley_yield_power_law 0.1 163.22530195496407
1.325 0.01
2.55 783.2685882762112
3.775 4636.797356564116
5 17027.483843506292

Excel formula:

=RHEOLOGY("herschel_bulkley_yield_power_law", {0.1;1.325;2.55;3.775;5}, {163.22530195496407;0.01;783.2685882762112;4636.797356564116;17027.483843506292})

Expected output:

y0 K n
62.32 9.239 4.67
65.07 1.35 0.09016

Example 5: Demo case 5

Inputs:

rheology_model xdata ydata
myega_glass_viscosity_temperature 600 34.01159491
650 24.17713078
700 18.18279646
750 14.28350911
800 11.61245996

Excel formula:

=RHEOLOGY("myega_glass_viscosity_temperature", {600;650;700;750;800}, {34.01159491;24.17713078;18.18279646;14.28350911;11.61245996})

Expected output:

y0 K C
0.8 800 120
8.325e-9 0.000005533 0.000002495

Example 6: Demo case 6

Inputs:

rheology_model xdata ydata
vogel_fulcher_tammann_glass_transition 300 199.5262315
350 125.89254118
400 95.4992586
450 79.43282347
500 69.6397403

Excel formula:

=RHEOLOGY("vogel_fulcher_tammann_glass_transition", {300;350;400;450;500}, {199.5262315;125.89254118;95.4992586;79.43282347;69.6397403})

Expected output:

A B x0
1.5 120 150
9.422e-11 3.947e-8 3.225e-8

Python Code

import numpy as np
from scipy.optimize import curve_fit as scipy_curve_fit
import math

def rheology(xdata, ydata, rheology_model):
    """
    Fits rheology models to data using scipy.optimize.curve_fit. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html for details.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        xdata (list[list]): The xdata value
        ydata (list[list]): The ydata value
        rheology_model (str): The rheology_model value Valid options: Bingham Viscoplastic Yield Stress, Carreau Shear Thinning Viscosity, Cross Pseudoplastic Viscosity, Herschel Bulkley Yield Power Law, Myega Glass Viscosity Temperature, Vogel Fulcher Tammann Glass Transition.

    Returns:
        list[list]: 2D list [param_names, fitted_values, std_errors], or error string.
    """
    def _validate_data(xdata, ydata):
        """Validate and convert both xdata and ydata to numpy arrays."""
        for name, arg in [("xdata", xdata), ("ydata", ydata)]:
            if not isinstance(arg, list) or len(arg) < 2:
                raise ValueError(f"{name}: must be a 2D list with at least two rows")
            vals = []
            for i, row in enumerate(arg):
                if not isinstance(row, list) or len(row) == 0:
                    raise ValueError(f"{name} row {i}: must be a non-empty list")
                try:
                    vals.append(float(row[0]))
                except Exception:
                    raise ValueError(f"{name} row {i}: non-numeric value")
            if name == "xdata":
                x_arr = np.asarray(vals, dtype=np.float64)
            else:
                y_arr = np.asarray(vals, dtype=np.float64)

        if x_arr.shape[0] != y_arr.shape[0]:
            raise ValueError("xdata and ydata must have the same number of rows")
        return x_arr, y_arr

    # Model definitions dictionary
    models = {
        'bingham_viscoplastic_yield_stress': {
            'params': ['y0', 'A'],
            'model': lambda x, y0, A: y0 + A * x,
            'guess': lambda xa, ya: (float(np.min(ya)), max(float(np.ptp(ya) / (np.ptp(xa) if np.ptp(xa) else 1.0)), 1e-3)),
            'bounds': (0.0, np.inf),
        },
        'carreau_shear_thinning_viscosity': {
            'params': ['A1', 'A2', 't', 'a', 'n'],
            'model': lambda x, A1, A2, t, a, n: A2 + (A1 - A2) * np.power(1.0 + np.power(t * x, a), (n - 1.0) / a),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.min(ya)), 1.0, 1.0, 0.5),
            'bounds': (0.0, np.inf),
        },
        'cross_pseudoplastic_viscosity': {
            'params': ['A1', 'A2', 't', 'm'],
            'model': lambda x, A1, A2, t, m: A2 + (A1 - A2) / (1.0 + np.power(t * x, m)),
            'guess': lambda xa, ya: (float(np.max(ya)), float(np.min(ya)), 1.0, 1.0),
            'bounds': (0.0, np.inf),
        },
        'herschel_bulkley_yield_power_law': {
            'params': ['y0', 'K', 'n'],
            'model': lambda x, y0, K, n: y0 + K * np.power(np.clip(x, 0.0, None), n),
            'guess': lambda xa, ya: (float(np.min(ya)), 1.0, 1.0),
            'bounds': (0.0, np.inf),
        },
        'myega_glass_viscosity_temperature': {
            'params': ['y0', 'K', 'C'],
            'model': lambda x, y0, K, C: y0 * np.power(10.0, (K / x) * np.exp(C / x)),
            'guess': lambda xa, ya: (
                float(max(np.median(np.clip(ya, 1e-9, None)), 1e-6)),
                float(max(1e-3, 0.1 * np.ptp(np.log10(np.clip(ya, 1e-9, None))) * max(np.min(np.clip(xa, 1e-3, None)), 1e-3))),
                float(max(1e-3, 0.1 * max(np.min(np.clip(xa, 1e-3, None)), 1e-3)))
            ),
            'bounds': (0.0, np.inf),
        },
        'vogel_fulcher_tammann_glass_transition': {
            'params': ['A', 'B', 'x0'],
            'model': lambda x, A, B, x0: np.power(10.0, A + B / (x - x0 + 1e-9)),
            'guess': lambda xa, ya: (
                float(np.log10(max(np.mean(np.clip(ya, 1e-9, None)), 1e-9))),
                float(max(0.1, 0.1 * np.ptp(np.log10(np.clip(ya, 1e-9, None))) * max(np.median(np.clip(xa, 1e-6, None)), 1e-3))),
                float(max(min(np.clip(xa, 1e-6, None)) * 0.8, 1e-6))
            ),
            'bounds': ([-np.inf, 0.0, 0.0], np.inf),
        }
    }

    # Validate model parameter
    if rheology_model not in models:
        return f"Invalid model: {str(rheology_model)}. Valid models are: {', '.join(models.keys())}"

    model_info = models[rheology_model]

    # Validate and convert input data
    try:
        x_arr, y_arr = _validate_data(xdata, ydata)
    except ValueError as e:
        return f"Invalid input: {e}"

    # Perform curve fitting
    try:
        p0 = model_info['guess'](x_arr, y_arr)
        bounds = model_info.get('bounds', (-np.inf, np.inf))
        if bounds == (-np.inf, np.inf):
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, maxfev=10000)
        else:
            popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, bounds=bounds, maxfev=10000)

        fitted_vals = [float(v) for v in popt]
        for v in fitted_vals:
            if math.isnan(v) or math.isinf(v):
                return "Fitting produced invalid numeric values (NaN or inf)."
    except ValueError as e:
        return f"Initial guess error: {e}"
    except Exception as e:
        return f"curve_fit error: {e}"

    # Calculate standard errors
    std_errors = None
    try:
        if pcov is not None and np.isfinite(pcov).all():
            std_errors = [float(v) for v in np.sqrt(np.diag(pcov))]
    except Exception:
        pass

    return [model_info['params'], fitted_vals, std_errors] if std_errors else [model_info['params'], fitted_vals]

Online Calculator