CN112749488B - 一种纤维素生物墨水黏度特性预测方法 - Google Patents

Abstract

本发明公开了一种纤维素生物墨水黏度特性预测方法，步骤包括S1、使用旋转流变仪对纤维素生物墨水进行剪切扫描试验获得剪切速率‑黏度数据；S2、根据步骤S1获得的剪切速率‑黏度曲线得知纤维素生物墨水的黏度特性提出候选黏度模型；S3、通过定义优化变量，确定优化问题的约束以及优化目标，将候选黏度模型参数确定问题转化为多目标优化问题，以确定各候选黏度模型的参数；S4、利用改进的非支配排序遗传算法求解转化的多目标优化问题得到Pareto最优解；S5、利用逼近理想解排序法从Pareto最优解选择单一最优解从而确定最适合的黏度模型及黏度模型参数；该预测方法相对于现有技术能够精确且可靠的对纤维素生物墨水黏度特性进行准确预测。

Description

一种纤维素生物墨水黏度特性预测方法

技术领域

本发明涉及增材制造用生物墨水技术领域，特别涉及一种纤维素生物墨水黏度特性预测方法。

背景技术

生物墨水被定义为包含生物组分或生物材料的混合物凝胶，适合采用增材制造技术进行3D打印。近年来，由植物组分(纤维素，木质素等)制备的生物墨水得益于价低、无毒无害、可持续性以及来源充足等优点受到了学术界、生物界和工业界的广泛关注被成功应用于组织培养、康复工程等医学领域以及电池、传感器、记忆材料等工业领域。纤维素生物墨水主要由来源于棉花的纤维素和纤维素衍生物制备，具备良好的打印性能和机械性能，是由植物组分制备的生物墨水的典型代表，也是目前使用最广泛的增材制造用生物墨水。

纤维素生物墨水的黏度特性影响其在增材制造过程中的成型并最终影响产品的精度。因此，纤维素生物墨水的黏度特性预测是其制备、打印性能评估、打印过程参数优化以及打印产品性能评估等工作的第一个关键步骤。纤维素生物墨水是一种剪切变稀流体，黏度随着剪切速率增大而减小。黏度特性表征了材料剪切速率和黏度的定量关系，黏度特性预测包括黏度模型选择和模型参数确定两个部分。根据公开的期刊文献资料(J.L.Dávila and M.A.

"Rheological evaluation of Laponite/alginate inks for 3Dextrusion-based printing,"The International Journal of Advanced ManufacturingTechnology,vol.101,no.1,pp.675-686,2019/03/01 2019.)，现有的纤维素生物墨水黏度特性预测方法可总结为：首先，对制备的纤维素生物墨水进行剪切速率扫描试验获得不同剪切速率下的黏度数据；然后，选择幂律模型作为纤维素生物墨水的黏度模型；最后，利用两边取对数将幂律模型转化为剪切速率和黏度的线性方程并使用线性拟合方法确定模型参数。

现有的纤维素生物墨水黏度特性预测方法存在以下两个问题：

(1)缺少黏度模型的选择：纤维素生物墨水是剪切变稀材料，目前有四种广泛使用的经典模型(幂律模型，Hershel-Bulkley模型，Bird-Carreau模型和Cross幂律模型)被应用于剪切变稀材料的黏度模型，现有方法仅简单选用幂律模型，缺少模型的对比和选用，最终造成预测结果不准确；

(2)模型参数确定中线性拟合方法的限制：现有方法使用线性拟合方法确定模型参数，但Bird-Carreau模型和Cross幂律模型不能转换为线性方程，线性拟合方法无法通用。

发明内容

本发明的目的是提供一种解决上述技术问题，实现对纤维素生物墨水黏度特性进行准确且可靠预测的纤维素生物墨水黏度特性预测方法。

为此，本发明技术方案如下：

一种纤维素生物墨水黏度特性预测方法，步骤如下：

S1、使用旋转流变仪对纤维素生物墨水进行剪切扫描试验获得剪切速率-黏度数据；

S2、根据步骤S1获得的剪切速率-黏度曲线得知纤维素生物墨水的黏度特性提出候选黏度模型；

S3、通过定义优化变量，确定优化问题的约束以及优化目标，将候选黏度模型参数确定问题转化为多目标优化问题，以确定各候选黏度模型的参数；

S4、利用改进的非支配排序遗传算法求解转化的多目标优化问题得到Pareto最优解，以在候选黏度模型中确定最适合纤维素生物墨水黏度特性的黏度模型；

S5、利用逼近理想解排序法从Pareto最优解选择单一最优解从而确定最适合的黏度模型及黏度模型参数。

进一步地，步骤S2的具体实施步骤为：

S201、根据步骤S1获得的剪切速率-黏度曲线得知纤维素生物墨水的黏度特性为剪切变稀特性；

S202、选取四种剪切变稀模型作为纤维素生物墨水的候选黏度模型，分别为：

(1)幂率模型，即PLE模型，其表达式为：

式中，μ为动力黏度，单位：Pa·s；

为剪切速率，单位：s^-1；k_P为PLE模型中的稠度系数，单位：Pa·s；n_P为PLE模型中的幂率指数，单位：无量纲；

(2)Hershel-Bulkley模型，即HBE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,H为HBE模型中的极限黏度，单位：Pa·s；τ_0,H为HBE模型中的屈服应力，单位：Pa；k_H为HBE模型中的稠度系数，单位：Pa·s；n_H为HBE模型中的流动系数，单位：无量纲；

(3)Bird-Carreau模型，即BCE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,B为BCE模型中的零切黏度，单位：Pa·s；μ_∞,B为BCE模型中的无限剪切黏度，单位：Pa·s；k_B为BCE模型中的驰豫时间，单位：s；n_B为BCE模型中的幂率指数，单位：无量纲；

(4)Cross幂律模型，即CPLE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,C为CPLE中的零切黏度，单位：Pa·s；μ_∞,C为CPLE模型中的无限剪切黏度，单位：Pa·s；k_C为CPLE模型中的驰豫时间，单位：s；n_C为CPLE模型中的幂率指数，单位：无量纲。

进一步地，步骤S3的具体实施步骤为：

S301、以每个候选黏度模型中的参数为优化变量：

(1)在PLE模型中，优化变量为PLE模型中的稠度系数k_P和幂率指数n_P；

(2)在HBE模型中，优化变量为HBE模型中的极限黏度μ_0,H，屈服应力τ_0,H，稠度系数k_H和流动系数n_H；

(3)在BCE模型中，优化变量为BCE模型中的零切黏度μ_0,B，无限剪切黏度μ_∞,B，驰豫时间k_B和幂率指数n_B；

(4)在CPLE模型中，优化变量为CPLE中的零切黏度μ_0,C，无限剪切黏度μ_∞,C，驰豫时间k_C和幂率指数n_C；

S302、确定优化问题的约束：

四个候选黏度模型的约束条件具体为：

(1)PLE模型的约束条件为：

(2)HBE模型的约束条件为：

(3)BCE模型的约束条件为：

(4)CPLE模型的约束条件为：

S303、确定优化目标：以决定系数R²的最大化以及绝对平均相对偏差AARD％的最小化为目标，即：

优化目标为：min[-R²,AARD％]；

其中，-R²和AARD％为两个目标函数，二者分别定义为：

式中，μ_i,e为黏度试验值，μ_i,p为黏度预测值，N为试验数据点数，

为黏度平均试验值。

进一步地，步骤S4的具体实施步骤为：

S401、确定种群大小N_p，并随机产生初始种群P(t)|t＝0作为父代种群；其中，采用双精度向量作为种群类型并规定种群大小N_p为200；

S402、使用遗传算法算子对父代种群P(t)依次进行选择操作，交叉操作和变异操作，以产生子代种群Q(t)；其中，选择操作采用锦标赛选择算法；交叉操作中交叉概率设为0.8；变异操作中变异概率设为0.02；

S403、将父代种群P(t)与子代种群Q(t)合并生成新的集合R(t)；

S404、对R(t)依次进行快速非支配排序和拥挤度计算以选出新的父代种群P(t+1)；

S405、判断是否达到终止条件：若达到，则最新的父代种群P(t+1)为Pareto最优解；若未达到，则重复S402～S404步骤，直至达到终止条件；其中，终止条件设置为最大代数为20000。

进一步地，步骤S5的具体实施步骤为：

S501、构造m行n列的矩阵(x_ij)_m×n；其中，m为步骤S4中获得的Pareto最优解中解的个数(即达到终止条件时最新的父代种群P(t+1)的大小)，n为目标函数个数，n＝2，x_ij为第i个解的第j个目标函数的目标函数值；

S502、将矩阵(x_ij)_m×n通过下式处理生成正则化矩阵(t_ij)_m×n：

S503、通过下式处理正则化矩阵(t_ij)_m×n得到权重正则化矩阵(a_ij)_m×n：

a_ij＝w_j×t_ij,i＝1,2,...,m,j＝1,2,...,n，

式中，w_j为第j个目标函数的相对权重，由下式获得：

S504、通过下两式确定正理想解A⁺和负理想解A^-：

A⁺＝[min(a₁₁,..,a_m1),min(a₁₂,..,a_m2),..,min(a_1n,..,a_mn)],

A^-＝[max(a₁₁,..,a_m1),max(a₁₂,..,a_m2),..,max(a_1n,..,a_mn)]；

S505、分别按下两式计算Pareto最优解中每个解距离正理想解A⁺和负理想解A-的距离：

其中，

为Pareto最优解中第i个解距离正理想解A⁺的距离；

为多目标优化问题的正理想解；

为Pareto最优解中第i个解距离负理想解A^-的距离；

为多目标优化问题的负理想解；

S506、计算Pareto最优解中每个解距离正、负理想解的相对距离C_i：

S507、将Pareto最优解中的所有解距离正、负理想解的相对距离C_i进行排序，取让C_i最大的解为单一最优解A_final：A_final＝A∈max(C_i)，该单一最优解即为最终确定的模型参数。

与现有技术相比，该纤维素生物墨水黏度特性预测方法针对目前现有的纤维素生物墨水黏度特性预测技术中存在的缺少黏度模型的选择，模型参数确定中线性拟合方法的限制这两个问题，提出了新的纤维素生物墨水黏度特性预测方法以获得准确的纤维素生物墨水的黏度模型及模型参数，最后通过以本发明得到的纤维素生物墨水黏度特性为输入条件的仿真结果与试验结果对比验证该方法相对于现有技术能够精确且可靠的对纤维素生物墨水黏度特性进行准确预测。

附图说明

图1为本发明的纤维素生物墨水黏度特性预测方法的流程图；

图2为本发明实施例中的某纤维素生物墨水在黏度特性预测中通过S1步骤获得的数据绘制的剪切速率-黏度曲线图；

图3为本发明的纤维素生物墨水黏度特性预测方法中的S4步骤中的改进的非支配排序遗传算法方法的流程图；

图4为本发明的实施例中的某纤维素生物墨水在黏度特性预测中通过S4步骤获得的Pareto最优解示意图；

图5(a)为本发明的实施例中的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对PLE模型的单一最优解示意图；

图5(b)为本发明的实施例中的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对HBE模型的单一最优解示意图；

图5(c)为本发明的实施例中的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对BCE模型的单一最优解示意图；

图5(d)为本发明的实施例中的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对CPLE模型的单一最优解示意图；

图6为本发明的实施例中黏度特性预测的验证中打印纤维素生物墨水用的挤出式3D打印机示意图；

图7为本发明的实施例中黏度特性预测的验证中所使用的挤出式3D打印机中的挤出装置示意图；

图8(a)为本发明的实施例中黏度特性预测的验证方法中的数值仿真中的挤出过程三维模型示意图；

图8(b)为本发明的实施例中黏度特性预测的验证方法中的数值仿真步骤中的挤出过程简化二维轴对称模型示意图；

图8(c)为本发明的实施例中黏度特性预测的验证方法中的数值仿真步骤中的挤出过程网格划分后的简化二维轴对称模型示意图；

图9为本发明的实施例在黏度特性预测后进行验证获得的挤出的在空气中的某纤维素生物墨水在试验与仿真中的外形图；

图10为本发明的实施例在黏度特性预测后进行验证获得的某纤维素生物墨水在试验与仿真中的喷嘴底面速度曲线图；

图11为本发明的实施例在黏度特性预测后进行验证获得的某纤维素生物墨水试验与仿真速度偏差曲线图。

其中，1、三轴运动平台，2、活塞，3、挤出装置，301、针管，302、喷嘴，4、基板，5、纤维素生物墨水。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

如图1所示，对某纤维素生物墨水进行黏度特性预测为例，对本申请的纤维素生物墨水黏度特性预测方法的具体实施步骤进行描述，步骤如下：

S1、使用旋转流变仪对纤维素生物墨水进行剪切扫描试验获得剪切速率-黏度数据；具体地，

在本步骤中，旋转流变仪采用旋转流变仪ARES(德州仪器，美国)，并选用直径为25mm的两块圆形平板结构作为黏度特性测试的几何结构，设置平行设置的两块圆形平板之间的间距设置为1mm，设置测试温度为25℃，设置流变仪为剪切扫描模式，设置剪切速率从0.0125s^-1到100s^-1，按对数间隔进行20组试验；为减少试验误差，共进行三次测试，取三次测试结果的平均值为剪切速率-黏度数据；

如图2所示为该某纤维素生物墨水经测试得到的剪切速率-黏度曲线；

S2、根据步骤S1得到的剪切速率-黏度曲线确定纤维素生物墨水的剪切变稀黏度特性提出候选黏度模型；

具体地，步骤S2的具体步骤如下：

S201、根据如图2所示的某纤维素生物墨水剪切速率-黏度曲线可知，当剪切速率很小时，黏度保持在903Pa·s的常数；而当剪切速率高于0.2s^-1时，黏度随着剪切速率的增大而减小，表现为剪切变稀特性；

S202、选取四种常用的剪切变稀模型作为纤维素生物墨水的候选黏度模型，分别为：(1)幂率模型(以下简称为PLE模型)；(2)Hershel-Bulkley模型(以下简称为HBE模型)；(3)Bird-Carreau模型(以下简称为BCE模型)；(4)Cross幂律模型(以下简称为CPLE模型)；其中，

(1)PLE模型的表达式为：

式中，μ为动力黏度，单位：Pa·s；

为剪切速率，单位：s^-1；k_P为PLE模型中的稠度系数，单位：Pa·s；n_P为PLE模型中的幂率指数，单位：无量纲；

(2)HBE模型的表达为：

式中，μ为动力黏度，单位：s^-1；μ_0,H为HBE模型中的极限黏度，单位：Pa·s；τ_0,H为HBE模型中的屈服应力，单位：Pa；k_H为HBE模型中的稠度系数，单位：Pa·s；n_H为HBE模型中的流动系数，单位：无量纲；

(3)BCE模型的表达为：

式中，μ为动力黏度，单位：s^-1；μ_0,B为BCE模型中的零切黏度，单位：Pa·s；μ_∞,B为BCE模型中的无限剪切黏度，单位：Pa·s；k_B为BCE模型中的驰豫时间，单位：s；n_B为BCE模型中的幂率指数(单位：无量纲)；

(4)CPLE模型的表达为：

式中，μ为动力黏度，单位：s^-1；μ_0,C为CPLE中的零切黏度，单位：Pa·s；μ_∞,C为CPLE模型中的无限剪切黏度，单位：Pa·s；k_C为CPLE模型中的驰豫时间，单位：s；n_C为CPLE模型中的幂率指数，单位：无量纲；

S3、为确定候选黏度模型参数，进一步通过定义优化变量，确定优化问题的约束以及优化目标，将候选黏度模型参数确定问题转化为多目标优化问题；

具体地，步骤S3的具体步骤为：

S301、以每个候选黏度模型中的参数为优化变量，具体地：

(1)在PLE模型中，优化变量为PLE模型中的稠度系数k_P和PLE中的幂率指数n_P；

(2)在HBE模型中，优化变量为HBE模型中的极限黏度μ_0,H，HBE模型中的屈服应力τ_0,H，HBE模型中的稠度系数k_H和HBE模型中的流动系数n_H；

(3)在BCE模型中，优化变量为BCE模型中的零切黏度μ_0,B，BCE模型中的无限剪切黏度μ_∞,B，BCE模型中的驰豫时间k_B和BCE模型中的幂率指数n_B；

(4)在CPLE模型中，优化变量为CPLE中的零切黏度μ_0,C，CPLE模型中的无限剪切黏度μ_∞,C，CPLE模型中的驰豫时间k_C和CPLE模型中的幂率指数n_C；

S302、确定优化问题的约束：

四个候选黏度模型的约束条件具体为：

(1)PLE模型的约束条件为：

(2)HBE模型的约束条件为：

(3)BCE模型的约束条件为：

(4)CPLE模型的约束条件为：

S303、确定优化目标：以决定系数R²的最大化以及绝对平均相对偏差(AARD％)的最小化为目标，即：

优化目标为：min[-R²,AARD％]；

其中，-R²和AARD％为两个目标函数，二者分别定义为：

上两式中，μ_i,e为黏度试验值，μ_i,p为黏度预测值，N为试验数据点数，

为黏度平均试验值；

S4、利用改进的非支配排序遗传算法求解转化的多目标优化问题得到Pareto最优解。

如图3所示，步骤S4的具体步骤为：

S401、确定改进的非支配排序遗传算法中种群大小N_p，并随机产生初始种群P(t)|t＝0作为父代种群；其中，采用双精度向量作为种群类型并规定种群大小N_p为200；

S402、使用遗传算法算子对父代种群P(t)依次进行选择操作，交叉操作和变异操作，以产生子代种群Q(t)；其中，选择操作采用锦标赛选择算法；交叉操作中交叉概率设为0.8；变异操作中变异概率设为0.02；

S403、将父代种群P(t)与子代种群Q(t)合并生成新的集合R(t)；

S404、对R(t)依次进行快速非支配排序和拥挤度计算以选出新的父代种群P(t+1)；

S405、判断是否达到终止条件：

若达到，则最新的父代种群P(t+1)为Pareto最优解；

若未达到，则重复S402～S404步骤，直至达到终止条件；

其中，终止条件设置为最大代数为20000。

如图4所示为按照S4步骤获得的本实施例的某纤维素生物墨水的四种候选黏度模型的Pareto最优解。

具体地，PLE模型的AARD％从19.00％到22.41％，R²从0.9100到0.9223；HBE模型的AARD％从4.79％到14.60％，R²从0.9832到0.9971；BCE模型的AARD％从3.51％到4.85％，R²从0.9987到0.9993；CPLE模型的AARD％从5.85％到14.70％，R²从0.9873到0.9954；从上述结果可知，BCE模型有最小的AARD％和最大的R²，因此，BCE模型初选为本实施例中的某纤维素生物墨水的黏度模型；

S5、利用逼近理想解排序法从Pareto最优解选择单一最优解从而确定最适合的黏度模型及黏度模型参数；

具体地，步骤S5的具体步骤为：

S501、构造m行n列的矩阵(x_ij)_m×n；其中，m为步骤S4中获得的Pareto最优解中解的个数(即达到终止条件时最新的父代种群P(t+1)的大小)，n为目标函数个数，x_ij为第i个解的第j个目标函数的目标函数值；其中，由于有两个目标函数，因此n的值为2；

S502、将矩阵(x_ij)_m×n通过下式处理生成正则化矩阵(t_ij)_m×n：

S503、通过下式处理正则化矩阵(t_ij)_m×n得到权重正则化矩阵(a_ij)_m×n：

a_ij＝w_j×t_ij,i＝1,2,...,m,j＝1,2,...,n，

式中，w_j为第j个目标函数的相对权重，由下式获得：

S504、通过下两式确定正理想解A⁺和负理想解A^-：

A⁺＝[min(a₁₁,..,a_m1),min(a₁₂,..,a_m2),..,min(a_1n,..,a_mn)],

A^-＝[max(a₁₁,..,a_m1),max(a₁₂,..,a_m2),..,max(a_1n,..,a_mn)]；

S505、分别按下两式计算Pareto最优解中每个解距离正理想解A⁺和负理想解A^-的距离：

其中，

为Pareto最优解中第i个解距离正理想解A⁺的距离；

为多目标优化问题的正理想解；

为Pareto最优解中第i个解距离负理想解A^-的距离；

为多目标优化问题的负理想解；

S506、计算Pareto最优解中每个解距离正、负理想解的相对距离C_i：

S507、将Pareto最优解中的所有解距离正、负理想解的相对距离C_i进行排序，取让C_i最大的解为单一最优解A_final：A_final＝A∈max(C_i)，该单一最优解即为最终确定的模型参数；最后，通过对比每个候选黏度模型及其参数的目标函数，选择目标函数最优的模型及其参数为纤维素生物墨水的黏度特性预测结果，完成纤维素生物墨水的黏度特性预测。

如图5(a)所示为该实施例的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对PLE模型的单一最优解示意图。如图5(b)所示为该实施例的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对HBE模型的单一最优解示意图。如图5(c)所示为该实施例的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对BCE模型的单一最优解示意图。如图5(d)所示为该实施例的某纤维素生物墨水在黏度特性预测中通过S5步骤获得的针对CPLE模型的单一最优解示意图。在上述四幅附图中，实心圆点代表正理想解，正理想解为一个构造解，其拥有在Pareto解中最小的目标函数值；空心圆点代表负理想解，负理想解为一个构造解，其拥有在Pareto解中最大的目标函数值；实心方点为经过S5步骤得到的单一最优解，其特点是到正理想解的距离最近且到负理想解的距离最远。

因此，根据步骤S5和上述四张附图可知，该实施例的某纤维素生物墨水的黏度特性预测中，各个候选黏度模型的模型参数和目标函数值为：

(1)PLE模型：模型参数为：k_P＝73.39Pa·s，n_P＝0.38，目标函数值为：AARD％＝19.95％，R²＝0.9171；

(2)HBE模型：模型参数为：μ_0,H＝865.57Pa·s，τ_0,H＝0.049Pa，k_H＝91.05Pa·s，n_H＝0.32，目标函数值为：AARD％＝7.80％，R²＝0.9945；

(3)BCE模型：模型参数为：μ_0,B＝956.81Pa·s，μ_∞,B＝0.0022Pa·s，k_B＝26.63s，n_B＝0.29，目标函数值为：AARD％＝3.87％，R²＝0.9993；

(4)CPLE模型：模型参数为：μ_0,C＝1105.81Pa·s，μ_∞,C＝2.51Pa·s，k_C＝13.94s，n_C＝0.91，目标函数值为：AARD％＝9.09％，R²＝0.9958。

基于上述结果可以确定该实施例的某纤维素生物墨水的黏度特性预测中，BCE模型及其确定的模型参数对应的AARD％最小，而R²最大；因此，该某纤维素生物墨水的黏度特性预测结果为：模型为BCE模型，模型参数为μ_0,B＝956.81Pa·s，μ_∞,B＝0.0022Pa·s，k_B＝26.63s，n_B＝0.29。

为验证本发明提供的纤维素生物墨水黏度特性预测方法的准确性和可靠性，采用以本发明方法得到的黏度特性预测为输入条件的计算机仿真结果与实际试验结果对比的方法对本申请的纤维素生物墨水黏度特性预测结果的准确性进行验证和判断。其具体验证方法的步骤如下：

(1)将该实施例的某纤维素生物墨水放入挤出式3D打印机内进行打印，用相机获得挤出的某纤维素生物墨水的形状，并通过质量守恒公式获得理论的挤出的某纤维素生物墨水在喷嘴底面的速度作为其在喷嘴底面的速度参考值；

挤出式3D打印机采用法国Tobeca公司和法国洛林大学联合制造的TM-081型挤出式3D打印机。如图6所示，挤出式3D打印机包括三轴运动平台1、活塞2、挤出装置3、基板4和纤维素生物墨水5；如图7所示，挤出装置3包括针管301和喷嘴302。

如图6和图7所示，该挤出式3D打印机进行打印过程为：将某纤维素生物墨水5放入挤出式3D打印机内部的挤出装置3的针管301内，将针管301和喷嘴302设置为固定模式(即将三轴运动平台1固定不动)，喷嘴302底部到基板4的距离设置为20mm，活塞2速度设置为1.06×10^-2mm/s；在挤出过程中，用相机记录挤出的某纤维素生物墨水5在空气中的挤出形状，并利用如下的质量守恒定律计算出喷嘴302底面某纤维素生物墨水5的理论速度值：

式中，d_s和d_n分别为针管301的内径和喷嘴302的内径，其值分别为21.6mm和0.84mm；ν_p和ν_f分别为活塞2的速度和喷嘴302底面某纤维素生物墨水5的理论速度值；由ν_p设置为1.06×10^-2mm/s计算得到ν_f为7mm/s。

(2)对步骤(1)中的打印过程，以本发明方法得到的黏度特性预测为输入条件，用开源仿真软件OpenFOAM进行数值仿真并利用开源后处理软件ParaView进行仿真结果后处理，获得挤出的某纤维素生物墨水5在空气中的挤出形状和其在喷嘴302底面的速度的仿真结果；

具体地，首先，参见图8(a)，根据针管301和喷嘴302的尺寸构建挤出过程的三维模型；然后，考虑到在某纤维素生物墨水5从针管301挤出到空气中的3D打印过程为二维轴对称流动，因此通过取三维模型周向的1/60，将图8(a)所示的三维模型简化为图8(b)所示的二维轴对称模型；接着，将图8(b)所示的二维轴对称模型采用楔形网格进行网格划分成如图8(c)所示的划分网格后的二维轴对称模型，并导入OpenFOAM软件中；之后，在OpenFOAM软件中，以黏度特性预测为输入条件利用OpenFOAM中的InterFOAM工具进行数值仿真计算计算结果，并利用ParaView软件对计算结果进行后处理，得到某纤维素生物墨水5在3D打印过程中挤出的形状和其在喷嘴302底面速度的仿真值。

(3)对步骤(1)得到的试验结果和步骤(2)得到的以本发明方法得到的黏度特性预测为输入条件的数值仿真结果进行对比验证黏度特性预测的正确性。具体地，参加附图9，以本发明方法得到的黏度特性预测为输入条件获得的纤维素生物墨水在3D打印过程中挤出的形状仿真结果与试验结果有很好的一致性；参加附图10，以本发明方法得到的黏度特性预测为输入条件获得的纤维素生物墨水在3D打印过程中喷嘴底面的速度的仿真结果为7.2mm/s，试验理论结果为7mm/s，参加附图11，以本发明方法得到的黏度特性预测为输入条件获得的纤维素生物墨水在3D打印过程中喷嘴底面的速度的仿真与试验理论结果的差值仅为0.2mm/s。该仿真结果与试验结果相比，在纤维素生物墨水在3D打印过程中挤出的形状和在喷嘴底面的速度上，以本发明方法得到的黏度特性预测为输入条件的仿真结果与试验结果均有很好的一致性。可见，该纤维素生物墨水黏度特性预测方法预测准确率高。

进一步，为了验证本发明的一种纤维素生物墨水黏度特性预测方法的准确性，采用现有方法对该某纤维素生物墨水黏度特性进行预测，以对本申请对纤维素生物墨水黏度特性预测结果的准确率进行比较。

具体地，现有方法获得某纤维素生物墨水黏度特性具体步骤为：首先，选用PLE模型作为黏度模型：即：

然后，利用两边对数化将

转化为线性方程：

最后，通过线性拟合方法获得参数如下：k_P＝83.97Pa·s，n_P＝0.33。

在本实施例中，同时将现有方法获得的黏度特性预测作为输入条件获得仿真结果，并与将本发明的黏度特性预测方法获得的黏度特性预测作为输入条件获得仿真结果一起与试验结果对比，以此验证本发明黏度特性预测方法的准确性。

下面具体描述本实施例中某纤维素生物墨水黏度特性预测的验证方法中步骤(3)的结果：首先，进行某纤维素生物墨水5在3D打印中挤出的形状的对比，参见图9，利用本发明的纤维素生物墨水黏度特性预测方法获得的黏度特性得到的挤出形状仿真结果明显比现有方法得到的黏度特性得到的挤出形状仿真结果更接近试验结果；然后，进行某纤维素生物墨水5在3D打印中在喷嘴302底部的速度对比，参见图10，利用本发明的纤维素生物墨水黏度特性预测方法获得的黏度特性得到的在喷嘴302底部的速度明显比现有方法得到的黏度特性得到的在喷嘴302底部的速度更接近理论参考值；为更好对比在喷嘴302底部的速度，将某纤维素墨水5在喷嘴302底面的速度仿真结果减去理论参考值得到速度偏差，参见图11，利用本发明的纤维素生物墨水黏度特性预测方法获得的黏度特性得到的在喷嘴302底部的速度与理论参考值的偏差仅为0.2mm/s，远优于现有方法得到的黏度特性得到的在喷嘴302底部的速度与理论参考值的0.7mm/s的偏差。

综上所述，本发明的一种纤维素生物墨水黏度特性预测方法具体较好的预测精度，并且通过本发明的纤维素生物墨水黏度特性预测的验证方法验证了本发明的一种纤维素生物墨水黏度特性预测方法的精确性，确保了黏度特性预测的可靠性。

本发明未详细公开的部分属于本领域的公知技术。另外，尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均为保护之列。

Claims (2)

1.一种纤维素生物墨水黏度特性预测方法，其特征在于，步骤如下：

S1、使用旋转流变仪对纤维素生物墨水进行剪切扫描试验获得剪切速率-黏度数据；

S2、根据步骤S1获得的剪切速率-黏度曲线得知纤维素生物墨水的黏度特性提出候选黏度模型；步骤S2的具体实施步骤为：

S201、根据步骤S1获得的剪切速率-黏度曲线得知纤维素生物墨水的黏度特性为剪切变稀特性；

S202、选取四种剪切变稀模型作为纤维素生物墨水的候选黏度模型，分别为：

(1)PLE模型，其表达式为：

式中，μ为动力黏度，单位：Pa·s；

为剪切速率，单位：s^-1；k_P为PLE模型中的稠度系数，单位：Pa·s；n_P为PLE模型中的幂率指数，单位：无量纲；

(2)HBE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,H为HBE模型中的极限黏度，单位：Pa·s；τ_0,H为HBE模型中的屈服应力，单位：Pa；k_H为HBE模型中的稠度系数，单位：Pa·s；n_H为HBE模型中的流动系数，单位：无量纲；

(3)BCE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,B为BCE模型中的零切黏度，单位：Pa·s；μ_∞,B为BCE模型中的无限剪切黏度，单位：Pa·s；k_B为BCE模型中的驰豫时间，单位：s；n_B为BCE模型中的幂率指数，单位：无量纲；

(4)CPLE模型，其表达为：

式中，μ为动力黏度，单位：Pa·s；μ_0,C为CPLE中的零切黏度，单位：Pa·s；μ_∞,C为CPLE模型中的无限剪切黏度，单位：Pa·s；k_C为CPLE模型中的驰豫时间，单位：s；n_C为CPLE模型中的幂率指数，单位：无量纲；

S3、通过定义优化变量，确定优化问题的约束以及优化目标，将候选黏度模型参数确定问题转化为多目标优化问题，以确定各候选黏度模型的参数；步骤S3的具体实施步骤为：

S301、以每个候选黏度模型中的参数为优化变量：

(1)在PLE模型中，优化变量为PLE模型中的稠度系数k_P和幂率指数n_P；

(2)在HBE模型中，优化变量为HBE模型中的极限黏度μ_0,H，屈服应力τ_0,H，稠度系数k_H和流动系数n_H；

(3)在BCE模型中，优化变量为BCE模型中的零切黏度μ_0,B，无限剪切黏度μ_∞,B，驰豫时间k_B和幂率指数n_B；

(4)在CPLE模型中，优化变量为CPLE中的零切黏度μ_0,C，无限剪切黏度μ_∞,C，驰豫时间k_C和幂率指数n_C；

S302、确定优化问题的约束：PLE模型的约束条件为：

(2)HBE模型的约束条件为：

BCE模型的约束条件为：

CPLE模型的约束条件为：

S303、确定优化目标：以决定系数R²的最大化以及绝对平均相对偏差AARD％的最小化为优化目标：min[-R²,AARD％]；其中，-R²和AARD％为两个目标函数，二者分别定义为：

式中，μ_i,e为黏度试验值，μ_i,p为黏度预测值，N为试验数据点数，

为黏度平均试验值；

S4、利用改进的非支配排序遗传算法求解转化的多目标优化问题得到Pareto最优解；

S5、利用逼近理想解排序法从Pareto最优解选择单一最优解从而确定最适合的黏度模型及黏度模型参数；步骤S5的具体实施步骤为：

S501、构造m行n列的矩阵(x_ij)_m×n；其中，m为步骤S4中获得的Pareto最优解中解的个数，即达到终止条件时最新的父代种群P(t+1)的大小，n为目标函数个数，n＝2，x_ij为第i个解的第j个目标函数的目标函数值；

S502、将矩阵(x_ij)_m×n通过下式处理生成正则化矩阵(t_ij)_m×n：

S503、通过下式处理正则化矩阵(t_ij)_m×n得到权重正则化矩阵(a_ij)_m×n：

a_ij＝w_j×t_ij,i＝1,2,...,m,j＝1,2,...,n，

式中，w_j为第j个目标函数的相对权重，由下式获得：

S504、通过下两式确定正理想解A⁺和负理想解A^-：

A⁺＝[min(a₁₁,..,a_m1),min(a₁₂,..,a_m2),..,min(a_1n,..,a_mn)],

A^-＝[max(a₁₁,..,a_m1),max(a₁₂,..,a_m2),..,max(a_1n,..,a_mn)]；

S505、分别按下两式计算Pareto最优解中每个解距离正理想解A⁺和负理想解A^-的距离：

其中，

为Pareto最优解中第i个解距离正理想解A⁺的距离；

为多目标优化问题的正理想解；

为Pareto最优解中第i个解距离负理想解A^-的距离；

为多目标优化问题的负理想解；

S506、计算Pareto最优解中每个解距离正、负理想解的相对距离C_i：

S507、将Pareto最优解中的所有解距离正、负理想解的相对距离C_i进行排序，取让C_i最大的解为单一最优解A_final：A_final＝A∈max(C_i)，该单一最优解即为最终确定的模型参数。

2.根据权利要求1所述的纤维素生物墨水黏度特性预测方法，其特征在于，步骤S4的具体实施步骤为：

S401、确定种群大小N_p，并随机产生初始种群P(t)|t＝0作为父代种群；其中，采用双精度向量作为种群类型并规定种群大小N_p为200；

S402、使用遗传算法算子对父代种群P(t)依次进行选择操作，交叉操作和变异操作，以产生子代种群Q(t)；其中，选择操作采用锦标赛选择算法；交叉操作中交叉概率设为0.8；变异操作中变异概率设为0.02；

S403、将父代种群P(t)与子代种群Q(t)合并生成新的集合R(t)；

S404、对R(t)依次进行快速非支配排序和拥挤度计算以选出新的父代种群P(t+1)；

S405、判断是否达到终止条件：若达到，则最新的父代种群P(t+1)为Pareto最优解；若未达到，则重复S402～S404步骤，直至达到终止条件；其中，终止条件设置为最大代数为20000。

Images