CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    OF5版本以后重心坐标与绝对坐标的相互转换的问题。

    OpenFOAM
    3
    9
    1603
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 上级
      上级 最后由 上级 编辑

      最近在研读OpenFOAM中particle类的barycentric coordinates 和global coordinates转换矩阵的程序,遇到一点问题,现总结出来和大家讨论一下。
      OF5.0以后对颗粒追踪采用了barycentric方法,目的是更好的捕捉颗粒轨迹与tetrahedra网格面的交点,四面体内任意一点$\mathbf{p}$的global coordinates 可以表示为barycentric coordinates ($\lambda_i$)与顶点的乘积和:
      \begin{equation}
      \label{eq:barycentric coordinates}
      \mathbf{p}=\sum_{i=1}^{4}\lambda_i\mathbf{v_i}
      \end{equation}
      其中$\mathbf{v_i}$是四面体顶点坐标。
      示意图如下:
      1.JPG
      如上图,四面体中有一点$\mathbf{p}$,则$\mathbf{p}$点的重心坐标($\lambda_1,\lambda_2,\lambda_3,\lambda_4$)的计算方法为点$\mathbf{p}$将四面体分成的四个小四面体与大四面体的体积之比:
      \begin{equation}
      \label{eq:lambda1}
      \lambda_{1}=\frac{V_{\rm{pbcd}}}{V_{\rm{abcd}}}=1-\lambda_{2}-\lambda_{3}-\lambda_{4},
      \end{equation}
      \begin{equation}
      \label{eq:lambda2}
      \lambda_{2}=\frac{V_{\rm{pacd}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pacd}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)};,
      \end{equation}
      \begin{equation}
      \label{eq:lambda3}
      \lambda_{3}=\frac{V_{\rm{pabd}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pabd}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
      \end{equation}
      \begin{equation}
      \label{eq:lambda4}
      \lambda_{4}=\frac{V_{\rm{pabc}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pabc}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
      \end{equation}
      其中 $V_{\rm{pbcd}}$, $V_{\rm{pacd}}$, $V_{\rm{pabd}}$, 和 $V_{\rm{pabc}}$ 是小四面体的体积, $\mathbf{ap}=\mathbf{p}-\mathbf{a}$, $\mathbf{ab}=\mathbf{b}-\mathbf{a}$, $\mathbf{ac}=\mathbf{c}-\mathbf{a}$, $\mathbf{ad}=\mathbf{d}-\mathbf{a}$.
      定义从barycentric coordinates到global coordinates的转换矩阵$\mathbf{A}$,则
      \begin{equation}
      \begin{aligned}
      \label{eq:barycentric to global}
      \mathbf{p}=\mathbf{A}\boldsymbol{\lambda}\left(\mathbf{p}\right)\rightarrow
      \begin{bmatrix}x_{\rm{p}}\\y_{\rm{p}}\\z_{\rm{p}}\end{bmatrix}=
      \begin{bmatrix}x_{\rm{a}}&x_{\rm{b}}&x_{\rm{c}}&x_{\rm{d}}\\y_{\rm{a}}&y_{\rm{b}}&y_{\rm{c}}&y_{\rm{d}}\\z_{\rm{a}}&z_{\rm{b}}&z_{\rm{c}}&z_{\rm{d}}
      \end{bmatrix}
      \begin{bmatrix}\lambda_1\\
      \lambda_2\\
      \lambda_3\\
      \lambda_4\end{bmatrix};.
      \end{aligned}
      \end{equation}
      其中$\lambda_i$就是$\mathbf{p}$的重心坐标,矩阵$\mathbf{A}$就是由四面体各个顶点的global coordinates组成的坐标矩阵。
      既然有了正向变换,同样需要逆向变换,定义一个逆向变换矩阵$\mathbf{T}=\left[\mathbf{bd}\times\mathbf{bc}\quad\mathbf{ac}\times\mathbf{ad}\quad\mathbf{ad}\times\mathbf{ab}\quad\mathbf{ab}\times\mathbf{ac}\right]$将global coordinates转换为barycentric coordinates:
      \begin{equation}
      \label{eq:global to barycentric}
      \boldsymbol{\lambda}=\frac{\mathbf{ap}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}=\frac{\left[\mathbf{ap}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)\right]}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
      \end{equation}
      公式\eqref{eq:global to barycentric}的四个组成成分其实就是上面的公式\eqref{eq:lambda1},\eqref{eq:lambda2},\eqref{eq:lambda3},\eqref{eq:lambda4}。$\rm{det}\mathbf{A}$是矩阵$\mathbf{A}$的行列式,表示为$\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)$。
      至此,barycentric coordinates与global coordinates的正向变换和逆向变化都已经定义完成了,上面的分析是基于我看OF6的代码以及查阅相关的文献总结出来的。
      但是有一个问题困扰了我很久了:zoule: ,那就是逆向变换矩阵的第一个元素是:
      \begin{equation}
      \label{eq:the first element of T}
      \boldsymbol{\lambda}_1=\frac{\mathbf{ap}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
      \end{equation}
      这个和公式\eqref{eq:lambda1}的定义不同,因为分子上算的不是图中小四面体pbcd的体积,就这儿我无法理解为什么OF里要这么定义,并且我看OF5-OF8里面都是这么定义的,我觉得有点问题,所有发出来和大家讨论一下,我觉得我自己是解决不了这个问题了:141:

      1 条回复 最后回复 回复 引用
      • 李东岳
        李东岳 管理员 最后由 编辑

        简直是大作!!我一直在想研究下这方面工作没倒出空来。我先留言学习一下,在讨论

        感谢分享!!!!!!!:146: :146: :146:

        线上CFD课程开始报名:http://www.dyfluid.com/class.html

        CFD高性能服务器 http://dyfluid.com/servers.html

        1 条回复 最后回复 回复 引用
        • 李东岳
          李东岳 管理员 最后由 李东岳 编辑

          首先跟随你的思想:按照理论,应该是:
          \begin{equation}
          \boldsymbol{\lambda}=\frac{\mathbf{ap}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}=\frac{\left[\mathbf{bp}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)\right]}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
          \end{equation}
          同时
          \begin{equation}
          \mathbf{bp}\neq\mathbf{ap}
          \end{equation}
          所以代码跟理论对不上?

          线上CFD课程开始报名:http://www.dyfluid.com/class.html

          CFD高性能服务器 http://dyfluid.com/servers.html

          上级 1 条回复 最后回复 回复 引用
          • 上级
            上级 @李东岳 最后由 编辑

            @东岳 是的李老师,我的疑惑就在这儿,$\mathbf{bp}\neq\mathbf{ap}$,和代码对不上,细读代码发现不一样。

            1 条回复 最后回复 回复 引用
            • 李东岳
              李东岳 管理员 最后由 编辑

              问题解决了么?

              线上CFD课程开始报名:http://www.dyfluid.com/class.html

              CFD高性能服务器 http://dyfluid.com/servers.html

              1 条回复 最后回复 回复 引用
              • 上级
                上级 最后由 编辑

                仔细读了代码,现在还剩最后一点小问题,其实说到底还是原来的问题,barycentric追踪是通过bary displacement元素的正负来判断颗粒运动轨迹与四面体网格面的交点的,例如这个图形
                捕获.JPG
                颗粒$\mathbf{p}$运动到$\mathbf{n}$,轨迹$\mathbf{pn}$与面$\underline{\rm{bcd}}$相交于点$\mathbf{s}$,那么矢量$\mathbf{pn}$的barycentric displacement的计算方法为:
                \begin{equation}
                \begin{aligned}
                \label{eq:barycentric displacement}
                \boldsymbol{\lambda}_{\mathbf{p}\rightarrow\mathbf{n}}
                &=\boldsymbol{\lambda}\left(\mathbf{n}\right)-\boldsymbol{\lambda}\left(\mathbf{p}\right)=\frac{\mathbf{an}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}-\frac{\mathbf{ap}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}=\frac{\mathbf{pn}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}\
                &=\frac{\left[\mathbf{pn}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)\right]}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)};.
                \end{aligned}
                \end{equation}
                此时barycentric dispalcement的符号为(-,+++),可以判断轨迹与面$\underline{\rm{bcd}}$相交,同根根据$\mathbf{p}$和$\mathbf{n}$的barycentric coordinates可以求出点$\mathbf{s}$的barycentric coordinates,再转换为点$\mathbf{s}的$gobal coordinates,然后进如下一个四面体网格追踪。
                上面我只是简述了颗粒在一个四面体网格内的追踪过程,从计算barycentric dispalcement的角度来看,逆向变换矩阵$\mathbf{T}$是没有问题的,但是问题还是在于从barycentric coordinates to global coordinates的过程,矩阵$\mathbf{T}$的第一个元素有问题。如果把\ref{eq: the first element of T}变为$\boldsymbol\lambda_1=1-\boldsymbol\lambda_2-\boldsymbol\lambda_3-\boldsymbol\lambda_4$就好了,但是程序里没有这么写,说到底还是原来的问题:jingya: 如果李老师有时间的话可不可以读一下particle类,我觉得您可以帮我释疑:xinxin:

                上级 1 条回复 最后回复 回复 引用
                • 上级
                  上级 @上级 最后由 编辑

                  上面那个公式引用是公式(8)

                  1 条回复 最后回复 回复 引用
                  • B
                    BlookCFD 最后由 编辑

                    (7)这个公式确实是错的,但是particle里面都用的是(11)这个公式,用来计算“位移”的barycentric坐标,而不是点的barycentric坐标,OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H有计算点的barycentric坐标,你可以比较下。

                    上级 1 条回复 最后回复 回复 引用
                    • 上级
                      上级 @BlookCFD 最后由 编辑

                      @BlookCFD 谢谢您的指点,我有看了一眼程序,我理解错了,OF里面是通过计算barycentric displacement来获得barycentric coordinates的,我一直被src/lagrangian/basic/particle/particle.C里面的一句代码给误导了,现贴出来:

                      particle.C 1062-1096 OpenFOAM6
                      void Foam::particle::correctAfterInteractionListReferral(const label celli)
                       {
                           // Get the position from the barycentric data
                           const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d());
                       
                           // Create some arbitrary topology for the supplied cell
                           celli_ = celli;
                           tetFacei_ = mesh_.cells()[celli_][0];
                           tetPti_ = 1;
                           facei_ = -1;
                       
                           // Get the reverse transform and directly set the coordinates from the
                           // position. This isn't likely to be correct; the particle is probably not
                           // in this tet. It will, however, generate the correct vector when the
                           // position method is called. A referred particle should never be tracked,
                           // so this approximate topology is good enough. By using the nearby cell we
                           // minimize the error associated with the incorrect topology.
                           coordinates_ = barycentric(1, 0, 0, 0);
                           if (mesh_.moving())
                           {
                               Pair<vector> centre;
                               FixedList<scalar, 4> detA;
                               FixedList<barycentricTensor, 3> T;
                               movingTetReverseTransform(0, centre, detA, T);
                               coordinates_ += (pos - centre[0]) & T[0]/detA[0];
                           }
                           else
                           {
                               vector centre;
                               scalar detA;
                               barycentricTensor T;
                               stationaryTetReverseTransform(centre, detA, T);
                               coordinates_ += (pos - centre) & T/detA;
                           }
                       }
                      

                      其中的这句代码

                      coordinates_ += (pos - centre) & T/detA;
                      

                      其实前面已经先定义了

                      coordinates_ = barycentric(1, 0, 0, 0);
                      

                      后面计算单点的barycentric coordinates其实还是用的barycentric displacement的概念来转换的。

                      困扰了我两周的问题终于解决了,现在异常开心,在此对@东岳 和@BlookCFD 表示万分感谢。

                      1 条回复 最后回复 回复 引用
                      • First post
                        Last post