1 力链的概念

1.1 强力与弱力

光弹试验和数值模拟研究均表明颗粒材料中的接触力网络并非呈均匀分布状态。通常，由相对稀疏的强力（大于平均接触力）网络来传递主要的荷载，而大多数的接触（或颗粒）则构成弱力网络（如图1所示）。即使是在外界荷载和内部结构都处于各向同性状态，部分接触仍传递着几倍于平均接触力的接触力。而在剪切过程中，接触力的非均匀分布现象尤为显著，较大接触力的方向将会逐渐趋向于大主应力方向，并在一定条件下出现柱状结构。

图1 颗粒材料中的接触力网络

1.2 力链 (Force chain)

相关研究表明强、弱力对颗粒材料系统差异迥然的贡献：应力张量的偏应力全部来源于强力子网络，而弱力子网络则只构成球应力部分。为了进一步描述强力子网络的特征，许多学者提出了力链（Force chain）的概念。力链不同于接触力分布、组构张量和配位数等信息，它是应力传递路径，能更深刻地揭示颗粒材料的结构层次与宏观力学性质的关系，并提出了颗粒-力链-体系的多尺度研究框架。

力链的定量描述仍非常困难，因此本推文仅按照依据于接触力大小绘制力链。

2 Paraview与vtk简介

2.1 Paraview

ParaView是一个开源的科学可视化软件，‌专门用于可视化大规模数据集。‌ 它支持多种数据格式和平台，‌并提供丰富的可视化功能。‌ParaView的主要特点包括：‌

数据格式支持：‌ParaView能够处理多种数据格式，‌如VTK、‌NetCDF、‌HDF5等，‌适用于各种类型的科学数据集。‌
平台支持：‌它可以在Windows、‌Linux和Mac OS等多个平台上运行，‌并支持集群计算，‌能够在分布式计算环境中进行大规模数据集的可视化分析。‌
可视化功能：‌提供2D和3D可视化、‌流线图、‌切片、‌等值面、‌矢量图等，‌用户可以通过交互式界面进行数据的探索和分析。‌
可扩展性：‌具有良好的可扩展性，‌可以通过编写插件来增加新的功能和数据源。‌

大家可以在官网免费下载Paraview软件，相关的入门视频可以参考B站。

2.2 vtk

Visualization Toolkit (VTK) 是一个开源的免费软件系统，‌主要用于三维计算机图形学、‌图像处理和科学计算可视化。‌VTK提供了一个丰富的功能集，‌包括数据处理、‌图形生成、‌可视化、‌交互式操作和导出等。‌该系统的核心构建块包括数据对象、‌过程对象和过滤器，‌这些组件共同协作以实现复杂的数据处理和可视化任务。‌

数据对象：VTK提供了多种数据对象来表示不同类型的数据，包括点数据、网格数据、图像数据和多边形数据等。这些数据对象可以存储和处理各种属性和标量值。
过程对象：VTK中的过程对象是用于处理和操作数据的基本单元。它们可以执行数据的转换、过滤和操作，如数据剖析、几何变换、数据插值和数据重采样等。
过滤器：过滤器是VTK中用于处理数据的高级组件。它们可以将多个过程对象连接在一起，形成一个数据处理的流水线。每个过滤器接收输入数据，并生成输出数据。通过连接不同的过滤器，可以实现复杂的数据处理和可视化任务。

2.3 力链的vtk文件格式

传统VTK文件格式包含五个基本部分：

第一部分是文件版本和标识符。这部分包含单行# vtk DataFile Version x.x ，该行必须与显示的完全相同，但版本号x.x除外，该版本号会因VTK的不同发行版而变化，本推文利用5.1版本
第二部分是标题。标题由以行尾字符\n终止的字符串组成。最大为256个字符，可用于描述数据并包括任何其他相关信息。
下一部分是文件格式。文件格式描述文件的类型，可以是ASCII或二进制。在此行上，必须出现单个单词ASCII或BINARY。一般推荐ASCII格式。
第四部分是数据集结构。几何部分描述了数据集的几何和拓扑。这部分以包含关键字DATASET的行开头，后跟描述数据集类型的关键字。然后，根据数据集的类型，其他关键字/数据组合将定义实际数据。
最后一部分描述了数据集属性。这部分以关键字POINT_DATA或CELL_DATA开头，后跟一个整数，分别指定点或单元的数量。(先出现POINT_DATA还是CELL_DATA都没关系)。然后，其他关键字/数据组合定义实际的数据集属性值（即标量，向量，张量，法线，纹理坐标或字段数据）。

文件格式的概述如图2所示。前三个部分是必需的，而另两个部分是可选的。需要注意的是关键字不区分大小写，可以用空格分隔。

图2 VTK文件格式

针对本文的力链绘制，vtk的格式可以表示为：

# vtk DataFile Version 5.1
vtk output
ASCII
DATASET POLYDATA
POINTS {number of contacts} float64
pos_x1 pos_y1 pos_z1
pos_x2 pos_y2 pos_z2
……记录接触的坐标
……
POINT_DATA {number of contacts}
VECTORS directions float64
force_x1 force_y1 force_z1
force_x2 force_y2 force_z2 
……记录接触力向量

3 代码实现

3.1 生成vtk文件

根据2.3节中介绍的vtk格式,可以首先书写一个write_file_ini()函数，将vtk前几行文本写入vtk文件，注意到需要知道接触的个数，因此可以利用get_contact_num函数得到接触的个数。

def write_file_ini(self):
    with open(self.vtk_file_n, 'a') as file:
        file.write("# vtk DataFile Version 5.1\n")
        file.write("vtk output\n")
        file.write("ASCII\n")
        file.write("DATASET POLYDATA\n")
    
def get_contact_num(self):
    self.contact_num = it.contact.count()

随后利用 write_vtk(filename = '')函数将剩余的接触位置以及接触力向量写进vtk文件。最后将这些函数封装为类，得到最后的代码为：

# -*- coding: utf-8 -*-
"""
Created on Tue Aug 13 15:57:18 2024

@author: mq
"""

import itasca as it


class Force_Chain():
    
    def __init__(self):
        self.vtk_file_n = 'force chain.vtk'
        self.write_file_ini()
    
    def write_file_ini(self):
        with open(self.vtk_file_n, 'a') as file:
            file.write("# vtk DataFile Version 5.1\n")
            file.write("vtk output\n")
            file.write("ASCII\n")
            file.write("DATASET POLYDATA\n")
    
    def get_contact_num(self):
        self.contact_num = it.contact.count()
    
    def write_vtk(self, savfile):
        it.command("model restore '{}'".format(savfile))
        self.get_contact_num()
        with open(self.vtk_file_n, 'a') as file:
            file.write("POINTS {} float64\n".format(self.contact_num))
            for cp in it.contact.list():
                file.write('{} {} {}\n'.format(cp.pos_x(), 
                                                 cp.pos_y(), 
                                                 cp.pos_z()))
            
            file.write("POINT_DATA {}\n".format(self.contact_num))
            file.write("VECTORS directions float64\n")
            
            for cp in it.contact.list():
                file.write('{} {} {}\n'.format(cp.force_global_x(), 
                                                 cp.force_global_y(), 
                                                 cp.force_global_z()))

举例来说，例如我有一个PFC运行的jieguo_5.sav文件，我想得到力链分布，则可以调用：

from force_chain import Force_Chain
test = Force_Chain()
test.write_vtk('jieguo_5')

3.2 Paraview后处理

利用Paraview打开force chain.vtk文件，按照如下顺序执行，则可以得到如图3所示的力链分布图。

上方菜单栏Filters -> Alphabetical -> Glyph
Orientation Array选择directions. Scale Array选择directions by Magnitude. scale factor 选择 0.002.
Glyph Mode 选择 All Points
Glyph Type 选择 Line
上方菜单栏Filters -> Alphabetical -> Tube
Vary Radius 设置为By Vector，Radius按需缩放，0.0005比较合适，”Radius Factor”设为1
美化图片

图3 Paraview得到的力链图

参考文献

【1】徐小敏. 砂土液化及其判别的微观机理研究[D]. 浙江大学, 2012.

【2】https://www.zhihu.com/tardis/bd/art/273522056

离散元后处理与Python实践—利用Paraview绘制力链