将流变数据曲线拟合至 Carreau-Yasuda 模型

Question

我正在尝试将一些流变数据拟合到模型中。使用一系列粘度和剪切速率，我想将其拟合到 Carreau 模型，但是 scipy 中的 curve_fit 函数没有给我特别好的值（我已经绘制了之后收到的数据，它只遵循大约一半的数据的）。这是代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
from scipy.optimize import curve_fit
from scipy.odr import ODR, Model, Data, RealData
import seaborn as sns
sns.set(style="whitegrid")


df = pd.read_excel('Biomass.xls', sheet_name=1,)

data_array = df.to_numpy()
data_array_numbers = np.delete(data_array, 0, 0)
data_array_numbers = np.delete(data_array_numbers, 0, 0)
transposed_array = data_array_numbers.transpose()
sorted_array = transposed_array[:, np.argsort( transposed_array[1])]
shear_rate = np.array(sorted_array[1])
viscosity = np.array(sorted_array[2])


density=1000
for i in range(0, len(viscosity)):
    viscosity[i] = viscosity[i]/density
    
    
    
log_shear = np.zeros(shape=(len(shear_rate)))
log_viscosity = np.zeros(shape=(len(shear_rate)))
log_shear[i] = np.log(shear_rate[i])
log_viscosity[i] = np.log(viscosity[i])


def log_powerlaw(x, log_k, n_p):
    nu = log_k + n_p*x
    return nu

def powerlaw(x, k, n):
    return k*(x)**(n-1)

def carreau(x, n_inf, n_0, lamda, a, n):
    return n_inf + (n_0-n_inf)*(1+(x*lamda)**a)**((n-1)/a)


params, covariance = curve_fit(carreau, xdata=shear_rate, ydata=viscosity, maxfev=200000, method='trf')

n_inf_fit, n_0_fit, lambda_fit, a_fit, n_fit = params
print(n_inf_fit)
print(n_0_fit)
print(lambda_fit)
print(a_fit)
print(n_fit)
ax = plt.axes()
plt.scatter(shear_rate, viscosity, c='g', marker='x')
plt.plot(shear_rate, carreau(shear_rate, n_inf_fit, n_0_fit, lambda_fit, a_fit, n_fit ), c='b')
#plt.plot(shear_rate, powerlaw(shear_rate, k=params[0], n=params[1]) )
ax.set_title("Biomass Rheology")
ax.set_ylabel('Viscosity')
ax.set_xlabel('Shear Rate')
plt.yscale('log')
plt.show()

尝试想出一种线性化问题的方法，就像我对幂律拟合所做的那样（取对数使其成为直线方程），尝试给出不同的初始猜测。我期望得到一个很好的近似值，因为该模型有 5 个参数。

Answer 1

TL；博士

可以为这样的模型拟合参数，但它需要一些帮助，例如初始猜测或约束（至少对粘度常数进行粗略估计）。

MCVE

让我们生成一些数据：

import numpy as np
from scipy import optimize

def model(x, n_0, n_inf, lamda, a, n):
        return n_inf + (n_0 - n_inf) * np.power((1 + np.power((x[:, 0] * lamda), a)), ((n - 1) / a))

parameters = (4e-2, 5e-4, 1.6, 2.1, 0.3)
X = np.logspace(-4, 4, 30, base=10).reshape(-1, 1)
y = carreau(X, *parameters)
sigma = 1e-6 * np.ones(y.size)
y += sigma * np.random.normal(size=y.size)

然后我们有两个选择：使用 LM 算法，对耦合常数（至少幅度）有一些帮助。

popt, pcov = optimize.curve_fit(
    carreau, X, y, sigma=sigma, absolute_sigma=True,
    p0=(1e-2, 1e-4, 1, 1, 1)
)
# array([3.99998407e-02, 4.99859627e-04, 1.59989329e+00, 2.10012984e+00, 2.99977836e-01]

或者对问题添加约束。

popt, pcov = optimize.curve_fit(
    carreau, X, y, sigma=sigma, absolute_sigma=True,
    method="trf", bounds=[(1e-3, 1e-5, 0, 0, 0), (1e-1, 1e-3, 10, 3, 1)]
)
# array([3.99998407e-02, 4.99859627e-04, 1.59989329e+00, 2.10012983e+00, 2.99977836e-01]

两种方法返回的初始参数具有相当好的一致性。