#!/usr/bin/env python
# -*- coding: utf-8 -*-
# This file is part of the Cortix toolkit environment
# https://cortix.org
import logging
import numpy as np
import scipy.constants as const
from collections import namedtuple
import matplotlib
matplotlib.use('Agg', warn=False)
import matplotlib.pyplot as plt
from cortix.src.module import Module
from cortix.support.phase import Phase
from cortix.support.specie import Specie
from cortix.support.quantity import Quantity
[docs]class Vortex(Module):
'''
Vortex module used to model fluid flow using Cortix.
Notes
-----
Any `port` name and any number of ports are allowed.
'''
[docs] def __init__(self):
'''
Attributes
----------
initial_time: float
end_time: float
time_step: float
show_time: tuple
Two-element tuple, `(bool,float)`, `True` will print to standard
output.
'''
super().__init__()
species = list()
quantities = list()
self.initial_time = 0.0
self.end_time = 5*const.minute
self.time_step = 0.1
self.show_time = (False,1*const.minute)
self.log = logging.getLogger('cortix')
air = Specie(name='air', formula_name='Air', phase='gas')
air.massCCUnit = 'g/cc'
air.molarCCUnit = 'mole/cc'
air.molarMass = 0.3 * 16 * 2 + 0.7 * 14 *2
species.append(air)
# Constant values for the vortex fluid.
self.mass_density = 0.1 * const.gram / const.centi**3 # [kg/m^3]
self.dyn_viscosity = 1.81e-5 # kg/(m s)
# Domain box dimensions: LxLxH m^3 box with given H.
# z coordinate pointing upwards. -L <= x <= L, -L <= y <= L,
# z component is positive => vortex is blowing upwards.
self.box_half_length = 250.0 # vortex box [m]
self.box_height = 500.0 # [m]
# Vortex parameters.
self.min_core_radius = 2.5 # [m]
self.outer_v_theta = 1.0 # m/s # angular speed
self.v_z_0 = 0.50 # [m/s]
self.period = 20 # wind change period
[docs] def run(self, *args):
# namedtuple does not pickle into send message; investigate later: vfda TODO
#Props = namedtuple('Props',['mass_density','dyn_viscosity'])
#fluid_props = Props( self.mass_density, self.dyn_viscosity )
fluid_props = ( self.mass_density, self.dyn_viscosity )
time = self.initial_time
while time < self.end_time:
if self.show_time[0] and abs(time%self.show_time[1]-0.0)<=1.e-1:
self.log.info('Vortex::time[min] = '+str(round(time/const.minute,1)))
# Interactions in all nameless ports (lower level port send/recv used)
#---------------------------------------------------------------------
for port in self.ports:
(message_time, position) = port.recv()
# Compute the vortex velocity using the given position
velocity = self.compute_velocity(message_time, position)
# Send the vortex velocity to caller
port.send( (message_time, velocity, fluid_props) )
time += self.time_step
return
[docs] def compute_velocity(self, time, position):
'''
Compute the vortex velocity at the given external
position using a vortex flow model
Parameters
----------
time: float
Time in SI unit.
position: numpy.ndarray(3)
Spatial position in SI unit.
Returns
-------
vortex_velocity: numpy.ndarray(3)
'''
import math
cycle_freq = 1/self.period
radian_freq = 2*math.pi*cycle_freq
outer_cylindrical_radius = np.hypot(self.box_half_length, self.box_half_length)
circulation = 2 * np.pi * outer_cylindrical_radius * self.outer_v_theta # m^2/s
core_radius = self.min_core_radius
x = position[0]
y = position[1]
z = position[2]
relax_length = self.box_height / 2.0
z_relax_factor = np.exp(-(self.box_height-z)/relax_length)
v_z = self.v_z_0 * z_relax_factor * abs(math.cos( radian_freq * time))
cylindrical_radius = np.hypot(x,y)
azimuth = np.arctan2(y,x)
v_theta = (1 - np.exp(-cylindrical_radius**2 / 8 / core_radius**2)) *\
circulation / 2 / np.pi /\
max(cylindrical_radius,self.min_core_radius) *\
z_relax_factor * abs(math.cos( radian_freq * time))
v_x = - v_theta * np.sin(azimuth)
v_y = v_theta * np.cos(azimuth)
return np.array([v_x,v_y,v_z])
[docs] def plot_velocity(self, time=None):
'''
Plot the vortex velocity as a function of height.
'''
if time is None:
time = self.initial_time
(fig,axs) = plt.subplots(2,2)
fig.subplots_adjust(hspace=0.5, wspace=0.5)
for z in np.flip(np.linspace(0, self.box_height,3), 0):
xval = list()
yval = list()
for x in np.linspace(0, self.box_half_length, 500):
xval.append(x)
y = 0.0
vortex_velocity = self.compute_velocity( time,np.array([x,y,z]) )
yval.append(vortex_velocity[1])
axs[0,0].plot(xval, yval, label='z =' + str(round(z,2))+' [m]')
axs[0,0].set_xlabel('Radial distance [m]')
axs[0,0].set_ylabel('Tangential speed [m/s]')
axs[0,0].legend(loc='best')
fig.suptitle('Vortex Flow at '+str(round(time,1))+' [s]')
axs[0,0].grid(True)
xval = list()
yval = list()
for z in np.linspace(0,self.box_height,50):
yval.append(z)
vortex_velocity = self.compute_velocity( time,np.array([0.0,0.0,z]) )
xval.append(vortex_velocity[2])
axs[0,1].plot(xval,yval)
axs[0,1].set_xlabel('Vertical speed [m/s]')
axs[0,1].set_ylabel('Height [m]')
axs[0,1].grid(True)
fig.savefig('vortex_velocity.png',dpi=200)
plt.close(fig)
return
