首页/文章/ 详情

三角形法计算柔性双轴体变

lobby

2年前浏览1747

引言

实际中，三轴实验的体积变化是通过排水的体积来确定的。在柔性双轴的数值模拟中肯定是没有水的，所以我们没办法直接获得式样的体积变化，所以只能去直接计算式样的体积了。

这里采用在CAD等软件中普遍采用的三角形法来计算多边形的面积，下面介绍一下这个方法的原理。

一、三角形法原理

对于凸多边形，我们是可以将其直接分成多个三角形，分别计算三角形的面积，累加起来就是我们需要的面积了。如下图所示。

比较难理解的是式样形状不一定是凸多边形，有可能在局部是凹的形状，就比如下图这样的轮廓，这时候三角形分割的方法是否还有效呢？

答案是肯定的，下图中pt1-pt2-pt3围成的三角形是可以直接算的。之后的pt1-pt3-pt4-pt5围成的多边形面积，我们可以用pt1-pt4-pt5三角形面积S3减去pt1-pt3-pt4三角形面积S2。

这个道理就好比一个人先往后退了1m，然后前进了3m，问他产生的位移是多少道理一样。而决定这个面积是前进还是后退，取决于pt1-->ptn 向量到pt1-->ptn 1向量是顺时针还是逆时针转。

对于凸多边形，pt1-->ptn 向量到pt1-->ptn 1向量都是一个方向在转，所以累加即可。但是对于下图的轮廓，有时候顺时针，有时候逆时针的。我们就把顺时针产生的面积认为其在“前进”，逆时针产生的面积认为其在“后退”。

判断顺时针和逆时针的方法在高中应当讲过，可以利用叉乘，叉乘为负的时候是顺时针，叉乘为正的时候为逆时针。然后我们又想到叉乘的模也是有意义的，即为两个向量形成的平行四边形面积，除以2就是形成的三角形面积了。

二、代码实现

我们就在PFC里面来进行三角形法的实现。这里的pts为一组按顺序排列的点机，数据类型为array(vec2d)。pt0、pt1、pt2为每个三角形的三个点。vec1、vec2为其两个向量。这里假设vec1=(x1,y1,0) vec2=(x2,y2,0)，因为是二维的。他的叉乘应当是(0,0,x1*y2-x2*y1)，这里的叉乘向量的长度其实就体现在z位置的大小上，而顺时针还是逆时针体现在z位置的正负上。所以我们直接对所有三角形向量的叉乘的z分量进行累加即可。

这个函数是比较通用的，只要是多边形，处理成这样的数据传进去，都可以获得面积。

然后我们再看一下在柔性双轴中怎么去获取pts，其实比较简单的方法是通过id识别，因为膜颗粒的id号是有规律的。

这里为了通用性，用函数去识别这些颗粒。

首先先计算需要计算的膜颗粒数量，因为上下固定的膜颗粒是不用参加计算的。

第二步就是先左后右来获取膜颗粒的position。注意顺序必须是首尾相接的，对于foreach遍历，实际上会根据id大小先后大小顺序获取bp。左边膜颗粒我们需要从下往上，但是对于右边膜颗粒我们是需要从上往下的，所以对于右边的膜颗粒我们需要从数组的最后往前进行插入。

def PtsCount    pts_sum_count=0    loop foreach bp ball.groupmap("left_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            continue        endif        pts_sum_count =1    endloop        PtsCount=pts_sum_count*2end

def MoPts    moptstemp= array.create(PtsCount)    pt_count=1    loop foreach bp ball.groupmap("left_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            exit loop        endif        moptstemp(pt_count)=ball.pos(bp)        pt_count =1    endloop    pt_count=PtsCount    loop foreach bp ball.groupmap("right_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            exit loop        endif        moptstemp(pt_count)=ball.pos(bp)        pt_count-=1    endloop    MoPts=moptstempend

三、计算结果

这里给出完整的vol_tools文件，就是包括上面的三个函数：

def PtsCount    pts_sum_count=0    loop foreach bp ball.groupmap("left_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            continue        endif        pts_sum_count =1    endloop        PtsCount=pts_sum_count*2end

def MoPts    moptstemp= array.create(PtsCount)    pt_count=1    loop foreach bp ball.groupmap("left_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            exit loop        endif        moptstemp(pt_count)=ball.pos(bp)        pt_count =1    endloop    pt_count=PtsCount    loop foreach bp ball.groupmap("right_bianjie")        if ball.group(bp,10)="downMo" then            continue        endif        if ball.group(bp,10)="upMo" then            exit loop        endif        moptstemp(pt_count)=ball.pos(bp)        pt_count-=1    endloop    MoPts=moptstempend





def AreaPts(pts)    areaSum=0    pt0=pts(1)    pt_size=array.size(pts,1)    loop pt_count(2,pt_size-1)        pt1=pts(pt_count)        pt2=pts(pt_count 1)        vec1=pt1-pt0        vec2=pt2-pt0        areaSum =0.5*(comp.x(vec1)*comp.y(vec2)-comp.y(vec1)*comp.x(vec2))            endloop    AreaPts=(-1)*areaSumend