考虑微观变形特征的水凝胶均匀和非均匀溶胀分析及其影响参数研究1)

刘 岩 王惠明,†,2)

∗(浙江大学工程力学系,杭州 310027)

†(浙江省软体机器人与智能器件研究重点实验室,杭州 310027)

引言

聚合物网络是由聚合物链交联构成的一种复杂三维网状结构,在外界激励下可能产生超大变形.由于凝胶具有良好弹性、生物相容性和生物降解性,被广泛应用于组织工程[1]、药物运输[2]和细胞培养[3]等领域.浸没在溶剂中的聚合物吸收大量溶剂分子形成凝胶并伴随发生溶胀[4].溶胀过程受化学势的影响[5],当凝胶内外溶剂分子化学势相等且满足力学平衡时凝胶达到平衡状态,此时网络不再吸纳溶剂分子,进入饱和状态.凝胶的溶胀变形同时也受环境温度、pH 值和网络交联密度等因素影响[6-7].

关于凝胶的溶胀响应分析,通常采用基于自由能理论的模型,包括两部分,一部分是聚合物网络被拉伸产生的自由能,另一部分是溶剂分子与网络混合产生的混合能.对于聚合物网络被拉伸产生的自由能部分,基于宏观唯象模型的有neo-Hookean 模型[8]、Money-Rivlin 模型[9]、Gent 模型[10]和Ogden 模型[11]等.从微观角度,认为聚合物变形响应主要源于单链构形熵的变化,其理论模型包括三链模型[12](认为空间内单链主要分布在3 个相互垂直的方向),四面体模型[13-14](4 条链交联在正四面体中心点),八链模型[15](单链分布在六面体对角线)以及全网络模型[16](考虑空间中单链具有随机取向),其中三链模型认为聚合物在变形过程中宏观网络与微观单链变形是仿射关系(即宏观网络和微观单链变形相等),而四面体模型和八链模型认为宏观与微观变形是非仿射关系.

采用仿射变形假设的全网络模型预测能力介于三链模型和八链模型之间[17].当凝胶进入大变形阶段后,主方向上单链的伸长率逐渐接近极限值,随着变形的进一步增大,主方向上的单链伸长率增大量很小,但是其他方向上的单链还没有接近极限伸长率,因此采用仿射变形假设预测的其他方向上的单链伸长率就会低于实际伸长率[18].Miehe 等[19]采用非仿射变形假设的全网络模型,同时考虑了单链的圆管约束,很好地模拟超弹性材料的单轴拉伸、等双轴拉伸和剪切变形的实验结果,Huang 等[20]系统地总结了各种智能凝胶的本构模型.Li 等[21]基于仿射变形假设的全网络模型,通过聚合物网络被拉伸时单链伸长率的变化预测了水凝胶的断裂行为.Lei 等[22]构建了一个介观尺度模型研究水凝胶的大变形和断裂机理.

利用超弹性材料自由能模型[23-24](如 neo-Hookean 模型和八链模型等),可开展软材料结构的变形预测[25],结合Flory 提出的混合能模型[26],可进行溶胀引起的变形分析和预测[27-29].用neo-Hookean 模型和Flory 混合能模型可以较准确地预测凝胶中小变形情况,当凝胶进入大变形阶段,neo-Hookean 模型预测能力受到限制并且不能分析单链变形情况.采用单链末端距服从Langevin 分布的模型可以解决这一问题,Cohen 等[30]采用仿射变形假设的全网络模型,根据Langevin 函数计算单链有限伸长,然后利用数值微球方法将单链集合成网络[31],最后获得凝胶发生溶胀变形时单链和网络的力学响应,该方法已经被广泛应用[32-37].

本研究利用非仿射变形的全网络模型,假设每个单链变形受到由周围链形成的类似圆管状的约束,开展凝胶发生溶胀变形时微观单链和宏观网络的变形分析.为了更直观地解释单链变形,分析凝胶在三种溶胀状态下网络和单链的变形情况,第一种溶胀状态为自由溶胀,凝胶没有外力和约束作用,放入给定压力溶剂中溶胀至平衡状态;第二种溶胀状态为单一方向溶胀,首先将立方体凝胶进行等双轴预拉伸,然后放入装有给定压力溶剂的容器中,使其沿单一方向溶胀.这两种状态下,凝胶和单链都产生均匀变形(在各个位置的单链变形相同).第三种溶胀状态为具有刚性核的球形凝胶溶胀,将球形凝胶放入具有一定压力的溶剂中使其发生溶胀.本文给出的模型可以模拟凝胶在溶胀变形情况下具有确定链段数的单链的有限伸长以及网络变形情况,更好地理解影响溶胀变形的微观机理,为凝胶性能的微观设计提供指导.

1 自由能函数

将聚合物网络中没有水分子进入的状态设置为参考状态.考虑一个不可压缩的聚合物网络,单位参考体积内含有N条具有随机取向且均匀分布的单链,每条单链由n个长度为l的链段组成[38].聚合物网络占据的体积为Vp=nNvp,其中vp是一个链段占据的体积.单链构形示意图如图1(a)所示.参考状态下,每个单链的平均末端距r0=若单链完全伸展,末端距L=nl.在图1(b)中,每个单链受到由周围单链构成的半径为d0的类似圆管状约束.

将参考构型中任一物质点的位置矢量用X表示,当网络发生变形时,物质点会移动到新的位置,相应的位置矢量为x(X).定义网络的变形梯度为F=∂x(X)/∂X.将聚合物放入压力为p的溶剂中,使其溶胀至平衡状态,假设有m个溶剂分子进入网络,此时凝胶的体积就等于Vp+mvl,其中vl是一个液体分子的体积.假设液体分子和链段均是不可压缩的,且在三维空间中占据m+nN个晶格位置,每个晶格的体积为vc.为了简化计算,认为链段、液体分子和晶格体积是相等的,记为vp=vl=vc[30],每个单元的体积变化率J=detF=1/cp,其中cp=nN/(m+nN)是聚合物网络的体积分数,溶剂分子的体积分数记为cl=1 −cp.

在变形状态下,定义每个单链的伸长率和相对伸长率分别为λ=r/r0和λrel=r/L,其中r为变形状态下单链的末端距,则有和λrel∈[0,1).每条单链受到由周围单链构成的类似圆管状约束[39],这样的几何形状约束了单链构形数.将每条单链受到的约束用无量纲量g=(d0/d)2表示,其中d是当前状态下圆管的直径,d0是参考状态下圆管的直径,关于长度为l的Kuhn 链段约束直径见文献[40].

当其发生溶胀变形时,网络的构形熵会降低,但网络和溶剂分子形成的混合物的熵增加[41].根据统计力学,每个单链的自由能由式(1)给出

式中,T是绝对温度,K是玻尔兹曼常数,Ω 是单链的构形数密度.整个网络的自由能等于所有单链自由能的总和,由式(2)计算

式中,Wm是网络与溶剂分子混合过程产生的自由能,Π 是拉格朗日乘子.第三项是拉格朗日乘子与(溶剂分子和聚合物)体积不可压缩性条件的乘积.拉格朗日乘子解释为网络和溶剂分子混合物的总压力,第四项是静水压力对凝胶做功.

由于neo-Hookean 模型没有考虑单链的有限拉伸,因此该模型普遍适用于中小变形的预测,当凝胶发生大变形时,neo-Hookean 模型的预测能力将会降低.本研究采用Langevin 模型分析单链有限伸长以及预测大变形情况.假设网络中的每个单链都服从Langevin 分布,其构形数概率密度[19]

式中,a0是无量纲参数,α 是与管截面形状有关的修正因子[42],其中β=ℓ−1(λr),ℓ(β)=coth β −1/β 是Langevin 函数,逆Langevin 函数可用Pad´e估计值[43]

由式(1),每个单链自由能可表示为

式中,c=KTlna0是常数,U=α(l/d0)2称为有效圆管几何参数[19].

聚合物网络发生溶胀变形时,其体积会逐渐增大,因此网络变形产生的自由能需要考虑其体积变化,将式(2)中添加体积变化引起的自由能,得到网络自由能函数[18]

网络和溶剂分子混合过程产生的混合能通过自由能函数表示[4]

式中,χ 是凝胶网络和溶剂分子相互作用参数.凝胶中溶剂分子化学势为

当n→∞时的结果与Hong 等[5]结果一致.当凝胶达到溶胀平衡时,凝胶内外溶剂分子的化学势相等,此时凝胶处于动态平衡状态,即溶剂分子进出凝胶时刻保持相等.

将外部溶剂的参考状态设置为液体,此时液体在室温下与自身气体状态达到平衡.如果将自身气体状态下的压力设置为p0,该压力可视为气液相转化的临界值.若当前状态外部溶剂压力p小于p0,溶剂分子的化学势等于µl=KTln(p/p0),若当前状态外部溶剂压力p大于p0,溶剂分子化学势µl=(p−p0)vc[44].本研究考虑外部溶剂压力p/p0≥1的情况.

2 宏观与微观变形关系

凝胶发生自由溶胀时,宏观网络与微观单链变形相同.由于初始状态下微观单链蜷缩在凝胶网络内部,凝胶在外力或位移约束情况下发生溶胀变形时,网络与单链变形将会不一致.为了改善仿射变形假设的模型对超弹性材料变形的预测能力,放松在仿射变形中认为单链伸长率λ 与网络伸长率¯λ 相等的假设,引入定义在单位球上波动场f,将单链伸长率和网络伸长率间的关系表示为[19]

式中,〈〉h是h次方根平均值算子,有

最终得到非仿射变形单链的伸长率与宏观网络之间的关系[19]

用于计算式(13)的数值微球方法介绍如下.

认为空间中所有的单链具有随机取向的特性且均匀分布在微球内,单链在球面上每个位置处的概率密度为1/4π,将式(13)的积分计算转化为在球面上m个离散点处函数值的代数和.即

式中,a是微球内状态变量,w(s)是离散点处的权重,满足

研究表明,通常采用42 个离散点以及相应权重进行计算可以得到足够精确的结果[31].

根据自由能函数(3)确定聚合物网络的名义应力为

如果式(15)中,h=2,且在不考虑单链圆管约束的情况下,即q=0 或U=0,则该应力表达式等价于八链模型的应力表达式[15].如果假设单链与网络变形为仿射关系(即λ=),应力表达式等价于Cohen等[30]采用模型的应力表达式.

3 数值结果

考虑3 种不同情况的凝胶溶胀时的宏观和微观变形:(1)凝胶发生自由溶胀,此时产生均匀各向同性变形;(2)凝胶被放置在容器中,只能沿某一个方向发生溶胀,此时凝胶产生均匀但在各个方向上不同的变形;(3)内表面固定的球形凝胶浸没在溶剂中的溶胀行为,此时凝胶产生非均匀的且在各个方向上不同的变形.

将聚合物放入水溶液中使其溶胀,在下面的计算中,材料及环境参数取为[30]:KT=4.1× 10−21J,每个水分子的体积vl=3.0×10−29m3.室温下,水的饱和蒸气压p0=3 kPa,相互作用参数χ=0.5.除非特殊说明,计算时取聚合物网络中单链密度N=2.5×1025chains/m3,每条单链含有链段数n=100.

3.1 自由溶胀

选择尺寸为L×L×L立方体凝胶,内部含有随机取向且均匀分布的单链,将其浸没在溶剂外部压力为p的溶剂中发生自由溶胀.此时,凝胶在3 个相互垂直方向上的伸长率相同λ1=λ2=λ3=λfree.达到平衡状态时,凝胶内外化学势相等,同时满足力学边界条件σ=−pI.当单链完全伸展开,其伸长率λmax=此时凝胶达到最大可能体积变形,因此可以控制单链的链段数调控凝胶最大体积变形.

不考虑单链约束时,发生自由溶胀时凝胶各个方向伸长率相等,由式(13)积分结果

发现单链与网络变形关系无论是非仿射还是仿射,采用八链模型和全网络模型(仿射或非仿射变形)得到的真实应力均为

Cohen 等[30]采用仿射变形假设得到的自由溶胀情形的真实应力表达式为

图2 给出了凝胶在不同压力溶剂中达到平衡时的伸长率和体积变化的曲线,本研究所采用的模型的结果与Cohen 等[30]采用仿射变形假设计算得到的结果相一致,并且在自由溶胀情形可以覆盖采用仿射变形假设得到的结果.从图2 中可以发现,随着外部溶剂压力增大,进入网络中的溶剂分子越来越多,导致凝胶体积逐渐增大并且呈非线性增大趋势.当压力p/p0≥40 时,单链趋近于极限伸长即完全伸展状态,此时凝胶网络容纳液体分子能力趋近于极限值,即接近饱和状态.

图2 单链无约束凝胶自由溶胀伸长率和体积变化Fig.2 Stretch and the normalized volumetric ratio as a function of the normalized external solvent pressure in a freely swelling hydrogel with unconstrained polymer chain

在真实凝胶中,单链会受到周围链的约束,使其不能在凝胶内部任意移动.考虑单链约束时,有效圆管几何参数U会对自由溶胀产生影响.采用圆管约束模型,可得凝胶真实应力表达式

若凝胶具有不同有效圆管几何参数U,发生自由溶胀时会产生不同的力学响应.图3 表示q=1 条件下,参数U对自由溶胀变形的影响,参数U越大,平衡时单链和网络的伸长率越大.由式(18)发现表达式中不含有参数h,因此参数h不会对自由溶胀产生影响.在不考虑单链约束影响时,浸没在溶剂中的所有材料达到平衡时具有相同溶胀率.当溶剂压力p/p0≥20 时,参数U的影响基本消失,主要原因是单链伸长已经接近完全伸展状态.

图3 具有不同有效圆管几何参数U 的凝胶在溶胀至平衡状态时的体积随溶剂压力的演化规律Fig.3 Normalized volumetric ratio as a function of the normalized external solvent pressure for different values of the effective tube geometry parameter U

凝胶发生自由溶胀变形时,单链的链段数和单链的密度对网络和单链的伸长率也有一定的影响.从图4 中可以看出,当单链的密度N一定时,随着单链的链段数n的增加,凝胶网络和单链的伸长率也随之增大,这是因为构成单链的链段数越多,其长度越长,导致单链的蜷缩程度越高.单链的密度对凝胶溶胀影响也比较明显,单链的链段数n一定时,随单链密度的增大,网络和单链的伸长率先快速下降,然后趋于恒定值,造成这一结果的原因是单链的密度越高(即凝胶网络交联点个数越多),网络中可以容纳溶剂分子的空间越小使得能够进入网络的溶剂分子越少,因此网络溶胀变形越小.

图4 溶胀至平衡状态时的伸长率随单链密度和单链链段数的变化Fig.4 Effects of chain density and the number of polymer segments on the swelling-induced stretch

3.2 单一方向溶胀

首先,将尺寸为L×L×L的立方体聚合物先沿X1和X3方向等双轴拉伸至伸长率λb,使其具有新的尺寸,如图5 所示.然后,将预拉伸后的聚合物放入尺寸为λbL×H×λbL容器中,这里.最后,将容器中放入压力为p的液体使聚合物发生溶胀变形,不考虑容器内表面的摩擦,凝胶只能沿X2方向溶胀.达到平衡时,有λ1=λ3=λb,σ2=−p.

图5 凝胶预拉伸模型Fig.5 Hydrogel cube with biaxial pre-stretch

图6 给出无预拉伸凝胶溶胀平衡时其网络和单链的伸长率和横向应力σ11随标准化压力p/p0的变化.这种情形可具体操作为,例如将横截面尺寸为L×L的凝胶置于内表面边长为L×L的容器中进行溶胀.当外部溶剂压力增加至p/p0=4000,单链基本达到完全伸展状态.注意到X1和X3方向应力一直都是压应力,这是由于X1和X3方向的溶胀变形受到容器壁的限制而导致的.

图6 无预拉伸(λb=1)凝胶在受限溶胀至平衡时网络和单链的伸长率以及应力σ11随压力的变化Fig.6 Stretches of network and chains and the normalized stress as a function of the normalized external solvent pressure in a constraint swelling hydrogel cube without pre-stretch(λb=1)

当对凝胶施加一定预拉伸时,如图7 所示,凝胶出现单链变形和网络变形相等的溶胀平衡状态,即λ2=λchain=λb,此时X1方向的应力(σ11+p)/p0=0.这是因为对凝胶施加预拉伸时,使网络中单链向三维空间伸展.同时,被拉伸的网络有收缩趋势,对溶剂分子产生压力.在外部溶剂压力较小时,网络中预拉伸产生的应力大于吸收溶剂分子发生体积膨胀引起的应力,X1方向为拉应力,即(σ11+p)/p0>0.随着溶剂压力增大,越来越多溶剂分子进入网络,达到平衡时网络中预拉伸产生的应力小于体积膨胀产生的应力,X1方向出现压应力,即(σ11+p)/p0<0.因此,在某一个溶剂压力作用下一定存在预拉伸产生的应力等于体积膨胀产生的应力,即图7 中A点,对应的压力p/p0≈7.9,伸长率为λ2=λb=8.

图7 预拉伸(λb=8)凝胶在溶胀至平衡时其伸长率和应力随压力的变化Fig.7 Stretch and the normalized stress as a function of the normalized external solvent pressure in a constraint swelling hydrogel cube with biaxial pre-stretch(λb=8)

从图6 和图7 中可以看到,单链极限伸长率不会超过10,极限伸长率主要取决于构成单链的链段数.根据图6 单链变形曲线发现,无预拉伸的凝胶中单链变形均小于网络变形,其主要原因是模型中引入了圆管约束,例如单链的相互缠结就是约束链运动的一种形式.而从图7 中看到,具有预拉伸的凝胶在溶剂压力p/p0≈7.9 时,单链变形等于网络变形;如果溶剂压力p/p0≤7.9,单链变形大于网络变形,主要原因是预拉伸使得网络中绻缩的单链具有初始伸长率,将凝胶放入溶剂中发生溶胀,单链的伸长率继续增加.

图8 给出了参数h对无预拉伸凝胶平衡状态下伸长率与溶剂压力关系曲线的影响.从图8 中可以发现,随着h的增大,单链伸长率与网络伸长率λ2间的差异逐渐减少.当h=50 时,若溶剂压力p/p0≤1500,溶胀平衡状态下单链伸长率接近于网络伸长率.这也表明了材料参数h越大的凝胶发生溶胀变形时,单链伸长率和网络伸长率差异越小.表明h是调控单链伸长率和网络伸长率间关系的重要参数.

3.3 球形凝胶溶胀

考虑含有刚性核的球形凝胶的溶胀,这是非均匀溶胀情形.在参考状态下,球形凝胶的外半径和内半径之比B/A=2,其中A和B分别是球形凝胶的内半径和外半径.在凝胶内部有一个半径为λ0A的刚性核,刚性核和凝胶在界面处完美结合,即位移和应力均连续,然后让其溶胀至平衡状态.

图8 参数h 对无预拉伸凝胶平衡状态下伸长率与溶剂压力关系曲线的影响Fig.8 Effect of material parameter h on the swelling-induced stretch in a constraint swelling hydrogel cube

球形凝胶发生球对称受限溶胀时,网络变形是非均匀的且径向和环向变形是不一致的,各物理场均是位置坐标的函数.当凝胶达到平衡时,凝胶内外溶剂化学势相等µl=(p−p0)vc,平衡方程简化为

名义应力的表达式为

式中λr=dr/dR和λθ=r/R分别为径向伸长率和环向伸长率.溶剂分子体积分数为边界条件为

下文计算中设定p/p0=1.5,n=100,h=1.5,q=1.5,U=2,Nvc=7.5×10−4,χ=0.5.图9 给出了凝胶溶胀至平衡状态下网络的径向和环向名义应力沿径向的分布.图中表明,径向名义应力均为拉应力,环向名义应力均为压应力.由内向外,网络的径向和环向名义应力的幅值逐渐减小.图10 给出了平衡时凝胶网络和单链伸长率的沿径向的分布.从图10 可以看出,网络的径向和环向伸长率均大于1,表明网络在径向和环向均是受到拉伸.凝胶沿径向各点处的网络径向伸长率大于环向伸长率.不同位置单链伸长率均有=10,即各单链均没有达到完全伸展.在内表面附近,单链伸长率高于自由溶胀时的伸长率,但小于网络的径向伸长率,这是由于内表面刚性核的固定约束使得出现径向拉应力的缘故.随着坐标位置远离内表面,内表面的位移约束影响逐渐减小,单链伸长率接近网络伸长率,且接近于自由溶胀时的网络伸长率.

图9 凝胶溶胀至平衡状态时的径向和环向应力沿径向分布Fig.9 Distributions of the stresses in a swollen spherical hydrogel with a rigid core

图10 凝胶溶胀平衡状态时的宏观和微观伸长率沿径向分布Fig.10 Distributions of the stretches of micro-chain and macro-network in a swollen spherical hydrogel with a rigid core

图11 和图12 分别给出了U=2,6 和12 时凝胶中溶剂分子的体积分数和网络中的渗透压沿径向的分布.由图11 可以发现靠近刚性核位置,溶剂分子体积分数减小,这是因为内表面固定使得靠近内表面网络渗透压增大(如图12 所示),溶剂分子不易进入网络.从图12 中可以看出,网络中的渗透压随参数U的增大而降低.当R/A>1.2 时,刚性核的影响逐渐消失,溶剂分子体积分数的分布趋于均匀,且趋于自由溶胀状态的数值.溶剂分子的体积分数随参数U的增大而增大,即单链的链段长度与圆管约束半径比值越大,网络中的溶剂分子体积分数越大,最终网络和单链的伸长率越大.网络的渗透压会影响网络中溶剂分子的体积分数,从图11 和图12 可以看出,二者是一种负相关关系.

图11 凝胶溶胀平衡状态溶剂分子体积分数沿径向的分布Fig.11 Effect of parameter U on the distributions of the volume fraction of solvent molecules in a swollen spherical hydrogel with a rigid core

图12 凝胶溶胀平衡状态渗透压沿径向的分布Fig.12 Effect of parameter U on the distributions of the osmotic pressure in a swollen spherical hydrogel with a rigid core

另取计算参数n=90,h=2.0,q=0,U=0,采用非仿射全网络模型计算得到的内含刚性核球壳凝胶溶胀至平衡状态时的径向和环向应力沿径向分布曲线如图13 的实线所示.图13 中的离散点是Zhao等[27]采用neo-Hookean 模型的计算结果,结果表明,选取合适的参数时,本研究的模型可以取得与采用neo-Hookean 模型基本一致的结果.

4 结论

基于全网络模型,假设凝胶网络和单链变形为非仿射关系,利用聚合物网络的单链受到由于周围链的作用而产生类似圆管状的约束的模型,建立了可以预测溶胀引发的凝胶网络和单链变形的力学分析模型.开展了三种不同情形的凝胶溶胀性能分析,即立方体凝胶的自由溶胀和两个方向受限的单一方向溶胀以及具有刚性核的球形凝胶的非均匀溶胀.

图13 全网络模型与neo-Hookean 模型计算结果的对比Fig.13 Comparison between the presented model and the neo-Hookean model

自由溶胀时,凝胶网络伸长率和单链伸长率一致.伸长率随着参数U的增大而增大.随着溶剂压力的增大,网络伸长率和单链伸长率均趋近于渐近值两个方向受限的单一方向均匀溶胀情形,此时凝胶网络伸长率和单链伸长率表现出不一致.网络伸长率和单链伸长率可通过预拉伸进行调整.特别是,预拉伸增大到一定值时,可出现单链变形和网络变形相等的溶胀平衡状态.具有刚性核的球形凝胶发生溶胀时,凝胶和单链均为非均匀变形,在靠近刚性核的区域,出现较大的径向拉应力和环向压应力,这个区域对应的单链伸长率也较大.靠近外表面区域的网络和单链伸长率接近自由溶胀伸长率.网络的渗透压随参数U的增大而降低,溶剂分子体积分数随参数U的增大而增大.

所提出的模型可用于分析和预测凝胶溶胀时微观单链的变形情况,进而可通过调整单链微观结构来控制凝胶网络的溶胀变形,为凝胶器件的设计和性能优化提供理论基础.