一种基于平面傅立叶轮廓分析的三维颗粒构形方法
技术领域
本发明涉及一种基于平面傅立叶轮廓分析的三维颗粒构形方法,属于岩土工程技术领域。
背景技术
在岩土介质中,宏观的变形破坏规律及力学特性(如破坏模式、裂纹扩展、承载能力等)很大程度上依赖于其内部细观结构特征(如粒度组成、颗粒排列方式等),而天然岩土介质多由尺度不等的颗粒混合构成,如粗粒土、土石混合介质等,颗粒在介质中起骨架支撑作用,粒径、级配、外轮廓、粗糙度等细观特征对颗粒的物理力学特性(如剪切强度、膨胀性、抗压性等)有重要影响。
目前,工程试验中常将碎块石简化为规则几何体,用圆形、三角形等简单凸边形进行考虑,但是这种简单的细观特征分析方法与实际情况相差较远,无法准确反映介质的破坏机制;数字图像处理技术可较为准确地反映原位颗粒细观特征,但数字图像处理技术易受外界各种因素影响,具有局限性;数值模拟方法(如离散单元法),虽然能很好地模拟介质的开裂、大变形,在岩土工程、散体力学等领域颇受关注,但是颗粒模型构建的准确性与合理性却一直制约着该方法的应用。
发明内容
本发明所要解决的技术问题是提供一种基于平面傅立叶轮廓分析的三维颗粒构形方法,能实现基于个别典型颗粒,随机生成表面接近任意不规则岩土颗粒的重构,可较好地逼近真实岩土颗粒。
本发明为了解决上述技术问题采用以下技术方案:本发明设计了一种基于平面傅立叶轮廓分析的三维颗粒构形方法,用于实现对目标颗粒的三维构形,包括如下步骤:
步骤A.以目标颗粒的形心作为原点,构建正交三维坐标系,首先以XOY平面切割目标颗粒,获得一个过形心的水平切面,作为目标颗粒主切面,然后过形心绕z轴相隔预设角度切割获得各个垂直切面,得到一系列垂直轮廓线,针对目标颗粒主切面轮廓、各个垂直切面轮廓,分别进行二维傅里叶分析,获取目标颗粒每一个切面轮廓对应预设阶数目的各阶傅里叶谱值,然后进入步骤B;
步骤B.根据所获目标颗粒切面的各阶傅里叶谱值,针对目标颗粒切面的各阶谱,分别赋予随机幅角,且各随机幅角相位服从[0,2π]均匀分布;接着利用傅里叶逆变换,分别获得预设各阶数目的过形心切面轮廓,作为预设阶数目的目标颗粒随机切面轮廓,然后进入步骤C;
步骤C.针对目标颗粒各切面轮廓,以及预设阶数目的目标颗粒随机切面轮廓进行调整,使得各个切面轮廓相互匹配,获得目标颗粒三维轮廓,然后进入步骤D;
步骤D.基于标准球面三角形网格,根据所获目标颗粒三维轮廓,推断获得目标颗粒的径向距离,进而获得目标颗粒的三维构形。
作为本发明的一种优选技术方案,所述步骤A包括如下步骤:
步骤A1.以目标颗粒的形心作为原点,构建三维坐标系,首先获取目标颗粒过形心的XOY切面,作为目标颗粒主切面,过形心绕z轴相隔预设角度进行竖向切面切割,获得M个竖向切面轮廓,即获得目标颗粒M+1个切面轮廓;遍历M+1个切面轮廓线,任取一个令该轮廓线在水平面内y为正方向的交线,作为X'正方向,则形成了X'OZ轮廓面,过切面中心预设角度在轮廓线上划分为预设N1个轮廓点,同时设置傅里叶变换的阶数为N,则每一个轮廓点均可以用极坐标(r,θ)表示,它可以用如下傅里叶系数来表示:
分别获得各个轮廓点的径向距离r(θi),i=1~N1,然后进入步骤A2;其中,θi表示所取颗粒切面轮廓上第i个轮廓点与平面内水平X'轴的夹角,逆时针为正;r0表示目标颗粒切面轮廓上轮廓点到颗粒质心的平均距离,An和Bn为目标颗粒切面轮廓线所对应各阶中第n阶的傅里叶系数,0≤n≤N;
步骤A2.从M+1个切面轮廓线中任取其中一个的二维傅里叶谱{An,Bn}n=1~N,将其矢量取模后除以平均距离r0得不考虑幅角的傅里叶谱:
获得目标颗粒各切面轮廓对应预设阶数目的各阶傅里叶谱值Dn。计算{An,Bn,Dn}n=1~N的过程称为傅里叶变换,利用这些计算计算(r,θ)过程称为傅里叶逆变换。
作为本发明的一种优选技术方案,所述步骤B包括如下步骤:
步骤B1.根据所获目标颗粒切面对应预设阶数目的各阶傅里叶谱值Dn,针对目标颗粒切面的各阶傅里叶谱,分别赋予随机幅角ηn,且各随机幅角ηn相位服从[0,2π]均匀分布,则按如下公式:
An=Dn·cosηn
Bn=Dn·sinηn
分别获得目标颗粒切面所对应各阶的傅里叶系数An和Bn,然后进入步骤B2;其中,ηn表示针对目标颗粒切面所对应的第n阶赋予的随机幅角;
步骤B2.根据目标颗粒主切面所对应各阶的傅里叶系数An和Bn,利用傅里叶逆变换,分别获得预设各阶数目的过形心切面轮廓,作为预设阶数目的目标颗粒随机切面轮廓。
作为本发明的一种优选技术方案,所述步骤C中,针对目标颗粒各切面轮廓,以及预设阶数目的目标颗粒随机切面轮廓进行调整,即针对个轮廓,执行如下步骤C1至步骤C2,使得各个切面轮廓相互匹配,更新作为各个目标颗粒待处理切面轮廓;
步骤C1.以水平XOY面上的轮廓线作为基准,对各竖向轮廓线进行调整,具体为:根据各纵向剖面轮廓线分别与Z轴的正向交点(0,z1)...(0,zi)...(0,zN1),以及分别与Z轴的负向交点分别(0,z'1)...(0,zi')...(0,z'N1),按如下公式:
求出正向纵坐标的平均数负向纵坐标的平均数并分别作为轮廓线Z轴正负方向顶点的基准点;
步骤C2.如图2(a)所示任一竖向轮廓线,引入校正因子δC,调整各纵向剖面切线Z轴的位置,采用下列公式对各轮廓的每个半径,根据局部角θC进行校正,θC表示轮廓点-原点连线与Z轴正向的夹角;
δ1,δ2,δ3和δ4分别是点Ci,B1,C'i和B2处的局部校正因子:
进而完成对各个轮廓的调整,使得各个切面轮廓相互匹配,获得目标颗粒三维轮廓;其中Ci,C'i为所选择竖向剖面位于Z轴正负方向的点,B1,B2分别为所选择竖向剖面与水平切面相交方向的轮廓点,A1,A2分别为所选择竖向剖面与水平切面轮廓线的交点。
作为本发明的一种优选技术方案,所述步骤D中,根据所获目标颗粒三维轮廓,利用径向函数插值法,获得目标颗粒的径向距离,进而获得目标颗粒的三维构形。
作为本发明的一种优选技术方案,所述步骤D包括如下步骤:
步骤D1.将目标颗粒三维轮廓形心移至坐标轴原点,并采用球坐标系对目标颗粒三维轮廓进行描述,球坐标上的各点可通过下式转化成笛卡尔坐标:
zij=rijsin(θi)
其中,rij表示目标颗粒三维轮廓形心到目标颗粒表面点的极半径,0≤θ≤2π,
步骤D2.将目标颗粒三维轮廓看作是由多个三角形面构成的,球面网格点服从均匀随机分布,采用Delaunay化,将目标颗粒三维轮廓化为闭合三角形构成的空间球面;
步骤D3.分别利用XOY面、XOZ面、YOZ面三个正交面,即竖向轮廓切面2个,总轮廓平面数目3,绘制目标颗粒在三个正交平面上的外轮廓,其中,若每个轮廓线所处的X'OZ平面内轮廓由极坐标r(θ),θ为轮廓点与原点连线与X轴夹角,逆时针为正,定义的各个点连接而成,O为颗粒形心所构的原点,A1,A2,A3,A4分别为XOY平面轮廓线与X轴正方向、Y轴正方向、X轴负方向、Y轴负方向的交点;B1,B2,B3,B4分别为XOZ平面轮廓线与X轴正方向、Z轴正方向、X轴负方向、Z轴负方向的交点;C1,C2,C3,C4分别为YOZ平面轮廓线与Y轴正方向、Z轴正方向、Y轴负方向、Z轴负方向的交点;则三个剖面轮廓线形状沿着他们的公共轴线满足以下条件:
OA1=OB1和OA3=OB3(x轴)
OA2=OC1和OA4=OC3(y轴)
OB2=OC2和OB4=OC4(z轴)
即三个轮廓切面满足几何匹配条件;利用目标颗粒的三个剖面外轮廓,以及步骤D2中标准球面网格点,获得目标颗粒的径向距离,构造目标颗粒的表面形状;
步骤D4.基于多个过原点的2D轮廓是相容的条件下,针对目标颗粒各切面轮廓,以及预设数目的目标颗粒随机切面轮廓,即采用M+1个轮廓面来表征目标颗粒的表面形状,各个竖向轮廓面会将水平XOY面切割为2M个象限,任取第Si、Si+1剖面象限,i为任意一个竖向切面编号,如i为M,则i+1取值为1,在两个切面包围的颗粒表面任意选取一点,向XOY面投影得到P0,即任选过P0点的Si剖面,选取剖切线上靠近赤道的任意一点P,P到原点距离为γ1,P与原点的连线与Si面的夹角为δ;原颗粒表面在XOY面上的点为P0,与原点的距离为γ0;在赤道附近的δ角范围内,选取任意一点Pi,与原点的距离为γi,Pi与原点连线与Si面的夹角为δ1,则P'i到原点的插值距离γ'i为:
遍历所有基准球上的轮廓点,不改变其θ,但是球心到表面的距离采用插值距离替换,即令rij取值为γ'i,最后获得目标颗粒的随机三维构形。
本发明所述一种基于平面傅立叶轮廓分析的三维颗粒构形方法采用以上技术方案与现有技术相比,具有以下技术效果:本发明设计的基于平面傅立叶轮廓分析的三维颗粒构形方法,通过单个颗粒的多切面轮廓分析与重构,建立了一种简易三维颗粒构造方法,能实现基于个别典型颗粒,随机生成表面接近任意不规则岩土颗粒的重构,可较好地逼近真实岩土颗粒,有效解决颗粒离散元数值模拟领域任意凹凸三维颗粒的随机构造问题,使得离散元计算更加逼近实际工程实践。
附图说明
图1a为傅里叶谱随机生成的多剖切面的平面图;
图1b为傅里叶谱随机生成的多剖切面的三维图;
图2a和图2b为任意纵向剖切面与Z轴交点图;
图3a为球坐标系中轮廓点表征方法;
图3b为球坐标构造的球面表征方法;
图4a为248个三角形网格化的单元球面示意;
图4b为4032个三角形网格化的单元球面示意;
图4c为25368个三角形网格化的单元球面示意;
图5a为颗粒原始示意图;
图5b为颗粒轮基于XOY面的示意图;
图5c为颗粒轮基于XOZ面的示意图;
图5d为颗粒轮基于YOZ面的示意图;
图5e为颗粒轮的3D剖切示意图;
图5f为颗粒轮廓线基于XOY面的示意图;
图5g为颗粒轮廓线基于XOZ面的示意图;
图5h为颗粒轮廓线基于YOZ面的示意图;
图6a为基于纵向点调整的颗粒径向距离推断示意图;
图6b为基于赤道附近点调整的颗粒径向距离推断示意图;
图6c为颗粒径向距离推断的多剖切面俯视示意图;
图7a为示例原颗粒视图,图7e为图7a示例的推断颗粒视图;
图7b为示例原颗粒视图,图7f为图7b示例的推断颗粒视图;
图7c为示例原颗粒视图,图7g为图7c示例的推断颗粒视图;
图7d为示例原颗粒视图,图7h为图7d示例的推断颗粒视图。
具体实施方式
下面结合说明书附图对本发明的具体实施方式作进一步详细的说明。
本发明设计了一种基于平面傅立叶轮廓分析的三维颗粒构形方法,用于实现对目标颗粒的三维构形,实际应用当中,具体包括如下步骤:
步骤A.如图1a和图1b所示,以目标颗粒的形心作为原点,构建正交三维坐标系,首先以XOY平面切割目标颗粒,获得一个过形心的水平切面,作为目标颗粒主切面,然后过形心绕z轴相隔预设角度切割获得各个垂直切面,得到一系列垂直轮廓线,针对目标颗粒主切面轮廓、各个垂直切面轮廓,分别进行二维傅里叶分析,获取目标颗粒每一个切面轮廓对应预设阶数目的各阶傅里叶谱值,然后进入步骤B。
上述步骤A具体包括如下步骤A1至步骤A2:
步骤A1.以目标颗粒的形心作为原点,构建三维坐标系(原点为o,三个坐标轴分别为X、Y、Z,其中Z轴数值向上,XOY面水平),首先获取目标颗粒过形心的XOY切面,作为目标颗粒主切面,过形心绕z轴相隔预设角度进行竖向切面切割,获得M个竖向切面轮廓,即获得目标颗粒M+1个切面轮廓;遍历M+1个切面轮廓线,任取一个令该轮廓线在水平面内y为正方向的交线,作为X'正方向,则形成了X'OZ轮廓面,过切面中心预设角度在轮廓线上划分为预设N1个轮廓点(实际应用中取为50),同时设置傅里叶变换的阶数为N(N应该取2的幂方,取为128),则每一个轮廓点均可以用极坐标(r,θ)表示,它可以用如下傅里叶系数来表示:
分别获得各个轮廓点的径向距离r(θi),i=1~N1,然后进入步骤A2;其中,θi表示所取颗粒切面轮廓上第i个轮廓点与平面内水平X'轴的夹角,逆时针为正;r0表示目标颗粒切面轮廓上轮廓点到颗粒质心的平均距离,An和Bn为目标颗粒切面轮廓线所对应各阶中第n阶的傅里叶系数,0≤n≤N。
步骤A2.从M+1个切面轮廓线中任取其中一个的二维傅里叶谱{An,Bn}n=1~N,将其矢量取模后除以平均距离r0得不考虑幅角的傅里叶谱:
获得目标颗粒各切面轮廓对应预设阶数目的各阶傅里叶谱值Dn。计算{An,Bn,Dn}n=1~N的过程称为傅里叶变换,利用这些计算计算(r,θ)过程称为傅里叶逆变换。
步骤B.根据所获目标颗粒切面的各阶傅里叶谱值,针对目标颗粒切面的各阶谱,分别赋予随机幅角,且各随机幅角相位服从[0,2π]均匀分布;接着利用傅里叶逆变换,分别获得预设各阶数目的过形心切面轮廓,作为预设阶数目的目标颗粒随机切面轮廓,如图1b所示,然后进入步骤C。
上述步骤B具体包括如下步骤B1至步骤B2:
步骤B1.根据所获目标颗粒切面对应预设阶数目的各阶傅里叶谱值Dn,针对目标颗粒切面的各阶傅里叶谱,分别赋予随机幅角ηn,且各随机幅角ηn相位服从[0,2π]均匀分布,则按如下公式:
An=Dn·cosηn
Bn=Dn·sinηn
分别获得目标颗粒切面所对应各阶的傅里叶系数An和Bn,然后进入步骤B2;其中,ηn表示针对目标颗粒切面所对应的第n阶赋予的随机幅角。
步骤B2.根据目标颗粒主切面所对应各阶的傅里叶系数An和Bn,利用傅里叶逆变换,分别获得预设各阶数目的过形心切面轮廓,作为预设阶数目的目标颗粒随机切面轮廓。
步骤C.针对目标颗粒各切面轮廓,以及预设阶数目的目标颗粒随机切面轮廓进行调整,使得各个切面轮廓相互匹配,如图2a至图2b所示,获得目标颗粒三维轮廓,然后进入步骤D。
上述步骤C中,针对目标颗粒各切面轮廓,以及预设阶数目的目标颗粒随机切面轮廓进行调整,即针对个轮廓,执行如下步骤C1至步骤C2,使得各个切面轮廓相互匹配,更新作为各个目标颗粒待处理切面轮廓。
步骤C1.以水平XOY面上的轮廓线作为基准,对各竖向轮廓线进行调整,具体为:根据各纵向剖面轮廓线分别与Z轴的正向交点(0,z1)...(0,zi)...(0,zN1),以及分别与Z轴的负向交点分别(0,z'1)...(0,zi')...(0,z'N1),按如下公式:
求出正向纵坐标的平均数负向纵坐标的平均数并分别作为轮廓线Z轴正负方向顶点的基准点。
步骤C2.如图2(a)所示任一竖向轮廓线,引入校正因子δC,调整各纵向剖面切线Z轴的位置,采用下列公式对各轮廓的每个半径,根据局部角θC进行校正,θC表示轮廓点-原点连线与Z轴正向的夹角;
δ1,δ2,δ3和δ4分别是如图2b上点Ci,B1,C'i和B2处的局部校正因子:
进而完成对各个轮廓的调整,使得各个切面轮廓相互匹配,获得目标颗粒三维轮廓;其中Ci,C'i为所选择竖向剖面位于Z轴正负方向的点,B1,B2分别为所选择竖向剖面与水平切面相交方向的轮廓点,A1,A2分别为所选择竖向剖面与水平切面轮廓线的交点。
步骤D.如图3a、图3b,以及图4a至图4c所示,基于标准球面三角形网格,根据所获目标颗粒三维轮廓,利用径向函数插值法,推断获得目标颗粒的径向距离,进而获得目标颗粒的三维构形,如图5a至图5h、图6a至图6c,以及图7a至图7h所示。
上述步骤D实际应用中,具体包括如下步骤:
步骤D1.将目标颗粒三维轮廓形心移至坐标轴原点,并采用球坐标系对目标颗粒三维轮廓进行描述,如图3a所示,球坐标上的各点可通过下式转化成笛卡尔坐标:
zij=rijsin(θi)
其中,rij表示目标颗粒三维轮廓形心到目标颗粒表面点的极半径,0≤θ≤2π,
步骤D2.如图3b所示,将目标颗粒三维轮廓看作是由多个三角形面构成的,球面网格点服从均匀随机分布,采用Delaunay化,将目标颗粒三维轮廓化为闭合三角形构成的空间球面。
步骤D3.分别利用XOY面、XOZ面、YOZ面三个正交面,即竖向轮廓切面2个,总轮廓平面数目3,绘制目标颗粒在三个正交平面上的外轮廓,其中,若每个轮廓线所处的X'OZ平面内轮廓由极坐标r(θ),θ为轮廓点与原点连线与X轴夹角,逆时针为正,定义的各个点连接而成,如图5a至图5h所示,O为颗粒形心所构的原点,A1,A2,A3,A4分别为XOY平面轮廓线与X轴正方向、Y轴正方向、X轴负方向、Y轴负方向的交点;B1,B2,B3,B4分别为XOZ平面轮廓线与X轴正方向、Z轴正方向、X轴负方向、Z轴负方向的交点;C1,C2,C3,C4分别为YOZ平面轮廓线与Y轴正方向、Z轴正方向、Y轴负方向、Z轴负方向的交点;则三个剖面轮廓线形状沿着他们的公共轴线满足以下条件:
OA1=OB1andOA3=OB3(x轴)
OA2=OC1andOA4=OC3(y轴)
OB2=OC2andOB4=OC4(z轴)
即三个轮廓切面满足几何匹配条件;利用目标颗粒的三个剖面外轮廓,以及步骤D2中标准球面网格点,获得目标颗粒的径向距离,构造目标颗粒的表面形状。
步骤D4.如图6a至图6c所示,基于多个过原点的2D轮廓是相容的条件下,针对目标颗粒各切面轮廓,以及预设数目的目标颗粒随机切面轮廓,即采用M+1个轮廓面来表征目标颗粒的表面形状,各个竖向轮廓面会将水平XOY面切割为2M个象限,任取第Si、Si+1剖面象限,i为任意一个竖向切面编号,如i为M,则i+1取值为1,在两个切面包围的颗粒表面任意选取一点,向XOY面投影得到P0,即任选过P0点的Si剖面,选取剖切线上靠近赤道的任意一点P,P到原点距离为γ1,P与原点的连线与Si面的夹角为δ;原颗粒表面在XOY面上的点为P0,与原点的距离为γ0;在赤道附近的δ角范围内,选取任意一点Pi,与原点的距离为γi,Pi与原点连线与Si面的夹角为δ1,则P'i到原点的插值距离γ'i为:
如图4a至图4c所示,遍历所有基准球上的轮廓点,不改变其θ,但是球心到表面的距离采用插值距离替换,即令rij取值为γ'i,最后获得目标颗粒的随机三维构形,如图7e至图7h所示。
将上述所设计基于平面傅立叶轮廓分析的三维颗粒构形方法,应用于实际当中,利用三个正交剖面,用38052个三角形构成的标准球面,推断出颗粒轮廓如图7e至图7h所示,这些颗粒基本参数如下:面积2.9306(原颗粒3.2049),体积0.4614(原颗粒0.5236),球度0.9445(原颗粒0.9802),径向形状参数为0.8559(原颗粒0.8406),二者相似度很高,表明采用少数几个切面轮廓线反推原颗粒的轮廓是可行的。本发明专利中,颗粒形状构造受多切面轮廓线的傅里叶谱控制,因此利用不同切面组合的傅里叶谱即可构造出大量服从同一统计特征的随机颗粒轮廓。
上面结合附图对本发明的实施方式作了详细说明,但是本发明并不限于上述实施方式,在本领域普通技术人员所具备的知识范围内,还可以在不脱离本发明宗旨的前提下做出各种变化。