表面应变计算（Surface strain calculation）

1. 背景介绍

应力-应变曲线能够有效地帮助我们分析材料的力学性能，应力指物体单位面积上所受到力的大小，应变指物体的相对变形量。正应变描述的是沿着某一方向上长度的变化量，而切应变描述的是角度的变化量。
根据力的大小可根据胡克定律理论上计算得到任意位置应变。对于现代的测量方式而言，主要有引伸计、应变片、光学测量等方式，光学测量因非接触、精度高的特点，在工业测量领域受到高度关注。

2 变形梯度张量

光学测量使用数字图像相关法(Digital Image Correlation, DIC)获得被测物体表面的变形情况，因此，计算的是表面应变(Surface Strain)。在连续介质力学中，描述物体运动的方式有两种，即Lagrange描述与Euler描述。前者以不同时刻下物体的位置为研究对象，而后者以固定空间内的变化为对象。本文中我们在Lagrange描述下进行分析。

上图中，参考位置$\overrightarrow{PQ}$，在$t$刻为$\overrightarrow{pq}$，P和Q的位移量分别为$u(\mathbf{X})$,$u(\mathbf{X}+d\mathbf{X})$，该过程的运动描述为$\mathbf{x}=\chi (\mathbf{X},t)$。根据应变的定义可知，要想得到应变，就需要建立$d\mathbf{x}$与$d\mathbf{X}$之间的关系，根据图示，很容易能够得到$d\mathbf{x}=d\mathbf{X}+d\mathbf{u}$，进一步可知

其中，$d\mathbf{u}$是相对位移，$\mathbf{F}$为变形梯度张量(Deformation Gradient Tensor)。根据运动描述$\mathbf{x}=\chi (\mathbf{X},t)$的全微分形式

由此前两式可知，变形梯度张量$\mathbf{F}$实际上是该物体运动描述函数的梯度。$d\mathbf{x}$如果理解为变形后$\overrightarrow{pq}$的长度，那么变形前长度就是$d\mathbf{X}$，因此，只要能够得到$\mathbf{F}$，就能计算出应变大小。关于变形梯度张量$\mathbf{F}$的计算，若是已知运动方程，则只需要根据方程，求$\mathbf{F}=\left[\begin{array}{cc}\partial x/\partial X& \partial x/\partial Y\\ \partial y/\partial X& \partial y/\partial Y\end{array}\right]$的偏导数矩阵即可，如下所示的计算例子。
前述内容是关于变形梯度张量的代数意义，而它在几何上的意义可以用公式表示为$$\mathbf{F}=\mathbf{R}\mathbf{U}$$它表示的是连续介质刚体旋转与变形的过程，如下图所示，立方体在变形矩阵$\mathbf{U}$作用下变形为长方体，其后做刚体旋转$\mathbf{R}$，最终为$t$时刻状态。或者也可以先旋转再变形，$\mathbf{F}=\mathbf{V}\mathbf{R}$。

3. 表面应变计方法一

实际情况运动方程是未知的，那应该如何求解$\mathbf{F}$。根据前面的式子可知，$$\mathbf{F}=\mathbf{I}+{\nabla }_x\mathbf{u}$$如果能得到相对位移的梯度，那就能得到$\mathbf{F}$。使用DIC设备匹配散斑图像，重建能够得到匹配点的三维坐标，同时计算得到位移。但是，这些这些数据都是离散的，如何来求相对位移的梯度呢？常用的方法就是拟合，根据这些离散点的坐标为自变量，位移为因变量，拟合得到函数$u(X,Y)$,$v(X,Y)$。然后，求函数梯度，即可得到$\mathbf{F}$。根据$\mathbf{F}$，计算得到格林-拉格朗日应变张量(Green-Lagrange tensor)，

3.1 局部坐标系建立

DIC结果得到的实际上是一个没有厚度的面片，Surface strain或者Plain strain就是说，应变计算是一个沿着表面的、二维的度量。利用$3\times 3$，$7\times 7$或者其他局部范围内的点来构造一个平面，以此建立局部坐标系，有两种方式：一是这些点的拟合平面；二是在局部中心点处的切平面。
采用构造切平面的方式，使用局部$3\times 3$内点坐标，拟合二次曲面$z=f(x,y)$，法向量为

切向量有无数条，关于切向量的选取，只需要固定$x$或$y$为一常数，从而确定一个切向量。它的选取最好能与沿着物体变形的方向。

根据向量叉乘得到另一个坐标轴，由此，获得$\vec{n}$, $\vec{t}$, $\vec{b}$单位向量构建的局部坐标系

3.2 位移函数拟合

值得注意的是，每个点的位移量值不是将点投影到局部坐标系，然后做差。这种做法只适合于小弯曲变形情况，在Andresen 1999年撰写的论文中就有提到在大弯曲变形情况下，这种做法准确性很低。Surface strain, 是沿着表面方向发生的变形量，准确的来讲，应该使用弧长而不是位移向量，如图所示。换句话说，相当于把面展开成平面。，准确的方法应该是在上一步建立的局部坐标下，将重建的三维点，根据拟合的二次曲面$z=f(x,y)$沿展到$tob$面内，最后在$tob$面内得到每一个点的“位移量”。

假设点的位移模式为如下二次形式

现在已知每一个点在$tob$面内的坐标$(x_1,x_2)$，以及该点的在$x_1$方向的“位移”$u(x_1,x_2)$，如果选取的是$3\times 3$的点，9个方程，6个未知系数，利用最小二乘法求解，得到函数的表达式。

3.3 位移函数偏导

由前述得到位移函数的表达式之后，也就能计算得到$$\partial u/\partial x_1={\beta }_2$$$$\partial u/\partial x_2={\beta }_3$$同理得到另一个方向的偏导数，带入上述变形梯度张量$\mathbf{F}$的计算中，最终求解格林应变张量。根据3.1~3.3的步骤，遍历DIC重建得到的每一个点，由此可以获得物体表面变形的应变分布。

4. 表面应变计算方法二

通过第3节中的计算过程，其实可以看出来，连续变形过程的应变计算本质上讲就是位移的变化率，用公式来说就是$$strains=\nabla displacements$$这里通过搜集的资料，给大家分享第二种表面应变计算方法，其本质没有变，只是策略是利用了有限元理论中，三角形单元的应变计算。
既然DIC设备重建得到的是点的信息，那就可以将这些点连接为一个个的三角形，三角形内任意点的位移可以根据如下位移模式得到

方程中含有6个未知数，三角形的每一个定点坐标与位移量大小都是DIC获取的，每一点能够建立两个方程，那么3个点正好6个方程，可以求解$a_1,a_2,...,a_6$。同样地，就可以得到位移的梯度
$$\partial u/\partial x=a_2,\partial u/\partial y=a_3$$$$\partial v/\partial x=a_5,\partial v/\partial y=a_6$$代入变形梯度张量公式与格林应变公式，即可得到结果。三角形单元是一次形式的单元，单元内任意一点的应变是一个常数，那么在四个三角形中心点处的应变，需要求平均值。
同样地，遍历所有点，构造周围四个三角形，得到该点格林应变。由上述过程可以看出，该方法的前提就是获取点之间的间距必须足够小，否则对于大变形的情况精度会很差，而方法一的误差则主要来自于拟合误差。