计算科学 - 在 FEA 中计算 8 个节点砖单元的向外法线和表面积 - 吾爱随笔录

在 FEA 中计算 8 个节点砖单元的向外法线和表面积

计算科学有限元计算几何正则

2021-12-27 11:09:27

我有一个立方体，它通过平分每个边被分成 8 个小立方体，我试图找出每个面的表面积和它们对应的外法线。此操作是在有限元网格上完成的，因此我使用 shape(basis) 函数将立方体转换为等参数形式，然后尝试提取面积和法线。

这是代码的一部分：

 program polyhedron 


DO INPT=1,8  !LOOP OVER 8 SMALL CUBES

! GAUSS POINTS
XII(1,1) = MONE/THREE**HALF
XII(1,2) = MONE/THREE**HALF
XII(1,3) = MONE/THREE**HALF
XII(2,1) = ONE/THREE**HALF
XII(2,2) = MONE/THREE**HALF
XII(2,3) = MONE/THREE**HALF
XII(3,1) = ONE/THREE**HALF
XII(3,2) = ONE/THREE**HALF
XII(3,3) = MONE/THREE**HALF
XII(4,1) = MONE/THREE**HALF
XII(4,2) = ONE/THREE**HALF
XII(4,3) = MONE/THREE**HALF
XII(5,1) = MONE/THREE**HALF
XII(5,2) = MONE/THREE**HALF
XII(5,3) = ONE/THREE**HALF
XII(6,1) = ONE/THREE**HALF
XII(6,2) = MONE/THREE**HALF
XII(6,3) = ONE/THREE**HALF
XII(7,1) = ONE/THREE**HALF
XII(7,2) = ONE/THREE**HALF
XII(7,3) = ONE/THREE**HALF
XII(8,1) = MONE/THREE**HALF
XII(8,2) = ONE/THREE**HALF
XII(8,3) = ONE/THREE**HALF

        XI(1) = XII(INPT,1)
        XI(2) = XII(INPT,2)
        XI(3) = XII(INPT,3) 


DO I=1,8
 DO J=1,3
    dNdxi(I,J) =  ZERO
 END DO
END DO

!HEXAHEDRAL SHAPE FUNCTION DERIVATIVES
dNdxi(1,1) =  MONE/EIGHT*(ONE-XI(2))*(ONE-XI(3))
dNdxi(1,2) =  MONE/EIGHT*(ONE-XI(1))*(ONE-XI(3))
dNdxi(1,3) =  MONE/EIGHT*(ONE-XI(1))*(ONE-XI(2))
dNdxi(2,1) =  ONE/EIGHT*(ONE-XI(2))*(ONE-XI(3))
dNdxi(2,2) =  MONE/EIGHT*(ONE+XI(1))*(ONE-XI(3))
dNdxi(2,3) =  MONE/EIGHT*(ONE+XI(1))*(ONE-XI(2))
dNdxi(3,1) =  ONE/EIGHT*(ONE+XI(2))*(ONE-XI(3))
dNdxi(3,2) =  ONE/EIGHT*(ONE+XI(1))*(ONE-XI(3))
dNdxi(3,3) =  MONE/EIGHT*(ONE+XI(1))*(ONE+XI(2))
dNdxi(4,1) =  MONE/EIGHT*(ONE+XI(2))*(ONE-XI(3))
dNdxi(4,2) =  ONE/EIGHT*(ONE-XI(1))*(ONE-XI(3))
dNdxi(4,3) =  MONE/EIGHT*(ONE-XI(1))*(ONE+XI(2))
dNdxi(5,1) =  MONE/EIGHT*(ONE-XI(2))*(ONE+XI(3))
dNdxi(5,2) =  MONE/EIGHT*(ONE-XI(1))*(ONE+XI(3))
dNdxi(5,3) =  ONE/EIGHT*(ONE-XI(1))*(ONE-XI(2))
dNdxi(6,1) =  ONE/EIGHT*(ONE-XI(2))*(ONE+XI(3))
dNdxi(6,2) =  MONE/EIGHT*(ONE+XI(1))*(ONE+XI(3))
dNdxi(6,3) =  ONE/EIGHT*(ONE+XI(1))*(ONE-XI(2))
dNdxi(7,1) =  ONE/EIGHT*(ONE+XI(2))*(ONE+XI(3))
dNdxi(7,2) =  ONE/EIGHT*(ONE+XI(1))*(ONE+XI(3))
dNdxi(7,3) =  ONE/EIGHT*(ONE+XI(1))*(ONE+XI(2))
dNdxi(8,1) =  MONE/EIGHT*(ONE+XI(2))*(ONE+XI(3))
dNdxi(8,2) =  ONE/EIGHT*(ONE-XI(1))*(ONE+XI(3))
dNdxi(8,3) =  ONE/EIGHT*(ONE-XI(1))*(ONE+XI(2))



      dXdXi = zero
      dXdEta = zero
      dXdZeta = zero
      dYdXi = zero
      dYdEta = zero
      dYdZeta = zero
      dZdXi = zero
      dZdEta = zero
      dZdZeta = zero

      do k=1,8
         dXdXi = dXdXi + dNdxi(k,1)*coords(1,k)
         dXdEta = dXdEta + dNdxi(k,2)*coords(1,k)
         dXdZeta = dXdZeta + dNdxi(k,3)*coords(1,k)
         dYdXi = dYdXi + dNdxi(k,1)*coords(2,k)
         dYdEta = dYdEta + dNdxi(k,2)*coords(2,k)
         dYdZeta = dYdZeta + dNdxi(k,3)*coords(2,k)
         dZdXi = dZdXi + dNdxi(k,1)*coords(3,k)
         dZdEta = dZdEta + dNdxi(k,2)*coords(3,k)
         dZdZeta = dZdZeta + dNdxi(k,3)*coords(3,k)
      enddo

      DO face=1,6  !LOOP OVER ALL 6 FACES OF THE SMALL CUBE         

      ! Jacobian of the mapping
      !
      if((face.eq.1).or.(face.eq.2)) then
         ! zeta = constant on this face
         dA = dsqrt((dYdXi*dZdEta - dYdEta*dZdXi)**two+ (dXdXi*dZdEta -&
         &dXdEta*dZdXi)**two+ (dXdXi*dYdEta - dXdEta*dYdXi)**two)
      elseif((face.eq.3).or.(face.eq.5)) then
         ! eta = constant on this face
         dA = dsqrt((dYdXi*dZdZeta - dYdZeta*dZdXi)**two+ (dXdXi*dZdZeta -&
         &dXdZeta*dZdXi)**two+ (dXdXi*dYdZeta - dXdZeta*dYdXi)**two)
      elseif((face.eq.4).or.(face.eq.6)) then
         ! xi = constant on this face
         dA = dsqrt((dYdEta*dZdZeta - dYdZeta*dZdEta)**two+ (dXdEta*dZdZeta -&
         dXdZeta*dZdEta)**two+ (dXdEta*dYdZeta - dXdZeta*dYdEta)**two)

      endif

  write(*,*) 'face',face,'area is dA', dA

        if((face.eq.1).or.(face.eq.2)) then
         ! zeta = constant on this face
         normal(1,1) = dYdXi*dZdEta - dYdEta*dZdXi
         normal(2,1) = dXdXi*dZdEta - dXdEta*dZdXi
         normal(3,1) = dXdXi*dYdEta - dXdEta*dYdXi
         if(face.eq.1) normal = -normal
      elseif((face.eq.3).or.(face.eq.4)) then
         ! eta = constant on this face
         normal(1,1) = dYdXi*dZdZeta - dYdZeta*dZdXi
         normal(2,1) = dXdXi*dZdZeta - dXdZeta*dZdXi
         normal(3,1) = dXdXi*dYdZeta - dXdZeta*dYdXi
         if(face.eq.3) normal = -normal
      elseif((face.eq.5).or.(face.eq.6)) then
         ! xi = constant on this face
         normal(1,1) = dYdEta*dZdZeta - dYdZeta*dZdEta
         normal(2,1) = dXdEta*dZdZeta - dXdZeta*dZdEta
         normal(3,1) = dXdEta*dYdZeta - dXdZeta*dYdEta
         if(face.eq.6) normal = -normal

      endif
      mag = dsqrt(normal(1,1)**two+normal(2,1)**two+normal(3,1)**two)

      normal(1,1) = normal(1,1)/mag
      normal(2,1) = normal(2,1)/mag
      normal(3,1) = normal(3,1)/mag

      write(*,*) 'face', face, 'normal is ',normal

     END DO

每个小立方体的每个表面的面积应该是 0.25，在这种情况下，但这个代码没有返回。法线也没有正确定向，这里立方体与全局轴对齐，所以它应该是 1 或 -1。

注：一。这可能永远不是立方体，它可以是任何不规则的六面体，因此面积不会总是相等，因此我们需要计算它们中的每一个。湾。面可能沿着不同的方向定向，因此等参数变换是必要的。

此外，我已经使用每个面的向量来解决这个问题，而不使用等参数变换，我使用立方体的面的对角线并计算它们的叉积来找到面积和外部单位法线，当元素是时，该方法也会失败不规律。这部分代码如下所示：

IF (J==1) THEN

        xsd1(1:4)=xsd(1:4)
        ysd1(1:4)=ysd(1:4)
        zsd1(1:4)=zsd(1:4)
        b(1)=xsd1(1)-xsd1(3)
        b(2)=ysd1(1)-ysd1(3)
        b(3)=zsd1(1)-zsd1(3)
        c(1)=xsd1(2)-xsd1(4)
        c(2)=ysd1(2)-ysd1(4)
        c(3)=zsd1(2)-zsd1(4)        
        call cross(c, b,crossproduct)
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2        
        ELSE IF (J==2) THEN

        xsd2(1:4)=xsd(5:8)
        ysd2(1:4)=ysd(5:8)
        zsd2(1:4)=zsd(5:8)
        b(1)=xsd2(4)-xsd2(2)
        b(2)=ysd2(4)-ysd2(2)
        b(3)=zsd2(4)-zsd2(2)
        c(1)=xsd2(3)-xsd2(1)
        c(2)=ysd2(3)-ysd2(1)
        c(3)=zsd2(3)-zsd2(1)
        call cross(c, b,crossproduct)   
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5        
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2            
        ELSE IF (J==3) THEN

        xsd3(1:2)=xsd(1:2)
        xsd3(3)=xsd(6)
        xsd3(4)=xsd(5)
        ysd3(1:2)=ysd(1:2)
        ysd3(3)=ysd(6)
        ysd3(4)=ysd(5)      
        zsd3(1:2)=zsd(1:2)
        zsd3(3)=zsd(6)
        zsd3(4)=zsd(5)
        b(1)=xsd3(4)-xsd3(2)
        b(2)=ysd3(4)-ysd3(2)
        b(3)=zsd3(4)-zsd3(2)
        c(1)=xsd3(3)-xsd3(1)
        c(2)=ysd3(3)-ysd3(1)
        c(3)=zsd3(3)-zsd3(1)
        call cross(c, b,crossproduct)   
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5        
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2            
        ELSE IF (J==4) THEN

        xsd4(1:2)=xsd(7:8)
        xsd4(3)=xsd(4)
        xsd4(4)=xsd(3)
        ysd4(1:2)=ysd(7:8)
        ysd4(3)=ysd(4)
        ysd4(4)=ysd(3)
        zsd4(1:2)=zsd(7:8)
        zsd4(3)=zsd(4)
        zsd4(4)=zsd(3)
        b(1)=xsd4(1)-xsd4(3)
        b(2)=ysd4(1)-ysd4(3)
        b(3)=zsd4(1)-zsd4(3)
        c(1)=xsd4(2)-xsd4(4)
        c(2)=ysd4(2)-ysd4(4)
        c(3)=zsd4(2)-zsd4(4)
        call cross(c, b,crossproduct)   
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5        
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2            
        ELSE IF (J==5) THEN

        xsd5(1:2)=xsd(2:3)
        xsd5(3)=xsd(7)
        xsd5(4)=xsd(6)
        ysd5(1:2)=ysd(2:3)
        ysd5(3)=ysd(7)
        ysd5(4)=ysd(6)
        zsd5(1:2)=zsd(2:3)
        zsd5(3)=zsd(7)
        zsd5(4)=zsd(6)
        b(1)=xsd5(4)-xsd5(2)
        b(2)=ysd5(4)-ysd5(2)
        b(3)=zsd5(4)-zsd5(2)
        c(1)=xsd5(3)-xsd5(1)
        c(2)=ysd5(3)-ysd5(1)
        c(3)=zsd5(3)-zsd5(1)
        call cross(c, b,crossproduct)
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5        
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2        
        ELSE IF (J==6) THEN

        xsd6(1)=xsd(1)
        xsd6(2)=xsd(4)
        xsd6(3)=xsd(8)
        xsd6(4)=xsd(5)
        ysd6(1)=ysd(1)
        ysd6(2)=ysd(4)
        ysd6(3)=ysd(8)
        ysd6(4)=ysd(5)
        zsd6(1)=zsd(1)
        zsd6(2)=zsd(4)
        zsd6(3)=zsd(8)
        zsd6(4)=zsd(5)
        b(1)=xsd6(3)-xsd6(1)
        b(2)=ysd6(3)-ysd6(1)
        b(3)=zsd6(3)-zsd6(1)
        c(1)=xsd6(4)-xsd6(2)
        c(2)=ysd6(4)-ysd6(2)
        c(3)=zsd6(4)-zsd6(2)
        call cross(c, b,crossproduct)   
        norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5        
        area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2            
        END IF

J 是面，我遍历各个面中的点并找到对角线的叉积。

我想知道这两种方法有什么问题以及纠正方法。如果有人能给我一条前进的道路，我会很高兴。

1个回答

要调试这样的情况，将问题分解为几个部分是最简单的。由于您实际上并不需要等参数图来计算面积和法线，因此我将忽略您方法的这一方面。

第 1 步： 找到面上节点的方向。他们应该遵循右手规则（逆时针顺序）。您可以使用https://en.wikipedia.org/wiki/Curve_orientation中的方法（在 xy、yz 或 zx 平面上进行适当的投影）。例如

function isCounterClockwise(x_coords, y_coords, z_coords) 
  real   :: isCounterClockwise
  real   :: x_coords(4)
  real   :: y_coords(4)
  real   :: z_coords(4)

  isCounterClockwise = ....
  return

步骤 2：如果点不是逆时针方向，则颠倒顺序。

subroutine reverseOrder(x_coords, y_coords, z_coords)
 ....

第 3 步：计算叉积

subroutine computeAreaAndNormal(x_coords, y_coords, z_coords, area, normal)

  u1 = x_coords(2) - x_coords(1)
  u2 = y_coords(2) - y_coords(1)
  u3 = z_coords(2) - z_coords(1)

  v1 = x_coords(4) - x_coords(1)
  v2 = y_coords(4) - y_coords(1)
  v3 = z_coords(4) - z_coords(1)

  normal = crossProduct(u1, u2, u3, v1, v2, v3)
  area = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
  normal = normal/area

如果您以这种方式编写代码，则调试起来会容易得多。 警告：请注意，上述内容不可编译 Fortran。

更新：

对于节点 (1,2,3,4) 顺时针排列的任意平面四边形，求三角形 (1,2,3) 和 (1,3,4) 的面积（叉积范数的一半）和添加它们。

其它你可能感兴趣的问题

上一篇稍微改变两个向量来满足一个约束下一篇如何在 Matlab 中找到默认的 ODE 求解器容差？