From fd53dbd98ef5155c5ee1322f0c0a09646237331b Mon Sep 17 00:00:00 2001 From: Semih Akkurt Date: Fri, 5 Jul 2024 10:55:11 +0100 Subject: [PATCH] Add actuator disc model. --- src/CMakeLists.txt | 1 + src/solver.f90 | 56 ++++++++++++- src/windturb_adm.f90 | 190 +++++++++++++++++++++++++++++++++++++++++++ src/xcompact.f90 | 7 +- 4 files changed, 249 insertions(+), 5 deletions(-) create mode 100644 src/windturb_adm.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b242ad21..29be2ea8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,6 +9,7 @@ set(SRC ordering.f90 mesh.f90 field.f90 + windturb_adm.f90 omp/backend.f90 omp/common.f90 omp/kernels/distributed.f90 diff --git a/src/solver.f90 b/src/solver.f90 index 0ce9a46a..81d831fc 100644 --- a/src/solver.f90 +++ b/src/solver.f90 @@ -11,6 +11,7 @@ module m_solver use m_tdsops, only: tdsops_t, dirps_t use m_time_integrator, only: time_intg_t use m_mesh, only: mesh_t + use m_windturb_adm, only: windturb_adm_t implicit none @@ -45,6 +46,7 @@ module m_solver real(dp) :: dt, nu integer :: n_iters, n_output integer :: ngrid + logical :: wind_turbine class(field_t), pointer :: u, v, w @@ -54,11 +56,13 @@ module m_solver type(allocator_t), pointer :: host_allocator class(dirps_t), pointer :: xdirps, ydirps, zdirps procedure(poisson_solver), pointer :: poisson => null() + type(windturb_adm_t) :: windturb_adm contains procedure :: transeq procedure :: divergence_v2p procedure :: gradient_p2v procedure :: curl + procedure :: compute_turbines procedure :: output procedure :: run end type solver_t @@ -99,6 +103,7 @@ function init(backend, mesh, time_integrator, host_allocator, & integer :: i, j, k integer, dimension(3) :: dims real(dp), dimension(3) :: xloc + real(dp), allocatable, dimension(:, :) :: disc_params solver%backend => backend solver%mesh => mesh @@ -132,8 +137,8 @@ function init(backend, mesh, time_integrator, host_allocator, & y = xloc(2) z = xloc(3) - u_init%data(i, j, k) = sin(x)*cos(y)*cos(z) - v_init%data(i, j, k) = -cos(x)*sin(y)*cos(z) + u_init%data(i, j, k) = 10._dp!sin(x)*cos(y)*cos(z) + v_init%data(i, j, k) = 0._dp!-cos(x)*sin(y)*cos(z) w_init%data(i, j, k) = 0 end do end do @@ -163,6 +168,19 @@ function init(backend, mesh, time_integrator, host_allocator, & solver%poisson => poisson_cg end select + solver%wind_turbine = .true. + if (solver%wind_turbine) then + allocate (disc_params(2, 8)) + ! coords x, y, z, yaw, tilt, D, C_T, alpha + disc_params(1, :) = [500._dp, 500._dp, 500._dp, 0._dp, 0._dp, & + 100._dp, 0.75_dp, 0.25_dp] + disc_params(2, :) = [1500._dp, 500._dp, 500._dp, 0._dp, 0._dp, & + 100._dp, 0.75_dp, 0.25_dp] + ! instantiate the actuator disc class + solver%windturb_adm = windturb_adm_t(backend, mesh, host_allocator, & + disc_params) + end if + end function init subroutine allocate_tdsops(dirps, dir, backend) @@ -542,6 +560,32 @@ subroutine curl(self, o_i_hat, o_j_hat, o_k_hat, u, v, w) end subroutine curl + subroutine compute_turbines(self, du, dv, dw, u, v, w) + implicit none + + class(solver_t) :: self + class(field_t), intent(inout) :: du, dv, dw + class(field_t), intent(in) :: u, v, w + + class(field_t), pointer :: Fx, Fy, Fz + + Fx => self%backend%allocator%get_block(DIR_X) + Fy => self%backend%allocator%get_block(DIR_X) + Fz => self%backend%allocator%get_block(DIR_X) + + call self%windturb_adm%compute_sources(Fx, Fy, Fz, u, v, w) + + ! add resulting forces into the rhs of the momentum equation + !call self%backend%vecadd(self%windturb_adm%rho_air, Fx, 1._dp, du) + !call self%backend%vecadd(self%windturb_adm%rho_air, Fy, 1._dp, dv) + !call self%backend%vecadd(self%windturb_adm%rho_air, Fz, 1._dp, dw) + + call self%backend%allocator%release_block(Fx) + call self%backend%allocator%release_block(Fy) + call self%backend%allocator%release_block(Fz) + + end subroutine compute_turbines + subroutine poisson_fft(self, pressure, div_u) implicit none @@ -656,6 +700,11 @@ subroutine run(self) call self%transeq(du, dv, dw, self%u, self%v, self%w) + ! calculate the source terms from turbines and add into du, dv, dw + if (self%wind_turbine) then + call self%compute_turbines(du, dv, dw, self%u, self%v, self%w) + end if + ! time integration call self%time_integrator%step(self%u, self%v, self%w, & du, dv, dw, self%dt) @@ -696,6 +745,9 @@ subroutine run(self) if (mod(i, self%n_output) == 0) then t = i*self%dt call self%output(t) + if (self%mesh%par%is_root()) then + print *, 'power', self%windturb_adm%disc(:)%power + end if end if end do diff --git a/src/windturb_adm.f90 b/src/windturb_adm.f90 new file mode 100644 index 00000000..f6c8c9fd --- /dev/null +++ b/src/windturb_adm.f90 @@ -0,0 +1,190 @@ +module m_windturb_adm + use iso_fortran_env, only: stderr => error_unit + use mpi + + use m_allocator, only: allocator_t, field_t + use m_base_backend, only: base_backend_t + use m_common, only: dp, pi, DIR_X, DIR_C, VERT + use m_field, only: field_t + use m_mesh, only: mesh_t + + implicit none + + type :: actuator_disc_t + integer :: id + real(dp) :: coords(3) ! center of rotation + real(dp) :: yaw ! rotor yaw angle, in radians w.r.t. y-axis + real(dp) :: tilt ! rotor tilt angle, in radians w.r.t. z-axis + real(dp) :: D, area ! actuator disk diameter and area + real(dp) :: C_T ! thrust coefficient + real(dp) :: alpha ! induction coefficient + real(dp) :: rot_N(3) ! axis of rotation + real(dp) :: U_disc, U_disc_filt ! disk-averaged speed and filtered speed + real(dp) :: power, thrust ! instantaneous quantities + real(dp) :: U_disc_tavg, power_tavg, thrust_tavg ! time averaged quantities + class(field_t), pointer :: gamma_disc + end type actuator_disc_t + + type :: windturb_adm_t + integer :: n_turb + real(dp) :: rho_air, cell_vol + type(actuator_disc_t), allocatable :: disc(:) + class(base_backend_t), pointer :: backend + class(mesh_t), pointer :: mesh + type(allocator_t), pointer :: host_allocator + contains + procedure :: compute_sources + end type windturb_adm_t + + interface windturb_adm_t + module procedure init + end interface windturb_adm_t + +contains + + function init(backend, mesh, host_allocator, disc_params) & + result(windturb_adm) + implicit none + + class(base_backend_t), target, intent(inout) :: backend + type(mesh_t), target, intent(inout) :: mesh + type(allocator_t), target, intent(inout) :: host_allocator + real(dp), dimension(:, :), intent(in) :: disc_params + type(windturb_adm_t) :: windturb_adm + + class(field_t), pointer :: gamma_host + real(dp) :: disc_coords(3), yaw, tilt, D, rot_N(3) + integer :: t, i, j, k, dims(3), ierr + real(dp) :: coords(3), delta, dx, dy, dz, delta_N, delta_R, disc_thick, & + gamma_val, gamma_tot + + windturb_adm%backend => backend + windturb_adm%mesh => mesh + windturb_adm%host_allocator => host_allocator + + windturb_adm%n_turb = size(disc_params, dim=1) + windturb_adm%rho_air = 1._dp + windturb_adm%cell_vol = product(mesh%geo%d) + + allocate (windturb_adm%disc(windturb_adm%n_turb)) + + dims = mesh%get_dims(VERT) + + gamma_host => windturb_adm%host_allocator%get_block(DIR_C) + + do t = 1, windturb_adm%n_turb + disc_coords(:) = disc_params(t, 1:3) + yaw = disc_params(t, 4)*pi/180._dp + tilt = disc_params(t, 5)*pi/180._dp + D = disc_params(t, 6) + rot_N(1) = cos(yaw)*cos(tilt) + rot_N(2) = sin(tilt) + rot_N(3) = sin(yaw) + windturb_adm%disc(t)%coords(:) = disc_coords(:) + windturb_adm%disc(t)%yaw = yaw + windturb_adm%disc(t)%tilt = tilt + windturb_adm%disc(t)%D = D + windturb_adm%disc(t)%C_T = disc_params(t, 7) + windturb_adm%disc(t)%alpha = disc_params(t, 8) + windturb_adm%disc(t)%rot_N(:) = rot_N(:) + windturb_adm%disc(t)%area = pi*(D**2)/4._dp + + gamma_host%data(:, :, :) = 0._dp + gamma_tot = 0._dp + + delta = sqrt((mesh%geo%d(1)*rot_N(1))**2 & + + (mesh%geo%d(2)*rot_N(2))**2 & + + (mesh%geo%d(3)*rot_N(3))**2) + disc_thick = max(D/8._dp, delta*1.5_dp) + do k = 1, dims(3) + do j = 1, dims(2) + do i = 1, dims(1) + coords = mesh%get_coordinates(i, j, k) + dx = coords(1) - disc_coords(1) + dy = coords(2) - disc_coords(2) + dz = coords(3) - disc_coords(3) + delta_N = dx*rot_N(1) + dy*rot_N(2) - dz*rot_N(3) + + dx = coords(1) - delta_N*rot_N(1) + dy = coords(2) - delta_N*rot_N(2) + dz = coords(3) + delta_N*rot_N(3) + delta_R = sqrt((dx - disc_coords(1))**2 & + + (dy - disc_coords(2))**2 & + + (dz - disc_coords(3))**2) + + gamma_val = exp(-((delta_N/(disc_thick/2._dp))**2 & + + (delta_R/(D/2._dp))**8)) + gamma_host%data(i, j, k) = gamma_val + gamma_tot = gamma_tot + gamma_val + end do + end do + end do + + ! normalise the gamma field + call MPI_Allreduce(MPI_IN_PLACE, gamma_tot, 1, MPI_DOUBLE_PRECISION, & + MPI_SUM, MPI_COMM_WORLD, ierr) + gamma_host%data(:, :, :) = gamma_host%data(:, :, :)/gamma_tot + + windturb_adm%disc(t)%gamma_disc & + => windturb_adm%backend%allocator%get_block(DIR_X) + call windturb_adm%backend%set_field_data( & + windturb_adm%disc(t)%gamma_disc, gamma_host%data & + ) + end do + + call windturb_adm%host_allocator%release_block(gamma_host) + + end function init + + subroutine compute_sources(self, Fx, Fy, Fz, u, v, w) + implicit none + + class(windturb_adm_t) :: self + class(field_t), intent(inout) :: Fx, Fy, Fz + class(field_t), intent(in) :: u, v, w + + real(dp) :: u_avg, v_avg, w_avg, rot_N(3) + real(dp) :: coeff_x, coeff_y, coeff_z, C_T_prime + real(dp), allocatable, dimension(:) :: vel_avg + integer :: t, ierr + + allocate (vel_avg(self%n_turb)) + + do t = 1, self%n_turb + rot_N(:) = self%disc(t)%rot_N(:) + u_avg = self%backend%scalar_product(self%disc(t)%gamma_disc, u) + v_avg = self%backend%scalar_product(self%disc(t)%gamma_disc, v) + w_avg = self%backend%scalar_product(self%disc(t)%gamma_disc, w) + vel_avg(t) = u_avg*rot_N(1) + v_avg*rot_N(2) - w_avg*rot_N(3) + end do + + ! obtain the average velocity combining all the subdomains + call MPI_Allreduce(vel_avg, self%disc%U_disc, self%n_turb, & + MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr) + + deallocate (vel_avg) + + ! no filtering for the moment + self%disc(:)%U_disc_filt = self%disc(:)%U_disc + + ! use a hack to zeroise the force fields + call self%backend%vecadd(0._dp, Fx, 0._dp, Fx) + call self%backend%vecadd(0._dp, Fy, 0._dp, Fy) + call self%backend%vecadd(0._dp, Fz, 0._dp, Fz) + + do t = 1, self%n_turb + C_T_prime = self%disc(t)%C_T/(1._dp - self%disc(t)%alpha)**2 + self%disc(t)%thrust = 0.5_dp*self%rho_air*C_T_prime & + *self%disc(t)%U_disc_filt**2*self%disc(t)%area + self%disc(t)%power = self%disc(t)%thrust*self%disc(t)%U_disc_filt + coeff_x = -self%disc(t)%thrust*self%disc(t)%rot_N(1)/self%cell_vol + coeff_y = -self%disc(t)%thrust*self%disc(t)%rot_N(2)/self%cell_vol + coeff_z = self%disc(t)%thrust*self%disc(t)%rot_N(3)/self%cell_vol + call self%backend%vecadd(coeff_x, self%disc(t)%gamma_disc, 1._dp, Fx) + call self%backend%vecadd(coeff_y, self%disc(t)%gamma_disc, 1._dp, Fy) + call self%backend%vecadd(coeff_z, self%disc(t)%gamma_disc, 1._dp, Fz) + end do + + end subroutine compute_sources + +end module m_windturb_adm diff --git a/src/xcompact.f90 b/src/xcompact.f90 index 272be196..9944b6ba 100644 --- a/src/xcompact.f90 +++ b/src/xcompact.f90 @@ -60,18 +60,19 @@ program xcompact #endif ! Global number of cells in each direction - dims_global = [256, 256, 256] + dims_global = [512, 256, 256] ! Global domain dimensions L_global = [2*pi, 2*pi, 2*pi] + L_global = [2000, 1000, 1000] ! Domain decomposition in each direction nproc_dir = [1, 1, nproc] mesh = mesh_t(dims_global, nproc_dir, L_global) - globs%dt = 0.001_dp - globs%nu = 1._dp/1600._dp + globs%dt = 0.1_dp + globs%nu = 1._dp/66000._dp globs%n_iters = 20000 globs%n_output = 100