CN111209657A - 考虑液体表面张力的固体变形界面计算方法 - Google Patents

Abstract

本发明提供了一种考虑液体表面张力的固体变形界面计算方法，直接考虑位移‑面力解耦数值方法，利用三次样条曲线对有限元模型中的流‑固边界面的离散点进行拟合，计算出每一点的曲率，得到液‑固界面的整体曲率数据，解决具有表面张力的流‑固耦合载荷问题。本发明的边界曲率计算与有限元分析交替进行，能够充分利用现有有限元程序的前处理、后处理、非线性本构关系等功能。

Description

考虑液体表面张力的固体变形界面计算方法

技术领域

本发明涉及固体变形界面计算领域，具体是一种考虑液体表面张力的固体变形界面计算方法。

背景技术

土壤、混凝土、生物组织属于含液多孔固体介质结构，这种具有流-固耦合属性的结构广泛应用于人们的生产实践中。在流体与固体的分界面，液体表面由于液体分子受力不均匀，促使液体表面收缩，形成表面张力。对于含液多孔固体介质结构，孔尺寸与固体的弹性模量都较小时，液体表面张力的作用对固体变形产生较为明显的影响。另外，液体表面张力对固体产生的面力大小与固体曲率相关。当固体变形时，固体的曲率发生变化，使液体对固体的面力变化，导致固体进一步改变形状，这是一个典型的力-位移耦合问题。

液体对固体的面力与表面张力、液-固界面的曲率有关，而对于同一种液体，其表面张力系数近似为常数。因此，只需要确定液-固界面的曲率就能够确定出液体对固体面力。采用分子动力学理论模拟液滴夹杂问题，能够在微观尺度下进行细致的分析，这种方法受到时间和空间尺度的限制，不适合计算大规模问题，不能够满足原子数目非常多的情况。利用现有的有限元商业软件，只能计算出每一个节点的位移，不能计算网格节点边界的曲率，不能够给出整体边界的曲率，无法计算流体对固体的面力。因此，需要进一步发展有限元计算方法，计算液-固界面的曲率，得到液体对固体面力，从而分析含液多孔结构介质结构的运动状态与力学性能。

发明内容

本发明为了解决现有技术的问题，提供了一种考虑液体表面张力的固体变形界面计算方法，直接考虑位移-面力解耦数值方法，利用三次样条曲线对有限元模型中的流-固边界面的离散点进行拟合，计算出每一点的曲率，得到液-固界面的整体曲率数据，解决具有表面张力的流-固耦合载荷问题。

一种考虑液体表面张力的固体变形界面拟合方法,其特征在于步骤如下：

步骤1：建立流-固耦合几何模型。

步骤2：利用有限元软件划分网格。

步骤3：将外载荷F，分为N步，初始载荷为0，载荷增量步△F＝F/N。

步骤4：提取流-固边界面节点坐标，设共有n个节点，坐标为(x₀,y₀)、(x₁,y₁)、……、(x_n,y_n)，根据公式(1)，进行三次样条曲线拟合。

其中

步骤4.1：在区间[0,1]中，生成n+1个点，为0,1/n，2/n，……，n-1/n,1。

步骤4.2：令t₀＝0,t₁＝1/n，t₂＝2/n，……，t_n-1＝n-1/n,t_n＝1。

步骤4.3：将x₀＝x(t₀),x₁＝x(t₁),…,x_n-1＝x(t_n-1),x_n＝x(t_n)和t₀,t₁,……,t_n-1,t_n带入公式(2)，得到n+1个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程。

其中

步骤4.4：将公式(2)求导，带入t₀与t_n,令x’(t₀)＝0,x’(t_n)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程。

步骤4.5：将步骤4.3与步骤4.4得到的n+3个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数，得到方程(2)的具体表达式。

步骤4.6：将y₀,y₁,…,y_n-1,y_n和t₀,t₁,……,t_n-1,t_n带入公式(3)，得到n+1个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程。

其中

步骤4.7：将公式(3)求导，带入t₀与t_n,令y’(t₀)＝0,y’(t_n)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程。

步骤4.8：将步骤4.6与步骤4.7得到的n个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数，得到方程(3)的具体表达式。

步骤4.9：三次样条曲线的参数方程为：

步骤4.10：按照公式(5)求解曲线的一阶导数，按照公式(6)求解曲线的二阶导数,按照公式(7)求解曲线的曲率k。

步骤5：按照公式(8)，求解各节点处的表面张力对固体产生的面力△P。

△P＝kγ (8)

γ为液体表面张力系数。

步骤6：将载荷△P与载荷增量的总载荷，施加在有限元模型上，进行分析计算，得到新的流-固边界面坐标。

步骤7：流-固耦合边界面的位置收敛，则外载荷F增加一个增量步；若不收敛，则返回步骤4，进行迭代。

步骤8：载荷大于总载荷F，则有限元计算结束，得出模型最终平衡状态位置；载荷小于或等于总载荷，返回步骤4，进行迭代。

本发明的有益效果在于：

1.通过三次样条曲线对流-固边界面离散点进行拟合，能够计算出液-固边界面的曲率，确定出液体对表面张力作用对固体的面力。

2.边界曲率计算与有限元分析交替进行，能够充分利用现有有限元程序的前处理、后处理、非线性本构关系等功能。

附图说明

图1是考虑液体表面张力的固体变形界面拟合方法流程图；

图2是加载前有限元网格模型；

图3是表面张力引起的面力与外载荷共同作用后，变形后网格图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步说明。

实施例1：

步骤1：建立流-固耦合几何模型。

步骤2：利用有限元软件划分网格，如图2。

步骤3：将外载荷F，分为N步，初始载荷为0，载荷增量步△F＝F/N。

步骤4：提取流-固边界面节点坐标，设共有40个节点，坐标为(x₀,y₀)、(x₁,y₁)、……、(x₄₀,y₄₀)，根据公式(1)，进行三次样条曲线拟合。

其中

步骤4.1：在区间[0,1]中，生成40+1个点，为0,1/40，2/40，……，39/40,1。

步骤4.2：令t₀＝0,t₁＝1/40，t₂＝2/40，……，t_39-1＝39/40,t₄₀＝1。

步骤4.3：将x₀＝x(t₀),x₁＝x(t₁),…,x₃₉＝x(t₃₉),x₄₀＝x(t₄₀)和t₀,t₁,……,t₃₉,t₄₀带入公式(9)，得到41个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₈，β₃₉参数的方程。

其中

步骤4.4：将公式(9)求导，带入t₀与t₄₀,令x’(t₀)＝0,x’(t₄₀)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₈，β₃₉参数的方程。

步骤4.5：将步骤4.3与步骤4.4得到的43个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₈，β₃₉参数值，从而得到方程(9)的具体表达式。

步骤4.6：将y₀,y₁,…,y₃₉,y₄₀和t₀,t₁,……,t₃₉,t₄₀带入公式(10)，得到41个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₉，β₄₀参数的方程。

步骤4.7：将公式(10)求导，带入t₀与t₄₀,令y’(t₀)＝0,y’(t₄₀)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₈，β₃₉参数的方程。

步骤4.8：将步骤4.6与步骤4.7得到的43个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β₃₈，β₃₉参数值，从而得到方程(10)的具体表达式。

步骤4.9：三次样条曲线的参数方程为：

步骤4.10：按照公式(5)求解曲线的一阶导数，按照公式(6)求解曲线的二阶导数,按照公式(7)求解曲线的曲率k。

步骤5：按照公式(8)，求解各节点处的表面张力对固体产生的面力△P。

△P＝kγ (8)

γ为液体表面张力系数。

步骤6：将载荷△P与载荷增量的总载荷，施加在有限元模型上，进行分析计算，得到新的流-固边界面坐标。

步骤7：流-固耦合边界面的位置收敛，则外载荷F增加一个增量步；若不收敛，则返回步骤4，进行迭代。

步骤8：载荷大于总载荷F，则有限元计算结束，得出模型最终平衡状态位置，如图3；载荷小于或等于总载荷，返回步骤4，进行迭代。流程图如图1。

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (2)

1.一种考虑液体表面张力的固体变形界面计算方法，其特征在于：步骤如下：

步骤1：建立流-固耦合几何模型；

步骤2：利用有限元软件划分网格；

步骤3：将外载荷F，分为N步，初始载荷为0，载荷增量步△F＝F/N；

步骤4：提取流-固边界面节点坐标，设共有n个节点，坐标为(x₀,y₀)、(x₁,y₁)、……、(x_n,y_n)，根据以下公式，进行三次样条曲线拟合，并计算节点处曲率k，

其中

步骤5：求解各节点处的表面张力对固体产生的面力△P，

△P＝kγ； (8)

γ为液体表面张力系数；

步骤6：将载荷△P与载荷增量的总载荷，施加在有限元模型上，进行分析计算，得到新的流-固边界面坐标；

步骤7：流-固耦合边界面的位置收敛，则外载荷F增加一个增量步；若不收敛，则返回步骤4)，进行迭代；

步骤8：载荷大于总载荷F，则有限元计算结束，得出模型最终平衡状态位置；载荷小于或等于总载荷，返回步骤4)，进行迭代。

2.根据权利要求1所述的考虑液体表面张力的固体变形界面计算方法，其特征在于，步骤4)所述的曲线拟合及节点处曲率k计算过程如下：

步骤4.1：在区间[0,1]中，生成n+1个点，为0,1/n，2/n，……，n-1/n,1；

步骤4.2：令t₀＝0,t₁＝1/n，t₂＝2/n，……，t_n-1＝n-1/n,t_n＝1；

步骤4.3：将x₀＝x(t₀),x₁＝x(t₁),…,x_n-1＝x(t_n-1),x_n＝x(t_n)和t₀,t₁,……,t_n-1,t_n带入以下公式，得到n+1个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程，

其中

步骤4.4：将公式(2)求导，带入t₀与t_n,令x’(t₀)＝0,x’(t_n)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程；

步骤4.5：将步骤4.3与步骤4.4得到的n+3个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数，得到方程(2)的具体表达式；

步骤4.6：将y₀,y₁,…,y_n-1,y_n和t₀,t₁,……,t_n-1,t_n带入以下公式，得到n+1个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程，

其中

步骤4.7：将步骤4.6)得到的参数方程进行求导，带入t₀与t_n,令y’(t₀)＝0,y’(t_n)＝0,得到2个关于α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数的方程；

步骤4.8：将步骤4.6)与步骤4.7)得到的n个方程进行求解，得到α₁，α₂，α₃，α₄，β₁，β₂，……，β_n-2，β_n-1参数，得到4.6)参数方程的具体表达式；

步骤4.9：三次样条曲线的参数方程为：

步骤4.10：按照公式(5)求解曲线的一阶导数，按照公式(6)求解曲线的二阶导数,按照公式(7)求解曲线的曲率k，

Images