scipy 用Carreau-Yasuda模型拟合流变数据

nhhxz33t  于 11个月前  发布在  其他
关注(0)|答案(1)|浏览(154)

我试图将一些流变学数据拟合到一个模型中。使用粘度和剪切速率的数组,我想将其拟合到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个参数。

3ks5zfa0

3ks5zfa01#

TL;DR

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

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))

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

字符串
然后我们有两个选择:使用LM算法,并对耦合常数(至少是幅度)提供一点帮助。

popt, pcov = optimize.curve_fit(
    model, 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(
    model, 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]


这两个方法返回的初始参数具有相当好的一致性。
100d 1xx 1c 1d 1x的字符串

相关问题