How to find the parameters of Herschel-Bulkley model through Nonlinear Regression!


(Maamar Dliouah) #1

Dear all,

First, I would like to thank everyone involved in this magnificent project, Math.NET saved my life!
I have few questions about the linear and nonlinear regression, I am a civil engineer and when I was working on my Master’s degree, I needed to develop a C# application that calculates the Rheological parameters of concrete based on data acquired from a test.
One of the models that describes the rheological behavior of concrete is the “Herschel-Bulkley model” and it has this formula :

y = T + K*x^n

x (the shear-rate), y (shear-stress) are the values obtained from the test, while T,K and N are the parameters I need to determine.
I know that the value of “T” is between 0 and Ymin (Ymin is the smallest data point from the test), so here is what I did:
Since it is nonlinear equation, I had to make it linear, like this : ln(y-T) = ln(K) + n*ln(x)
creat an array of possible values of T, from 0 to Ymin, and try each value in the equation, then through linear regression I find the values of K and N, then calculate the SSD, and store the results in an array, after I finish all the possible values of T, I see which one had the smallest SSD, and use it to find the K and N . this method works, but I feel it is not as smart or elegant as it should be, there must be a better way to do it, and I was hoping to find it here, it is also very slow.
here is the code that I used:

public static double HerschelBulkley(double shearRate, double tau0, double k, double n)
{
    var t = tau0 + k * Math.Pow(shearRate, n);
    return t;
}

public static (double Tau0, double K, double N, double DeltaMin, double RSquared) HerschelBulkleyModel(double[] shear, double[] shearRate, double step = 1000.0)
{
    // Calculate the number values from 0.0 to Shear.Min;
    var sm = (int) Math.Floor(shear.Min() * step);
    // Populate the Array of Tau0 with the values from 0 to sm
    var tau0Array = Enumerable.Range(0, sm).Select(t => t / step).ToArray();
    var kArray = new double[sm];
    var nArray = new double[sm];
    var deltaArray = new double[sm];
    var rSquaredArray = new double[sm];
    var shearRateLn = shearRate.Select(s => Math.Log(s)).ToArray();
    for (var i = 0; i < sm; i++)
    {
        var shearLn = shear.Select(s => Math.Log(s - tau0Array[i])).ToArray();
        var param = Fit.Line(shearRateLn, shearLn);
        kArray[i] = Math.Exp(param.Item1);
        nArray[i] = param.Item2;
        var shearHerschel = shearRate.Select(sr => HerschelBulkley(sr, tau0Array[i], kArray[i], nArray[i])).ToArray();
        deltaArray[i] = Distance.SSD(shearHerschel, shear);
        rSquaredArray[i] = GoodnessOfFit.RSquared(shearHerschel, shear);
    }
    var deltaMin = deltaArray.Min();
    var index = Array.IndexOf(deltaArray, deltaMin);
    var tau0 = tau0Array[index];
    var k = kArray[index];
    var n = nArray[index];
    var rSquared = rSquaredArray[index];
    return (tau0, k, n, deltaMin, rSquared);
}