From 556e7a6afbc61bda39d0d411f1fa9ae43ffb0512 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Tue, 20 Aug 2024 14:23:58 +0100
Subject: [PATCH] fixed signs in torque face 6 integrals. Changed switch to
rotating frame
---
src/forces.f90 | 98 +++++++++++++++++++++++++-------------------------
1 file changed, 49 insertions(+), 49 deletions(-)
diff --git a/src/forces.f90 b/src/forces.f90
index bb68e9f0..6e79d157 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -643,9 +643,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx -uxmid*uymid*dx*dz
fcvy = fcvy -uymid*uymid*dx*dz
@@ -702,9 +702,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx +uxmid*uymid*dx*dz
fcvy = fcvy +uymid*uymid*dx*dz
@@ -761,9 +761,9 @@ contains
! write(*,*) 'Calculating force at left x boundary', [xm,ym,zm]
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx -uxmid*uxmid*del_y(j)*dz
@@ -819,9 +819,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx + uxmid*uxmid*del_y(j)*dz
@@ -883,9 +883,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx= fcvx +uxmid*uzmid*dx*dy
fcvy= fcvy +uymid*uzmid*dx*dy
@@ -944,9 +944,9 @@ contains
call crossProduct(angularVelocity,[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)
- 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)
+ 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)
fcvx= fcvx -uxmid*uzmid*dx*dy
fcvy= fcvy -uymid*uzmid*dx*dy
@@ -1269,8 +1269,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,[fac1,fac2,fac3],coriolis)
+ ! call centrifugal_force(angularVelocity, [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
@@ -1364,9 +1364,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uymid*dx*dz
fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uymid*dx*dz
@@ -1437,9 +1437,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
- uzmid = half*(uz1(i,j,k)+uz1(i+1,j,k)) - linearVelocity(3) - rotationalComponent(3)
+ 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)
fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uymid*dx*dz
fcvy = fcvy +(uzmid*radial(1)-uxmid*radial(3))*uymid*dx*dz
@@ -1509,9 +1509,9 @@ contains
radial = [xm,ym,zm]-position
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uxmid*del_y(j)*dz
fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uxmid*del_y(j)*dz
@@ -1581,9 +1581,9 @@ contains
!momentum flux
call crossProduct(angularVelocity,[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)
+ 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)
fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uxmid*del_y(j)*dz
@@ -1596,8 +1596,8 @@ contains
prmid = half*(ppi2(i,j,k)+ppi2(i,j+1,k))
! fprx = fprx -prmid*del_y(j)*dz*(radial(3)-radial(2))
- fpry = fpry -prmid*del_y(j)*dz*(radial(3))
- fprz = fprz -prmid*del_y(j)*dz*(-radial(2))
+ fpry = fpry - prmid*del_y(j)*dz*(radial(3))
+ fprz = fprz - prmid*del_y(j)*dz*(-radial(2))
!viscous term
dudxmid = half*(ta2(i,j,k)+ta2(i,j+1,k))
@@ -1660,13 +1660,13 @@ contains
radial = [xm,ym,zm]-position
call crossProduct(angularVelocity,[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)
+ 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)
- fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j)
- fcvy = fcvy +(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
- fcvz = fcvz +(uxmid*radial(2)-uymid*radial(1))*uzmid*dx*del_y(j)
+ fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j) !!!CHANGE
+ fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
+ fcvz = fcvz -(uxmid*radial(2)-uymid*radial(1))*uzmid*dx*del_y(j)
!pressure
prmid = half*(ppi1(i,j,k)+ppi1(i+1,j,k))
@@ -1681,9 +1681,9 @@ contains
dwdymid = half*(tf1(i,j,k)+tf1(i+1,j,k))
dwdzmid = half*(ti1(i,j,k)+ti1(i+1,j,k))
- fdix = fdix + (xnu*(dvdzmid+dwdymid)*radial(3)-xnu*(two*dwdzmid)*radial(2))*dx*dy
- fdiy = fdiy + (xnu*(two*dwdzmid)*radial(1)-xnu*(dudzmid+dwdxmid)*radial(3))*dx*dy
- fdiz = fdiz + (xnu*(dudzmid+dwdxmid)*radial(2)-xnu*(dvdzmid+dwdymid)*radial(1))*dx*dy
+ fdix = fdix - (xnu*(dvdzmid+dwdymid)*radial(3)-xnu*(two*dwdzmid)*radial(2))*dx*dy
+ fdiy = fdiy - (xnu*(two*dwdzmid)*radial(1)-xnu*(dudzmid+dwdxmid)*radial(3))*dx*dy
+ fdiz = fdiz - (xnu*(dudzmid+dwdxmid)*radial(2)-xnu*(dvdzmid+dwdymid)*radial(1))*dx*dy
! fdix = fdix +(xnu*(dudzmid+dwdxmid)*dx*dy)
@@ -1734,13 +1734,13 @@ contains
call crossProduct(angularVelocity,[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)
- 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)
+ 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)
- fcvx = fcvx -(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j)
- fcvy = fcvy -(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
- fcvz = fcvz -(uxmid*radial(2)-uymid*radial(1))*uzmid*dx*del_y(j)
+ fcvx = fcvx +(uymid*radial(3)-uzmid*radial(2))*uzmid*dx*del_y(j)
+ fcvy = fcvy +(uzmid*radial(1)-uxmid*radial(3))*uzmid*dx*del_y(j)
+ fcvz = fcvz +(uxmid*radial(2)-uymid*radial(1))*uzmid*dx*del_y(j)
!pressure
--
GitLab