From 4a8db0fd1d0280980b584f05e77547147ed0fb01 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Mon, 9 Sep 2024 12:30:45 +0100
Subject: [PATCH] fixed ibm points and angularVelocity handling
---
src/ellip_utils.f90 | 7 +--
src/forces.f90 | 28 +++++-----
src/ibm.f90 | 125 ++++++++++++++++++++++++++++++++++++++++----
3 files changed, 132 insertions(+), 28 deletions(-)
diff --git a/src/ellip_utils.f90 b/src/ellip_utils.f90
index 2a517075..db683465 100644
--- a/src/ellip_utils.f90
+++ b/src/ellip_utils.f90
@@ -241,7 +241,7 @@ contains
end subroutine CrossProduct
subroutine CalculatePointVelocity(point, center, angularVelocity, linearVelocity, pointVelocity)
- real(mytype), intent(in) :: point(3), center(3), linearVelocity(3), angularVelocity(3)
+ real(mytype), intent(in) :: point(3), center(3), linearVelocity(3), angularVelocity(4)
real(mytype), intent(out) :: pointVelocity(3)
real(mytype) :: crossed(3)
! Compute the distance vector from the center to the point
@@ -250,7 +250,8 @@ contains
! Compute the cross product of angular velocity and distance vector
- call CrossProduct(distance, angularVelocity, crossed)
+ call CrossProduct(distance, angularVelocity(2:4), crossed)
+
! Calculate the velocity at the point
pointVelocity = crossed + linearVelocity
@@ -271,7 +272,7 @@ contains
subroutine navierFieldGen(center, linearVelocity, angularVelocity, ep1, ep1_x, ep1_y, ep1_z)
use param
use decomp_2d
- real(mytype), intent(in) :: center(3), linearVelocity(3), angularVelocity(3)
+ real(mytype), intent(in) :: center(3), linearVelocity(3), angularVelocity(4)
real(mytype), dimension(xsize(1),xsize(2),xsize(3)), intent(in) :: ep1
real(mytype),dimension(xsize(1),xsize(2),xsize(3)),intent(out) :: ep1_x, ep1_y, ep1_z
real(mytype) :: xm, ym, zm, point(3), x_pv, y_pv, z_pv, pointVelocity(3)
diff --git a/src/forces.f90 b/src/forces.f90
index 2895f3e6..74f60675 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -552,8 +552,8 @@ contains
fac2 = (onepfive*uy1(i,j,k)-two*uy01(i,j,k)+half*uy11(i,j,k))*(one-ep1(i,j,k))
fac3 = (onepfive*uz1(i,j,k)-two*uz01(i,j,k)+half*uz11(i,j,k))*(one-ep1(i,j,k))
- call coriolis_force(angularVelocity,[fac1,fac2,fac3],coriolis)
- call centrifugal_force(angularVelocity, [xm,ym,zm]-position,centrifugal)
+ call coriolis_force(angularVelocity(2:4),[fac1,fac2,fac3],coriolis)
+ call centrifugal_force(angularVelocity(2:4), [xm,ym,zm]-position,centrifugal)
! The velocity time rate has to be relative to the cell center,
! and not to the nodes, because, here, we have an integral
! relative to the volume, and, therefore, this has a sense
@@ -642,7 +642,7 @@ contains
! write(*,*) 'Calculating force at upper y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -701,7 +701,7 @@ contains
! write(*,*) 'Calculating force at lower y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -760,7 +760,7 @@ contains
ym=real(jj,mytype)*dz
! write(*,*) 'Calculating force at left x boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -818,7 +818,7 @@ contains
! write(*,*) 'Calculating force at right x boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -882,7 +882,7 @@ contains
! write(*,*) 'Calculating force at left z boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -941,7 +941,7 @@ contains
ii=xstart(1)+i-1
xm=real(ii,mytype)*dx
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
! write(*,*) 'Calculating force at right z boundary', [xm,ym,zm]
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
@@ -1363,7 +1363,7 @@ contains
radial = [xm,ym,zm]-position
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -1435,7 +1435,7 @@ contains
! write(*,*) 'Calculating force at lower y boundary', [xm,ym,zm]
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
radial = [xm,ym,zm]-position
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
@@ -1508,7 +1508,7 @@ contains
!momentum flux
radial = [xm,ym,zm]-position
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -1580,7 +1580,7 @@ contains
radial = [xm,ym,zm]-position
!momentum flux
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux2(i,j,k)+ux2(i,j+1,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy2(i,j,k)+uy2(i,j+1,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz2(i,j,k)+uz2(i,j+1,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -1659,7 +1659,7 @@ contains
!momentum flux
radial = [xm,ym,zm]-position
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
uymid = half*(uy1(i,j,k)+uy1(i+1,j,k)) - linearVelocity(2) + rotationalComponent(2)
uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) + rotationalComponent(3)
@@ -1731,7 +1731,7 @@ contains
!momentum flux
radial = [xm,ym,zm]-position
- call crossProduct(angularVelocity,[xm,ym,zm]-position,rotationalComponent)
+ call crossProduct(angularVelocity(2:4),[xm,ym,zm]-position,rotationalComponent)
! write(*,*) 'Calculating force at right z boundary', [xm,ym,zm]
uxmid = half*(ux1(i,j,k)+ux1(i+1,j,k)) - linearVelocity(1) + rotationalComponent(1)
diff --git a/src/ibm.f90 b/src/ibm.f90
index 3f4cc296..8a389b57 100644
--- a/src/ibm.f90
+++ b/src/ibm.f90
@@ -402,7 +402,7 @@ subroutine cubsplx(u,lind)
USE decomp_2d
USE variables
USE ibm_param
- USE ellipsoid_utils, ONLY: CalculatePointVelocity
+ USE ellipsoid_utils, ONLY: CalculatePointVelocity, EllipsoidalRadius
!
implicit none
!
@@ -420,7 +420,7 @@ subroutine cubsplx(u,lind)
integer :: inxi,inxf
real(mytype) :: ana_resi,ana_resf ! Position of Boundary (Analytically)
real(mytype) :: point(3),pointVelocity(3)
- real(mytype) :: xm,ym,zm,x_pv,y_pv,z_pv
+ real(mytype) :: xm,ym,zm,x_pv,y_pv,z_pv, dummy
!
! Initialise Arrays
xa(:)=0.
@@ -430,15 +430,18 @@ subroutine cubsplx(u,lind)
! bcimp=lind
! write(*,*) lind
!
+ dummy=0.0
do k=1,xsize(3)
- zm=real(xstart(3)+k-2,mytype)*dz
+ zm=real(xstart(3)+k-1,mytype)*dz
do j=1,xsize(2)
- ym=real(xstart(2)+j-2,mytype)*dy
+ ym=real(xstart(2)+j-1,mytype)*dy
if(nobjx(j,k).ne.0)then
ia=0
do i=1,nobjx(j,k)
! 1st Boundary - I DON'T UNDERSTAND THIS XM CONVERSION.
- xm=real(xstart(1)+i-2,mytype)*dx
+ ! write(*,*) "Nobjx = ", nobjx(j,k)
+ ! xm=real(xstart(1)+i-2,mytype)*dx
+ xm = xi(i,j,k)
nxpif=npif
ia=ia+1
if (ianal.eq.0) then
@@ -450,6 +453,14 @@ subroutine cubsplx(u,lind)
endif
point=[xm,ym,zm]
call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
+ ! write(*,*) "Called CPV at ", point, "returned", pointVelocity
+ ! call EllipsoidalRadius(point, position, orientation, shape, dummy)
+ ! dummy = maxval(abs((point(2:3)-position(2:3))))
+ ! if ((dummy.gt.(1.01)).or.(dummy.lt.0.99)) then
+ ! write(*,*) "At ", point, "r = ", dummy
+ ! endif
+ ! write(*,*) "Radius = ", dummy
+ ! write(*,*) "radius = ", dummy, "At point", point
x_pv=pointVelocity(1)
y_pv=pointVelocity(2)
z_pv=pointVelocity(3)
@@ -519,6 +530,42 @@ subroutine cubsplx(u,lind)
call analitic_x(j,xf(i,j,k),ana_resf,k) ! Calculate the position of BC analytically
xa(ia)=ana_resf
endif
+ xm = xf(i,j,k)
+ point=[xm,ym,zm]
+ call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
+ ! write(*,*) "Called CPV at ", point, "returned", pointVelocity
+ ! call EllipsoidalRadius(point, position, orientation, shape, dummy)
+ ! dummy = maxval(abs(point-position))
+ ! write(*,*) "Radius = ", dummy
+ ! dummy = maxval(abs(point-position))
+ ! if ((dummy.gt.(1.01)).or.(dummy.lt.0.99)) then
+ ! write(*,*) "At ", point, "r = ", dummy
+ ! endif
+ x_pv=pointVelocity(1)
+ y_pv=pointVelocity(2)
+ z_pv=pointVelocity(3)
+ if (lind.eq.0) then
+ bcimp=zero
+ elseif (lind.eq.1) then
+ bcimp=x_pv
+ ! write(*,*) bcimp
+ elseif (lind.eq.2) then
+ bcimp=y_pv
+ elseif (lind.eq.3) then
+ bcimp=z_pv
+ elseif (lind.eq.4) then
+ bcimp=x_pv*x_pv
+ elseif (lind.eq.5) then
+ bcimp=y_pv*y_pv
+ elseif (lind.eq.6) then
+ bcimp=z_pv*z_pv
+ elseif (lind.eq.7) then
+ bcimp=x_pv*y_pv
+ elseif (lind.eq.8) then
+ bcimp=x_pv*z_pv
+ elseif (lind.eq.9) then
+ bcimp=y_pv*z_pv
+ endif
ya(ia)=bcimp
if(xf(i,j,k).lt.xlx)then ! Immersed Object
inxf=0
@@ -624,13 +671,13 @@ subroutine cubsply(u,lind)
!
do k=1,ysize(3)
- zm=real(ystart(3)+k-2,mytype)*dz
+ zm=real(ystart(3)+k-1,mytype)*dz
do i=1,ysize(1)
- xm=real(ystart(1)+i-2,mytype)*dx
+ xm=real(ystart(1)+i-1,mytype)*dx
if(nobjy(i,k).ne.0)then
ia=0
do j=1,nobjy(i,k)
- ym=real(ystart(2)+j-2,mytype)*dy
+ ! ym=real(ystart(2)+j-2,mytype)*dy
! 1st Boundary
nypif=npif
ia=ia+1
@@ -641,6 +688,7 @@ subroutine cubsply(u,lind)
call analitic_y(i,yi(j,i,k),ana_resi,k) ! Calculate the position of BC analytically
xa(ia)=ana_resi
endif
+ ym = yi(j,i,k)
point=[xm,ym,zm]
call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
x_pv=pointVelocity(1)
@@ -716,6 +764,33 @@ subroutine cubsply(u,lind)
call analitic_y(i,yf(j,i,k),ana_resf,k) ! Calculate the position of BC analytically
xa(ia)=ana_resf
endif
+ ym = yf(j,i,k)
+ point=[xm,ym,zm]
+ call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
+ x_pv=pointVelocity(1)
+ y_pv=pointVelocity(2)
+ z_pv=pointVelocity(3)
+ if (lind.eq.0) then
+ bcimp=zero
+ elseif (lind.eq.1) then
+ bcimp=x_pv
+ elseif (lind.eq.2) then
+ bcimp=y_pv
+ elseif (lind.eq.3) then
+ bcimp=z_pv
+ elseif (lind.eq.4) then
+ bcimp=x_pv*x_pv
+ elseif (lind.eq.5) then
+ bcimp=y_pv*y_pv
+ elseif (lind.eq.6) then
+ bcimp=z_pv*z_pv
+ elseif (lind.eq.7) then
+ bcimp=x_pv*y_pv
+ elseif (lind.eq.8) then
+ bcimp=x_pv*z_pv
+ elseif (lind.eq.9) then
+ bcimp=y_pv*z_pv
+ endif
ya(ia)=bcimp
if(yf(j,i,k).lt.yly)then ! Immersed Object
jy=1
@@ -822,14 +897,14 @@ subroutine cubsplz(u,lind)
!
do j=1,zsize(2)
- ym=real(zstart(2)+j-2,mytype)*dy
+ ym=real(zstart(2)+j-1,mytype)*dy
do i=1,zsize(1)
- xm=real(zstart(1)+i-2,mytype)*dx
+ xm=real(zstart(1)+i-1,mytype)*dx
if(nobjz(i,j).ne.0)then
ia=0
do k=1,nobjz(i,j)
! 1st Boundary
- zm=real(zstart(3)+k-2,mytype)*dz
+ ! zm=real(zstart(3)+k-2,mytype)*dz
nzpif=npif
ia=ia+1
if (ianal.eq.0) then
@@ -839,6 +914,7 @@ subroutine cubsplz(u,lind)
! call analitic_z(i,zi(k,i,j),ana_resi,j) ! Calculate the position of BC analytically
xa(ia)=ana_resi
endif
+ zm = zi(k,i,j)
point=[xm,ym,zm]
call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
x_pv=pointVelocity(1)
@@ -907,6 +983,33 @@ subroutine cubsplz(u,lind)
!call analitic_z(i,zf(k,i,j),ana_resf,j) ! Calculate the position of BC analytically
xa(ia)=ana_resf
endif
+ zm = zi(k,i,j)
+ point=[xm,ym,zm]
+ call CalculatePointVelocity(point, position, angularVelocity, linearVelocity, pointVelocity)
+ x_pv=pointVelocity(1)
+ y_pv=pointVelocity(2)
+ z_pv=pointVelocity(3)
+ if (lind.eq.0) then
+ bcimp=zero
+ elseif (lind.eq.1) then
+ bcimp=x_pv
+ elseif (lind.eq.2) then
+ bcimp=y_pv
+ elseif (lind.eq.3) then
+ bcimp=z_pv
+ elseif (lind.eq.4) then
+ bcimp=x_pv*x_pv
+ elseif (lind.eq.5) then
+ bcimp=y_pv*y_pv
+ elseif (lind.eq.6) then
+ bcimp=z_pv*z_pv
+ elseif (lind.eq.7) then
+ bcimp=x_pv*y_pv
+ elseif (lind.eq.8) then
+ bcimp=x_pv*z_pv
+ elseif (lind.eq.9) then
+ bcimp=y_pv*z_pv
+ endif
ya(ia)=bcimp
if(zf(k,i,j).lt.zlz)then ! Immersed Object
inxf=0
--
GitLab