This tutorial is intended for the full version of the toolbox.
Problem Definition
In this tutorial we are going to study an environmental aerodynamics case. We will simulate turbulent air flow with an algebraic turbulence model designed and calibrated for outdoor aerodynamics problems. For details of the turbulence model see reference [1]. We are going to start with a simple quasi two-dimensional flow, where two long buildings standing parallel to each other are located perpendicularly to the direction of incoming wind. Hence, we only consider a 2-D flow case in the cross section plane and study, how the flow recirculates between the buildings and on the leeward side of the downstream located building. This simple model can seem to be oversimplified, however it is a good approximation of flow patterns in streets surrounded by two long facades of buildings on each side of the street and pavement. Similar situation is present on motorways along which wind or noise shields have been installed. In this tutorial we only study the flow field, but Tutorial 10 is focused on modelling contamitant transport in urban flows, such as ventilation of exhaust gases from the street level.
Geometry
Both geometry and mesh have already been prepared for this tutorial. To obtain good resolution of the flow solution in the turbulent boundary layer close to building walls, where big velocity gradients are present, additional mesh layers for boundary layers will be generated in QuickerSim CFD Toolbox. To make it clear: The initial mesh created in Gmsh is a pure unstructured, triangular mesh. The boundary layer mesh for 2-D problems can automatically be generated in this toolbox.
Code
Clear MATLAB Command Window and Workspace
clc; clear; % Read mesh [p,e,t] = importMeshGmsh('urbanMesh.msh'); % Add boundary layer mesh layers.h1 = 0.02; layers.q = 1.5; layers.n = 4; wallIds = [18 20 19 21 22 23]; [p,e,t] = extrudeLayers2D(p,e,t,wallIds,layers,40,0.8); % Convert mesh to second order and display mesh [p,e,t,nVnodes,nPnodes,indices] = convertMeshToSecondOrder(p,e,t); displayMesh2D(p,t);
% Define fluid nu = 1.5e-5; % Init solution and define convergence and under-relaxation parameters [u, convergence] = initSolution(p,t,[0 0],0); maxres = 1e-3; maxiter = 20; alpha = 0.6; % Compute wall distance and duplicate velocity vector d = computeWallDistance(p,e,t,wallIds); us = u; % Define inlet velocity profile vel = [0.43 0]; % Iterate nonlinear terms for iter = 1:maxiter % Compute turbulent viscosity nut = turbulentViscosity(p,t,u,d,nu,indices,'Srebric',0.1,5700); % Under-relax velocity and assemble matrix and right hand side u = alpha*u+(1-alpha)*us; % [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu+nut,u(indices.indu),u(indices.indv),'nosupg'); [NS, F] = assembleNavierStokesMatrix2D(p,e,t,nu+nut,u(indices.indu),u(indices.indv),'supgDoublyAsymptotic'); % Apply boundary conditions [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,15,'inlet', vel); [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,wallIds,'wall'); [NS, F] = imposeCfdBoundaryCondition2D(p,e,t,NS,F,16,'slipAlongX'); % Compute and plot residuals [stop, convergence] = computeResiduals(NS,F,u,size(p),convergence,maxres); plotResiduals(convergence,2); % Break if solution converged if(stop) break; end % Solve equations us = u; u = NS\F; end % Display Solution figure(5) displaySolution2D(p,t,u(indices.indu),'x-velocity');
figure(6) displaySolution2D(p,t,nut,'Turbulent viscosity [m^2/s^2]');
figure(7) displaySolution2D(p,t,nut/nu,'Turbulent viscosity ratio [-]');
References
[1] Y. Qian, J. Srebric, Development and Validation of an Algebraic Turbulence Model for Outdoor Airflow and Contaminant Simulations around Buildings, International Journal in Building, Urban, Interior and Landscape Technology, 1(1), 2010.