From 6e7a80c1ef495411e5339232f691aefad9e5a1a7 Mon Sep 17 00:00:00 2001
From: s2006749 <s2006749@ed.ac.uk>
Date: Wed, 26 Jun 2024 14:01:20 +0100
Subject: [PATCH] added shear flow
---
src/BC-Ellipsoid.f90 | 25 +++++++++++++++++++++++--
src/module_param.f90 | 4 ++--
src/parameters.f90 | 2 +-
3 files changed, 26 insertions(+), 5 deletions(-)
diff --git a/src/BC-Ellipsoid.f90 b/src/BC-Ellipsoid.f90
index 022acf80..21847075 100644
--- a/src/BC-Ellipsoid.f90
+++ b/src/BC-Ellipsoid.f90
@@ -130,7 +130,7 @@ subroutine inflow (phi)
implicit none
- integer :: j,k,is
+ integer :: i,j,k,is
real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi
!call random_number(bxo)
@@ -144,6 +144,19 @@ subroutine inflow (phi)
enddo
enddo
+ if (shear_flow_ybc.eq.1) then
+ do k=1,xsize(3)
+ do i=1,xsize(1)
+ byxn(i,k)=+shear_velocity
+ enddo
+ enddo
+ do k=1,xsize(3)
+ do i=1,xsize(1)
+ byx1(i,k)=-shear_velocity
+ enddo
+ enddo
+ endif
+
if (iscalar.eq.1) then
do is=1, numscalar
do k=1,xsize(3)
@@ -245,7 +258,7 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ux1,uy1,uz1
real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
- real(mytype) :: y,um,eqr
+ real(mytype) :: y,um,eqr,temp_vel
integer :: k,j,i,ii,is,code
! write(*,*) 'INSIDE INIT ELLIP'
@@ -279,6 +292,8 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
ux1=zero; uy1=zero; uz1=zero
+
+
if (iin.ne.0) then
call system_clock(count=code)
if (iin.eq.2) code=0
@@ -293,6 +308,7 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
do j=1,xsize(2)
do i=1,xsize(1)
ux1(i,j,k)=init_noise*(ux1(i,j,k)-0.5)
+
uy1(i,j,k)=init_noise*(uy1(i,j,k)-0.5)
uz1(i,j,k)=init_noise*(uz1(i,j,k)-0.5)
enddo
@@ -319,6 +335,11 @@ subroutine init_ellip (ux1,uy1,uz1,phi1)
do j=1,xsize(2)
do i=1,xsize(1)
ux1(i,j,k)=ux1(i,j,k)+u1
+ if (shear_flow_ybc.eq.1) then
+ temp_vel=(real((j+xstart(2)-1-1),mytype)*dy-yly/2.)/(yly/2.0)*shear_velocity
+ ux1(i,j,k)=ux1(i,j,k)+temp_vel
+ ! write(*,*) "shear velocity at y = ", real((j+xstart(2)-1-1),mytype)*dy, " = ", temp_vel
+ endif
uy1(i,j,k)=uy1(i,j,k)
uz1(i,j,k)=uz1(i,j,k)
enddo
diff --git a/src/module_param.f90 b/src/module_param.f90
index 9d16bbe8..0f700ca2 100644
--- a/src/module_param.f90
+++ b/src/module_param.f90
@@ -620,8 +620,8 @@ module ibm_param
real(mytype) :: cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra,ubcx,ubcy,ubcz,rads,rho_s,ellip_m,c_air,cvl_scalar,grav_y,grav_x,grav_z
real(mytype) :: position(3),orientation(4),linearVelocity(3),angularVelocity(4),linearAcceleration(3),linearForce(3),torque(3),shape(3),inertia(3,3)
real(mytype) :: position_1(3),linearVelocity_1(3),orientation_1(4),angularVelocity_1(4)
- real(mytype) :: chord,thickness,omega, tconv2_sign
+ real(mytype) :: chord,thickness,omega, tconv2_sign, shear_velocity
integer :: inana ! Analytical BC as Input
- integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag
+ integer :: imove, nozdrift, force_csv, bodies_fixed, cube_flag, shear_flow_ybc
end module ibm_param
!############################################################################
diff --git a/src/parameters.f90 b/src/parameters.f90
index 2e50404e..a2e0e803 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -62,7 +62,7 @@ subroutine parameter(input_i3d)
NAMELIST /Tripping/ itrip,A_tr,xs_tr_tbl,ys_tr_tbl,ts_tr_tbl,x0_tr_tbl
NAMELIST /ibmstuff/ cex,cey,cez,shx,shy,shz,oriw,orii,orij,orik,lvx,lvy,lvz,avx,avy,avz,ra, &
nobjmax,nraf,nvol,iforces, cvl_scalar, npif, izap, ianal, imove, thickness, chord, omega , &
- ubcx,ubcy,ubcz,rads,rho_s, c_air, grav_x,grav_y,grav_z, nozdrift, force_csv, bodies_fixed, cube_flag, tconv2_sign
+ ubcx,ubcy,ubcz,rads,rho_s, c_air, grav_x,grav_y,grav_z, nozdrift, force_csv, bodies_fixed, cube_flag, tconv2_sign, shear_flow_ybc, shear_velocity
NAMELIST /ForceCVs/ xld, xrd, yld, yud, zld, zrd
NAMELIST /LMN/ dens1, dens2, prandtl, ilmn_bound, ivarcoeff, ilmn_solve_temp, &
massfrac, mol_weight, imultispecies, primary_species, &
--
GitLab