Mathematica Code for streamfunction--vorticity CFDxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxstartinglength Xsyxinflowoutflowsurface lengthL = 1Hi,0i,1i,NyFigure 1: External boundary layer control volume.4 Mathematica Code for streamfunction–vorticity CFDListed below is a Mathematica co de which implements the solution procedure for the flat plate, externalboundary layer problem. The domain consists of a rectangular region show in Fig. (1). The mathematicalformulation of the boundary conditions for this flow configuration was discussed in the previous set of notes.The objective here is to show how it is coded into Mathematica.The code consists of the following basic modules (or elements):• A module to numerically solve for stream function values, as a function of the vorticity values and thestream function boundary conditions. The system of equations that are solved have the formψi+1,j+ ψi−1,j− 2ψi,j(∆x)2+ψi,j+1+ ψi,j−1− 2ψi,j(∆y)2= −ωi,j(1)• A module to explicitly calculate the values of either ωi,jor Ti,jat the k time step from the flowinformation at the k − 1 time step. This procedure utilizes the explicit time integration, upwind–differencing scheme.• Associated modules and functions to interpolate the numerical results, compute derivatives of results,produce plots, and so on.The code is presented as follows• The first lines set up global definitions and define assorted utility functions.Remove["Global‘*"]Off[General::"spell1",General::"spell"]sp[x_]:=Simplify[PowerExpand[x]]$TextStyle={FontFamily->"Times",FontSize->16};zeroprint[n_]:=Module[{s},s="";Do[s=s<>"0",{i,1,n}];s];index[l_]:=Which[l<10,"00",l<100,"0",1==1,""]<>ToString[l];1printtodigits[x_,nd_]:=Module[{t,lt,ndot,ndrop,t2,nzero,t3},t=ToString[x];lt=SyntaxLength[t];ndot=StringPosition[t,"."][[1,1]];ndrop=Max[0,lt-ndot-nd];t2=StringDrop[t,-ndrop];nzero=Max[0,nd+ndot-lt];t3=t2<>zeroprint[nzero];t3]• The system consists of a rectangular region, with length 1 + Xs, where Xsis the starting region, and aheight of H. The number of interior control volumes in the x and y directions is nx-1 and ny-1. Thesizes of the control volume, delx and dely, are∆x =1 + XsNx, ∆y =HNy• The module (or subroutine) psisoln[omega] generates the numerical solution of the stream functionDE for ψ as a function of ω. The vorticity ω is stored in the variable omega[i,j], with i = 0, 1, . . . , Nxand j = 0, 1, . . . , Ny. The module stores the calculated values of ψ in psi[i,j].psisoln[omega_] := Module[{vars, eqns, psit, soln},vars = Flatten[Table[psit[i, j], {i, 1, nx - 1}, {j, 1, ny - 1}]];eqns =Flatten[Table[(psit[i - 1, j] + psit[i + 1, j] - 2 psit[i, j])/(delx)^2+ (psit[i, j - 1] + psit[i, j + 1] - 2 psit[i, j])/(dely)^2== -omega[i, j], {i, 1, nx - 1}, {j, 1, ny - 1}]];Do[psit[0, j] = j dely;psit[nx, j] = 2 psit[nx - 1, j] - psit[nx - 2, j];, {j, 0, ny}];Do[psit[i, 0] = 0;psit[i, ny] = (dely + 2 psit[i, ny - 1] - 1/2 psit[i, ny - 2]) 2/3;, {i, 0, nx}];soln = NSolve[eqns, vars];Do[psi[i, j] = psit[i, j] /. soln[[1]], {i, 0, nx}, {j , 0, ny}];]The first line defines the local variables: vars, eqns, psit, soln. Note the symbolic set–up of thenumerical procedure; vars represents the variables to be solved, and eqns are the symbolic forms of theequations. The command Flatten turns a 2–D list – which would be produced by the double iterator{i,1,nx-1},{j,1,ny-1} on the Table commands – into a 1–D list, i.e., it collapses the matrix intoa vector. The two Do loops impose the boundary conditions, and the command NSolve performs anumerical solution of the equations. Finally, the solution is transferred into the global variable psi.2• The module deltmax calculates the largest allowable time step size, as a function of the Reynoldsnumber Re and the current values of the stream function.deltmax[re_]:=Module[{ull,uc,urr,vbb,vc,vtt,t1,tmax,ul,ur,vt,vb},tmax=1000;Do[ull=(psi[i-1,j+1]-psi[i-1,j-1])/(2 dely);uc=(psi[i,j+1]-psi[i,j-1])/(2 dely);urr=(psi[i+1,j+1]-psi[i+1,j-1])/(2 dely);vbb=-(psi[i+1,j-1]-psi[i-1,j-1])/(2 delx);vc=-(psi[i+1,j]-psi[i-1,j])/(2 delx);vtt=-(psi[i+1,j+1]-psi[i-1,j+1])/(2 delx);ul=(ull+uc)/2;ur=(urr+uc)/2;vb=(vbb+vc)/2;vt=(vtt+vc)/2;t1=1/(2/re/delx^2+2/re/dely^2+(ur+Abs[ur]-ul+Abs[ul])/2/delx+(vt+Abs[vt]-vb+Abs[vb])/2/dely);tmax=Min[tmax,t1];,{i,1,nx-1},{j,1,ny-1}];tmax]• The module explicit applies the explicit time integration formula to calculate the updated value ofthe transported quantity (ω or T ) at interior node i, j. Input variables are (in the order that they arelisted): the node index, Reynolds number, ∆t, stream function, the transported quantity (ω in thiscase), and the source array.explicit[i_,j_,re_,deltime_,psi_,omega_,source_]:=Module[{ull,uc,urr,vbb,vc,vtt,a,b,c,d,e,ul,ur,vt,vb},ull=(psi[i-1,j+1]-psi[i-1,j-1])/(2 dely);uc=(psi[i,j+1]-psi[i,j-1])/(2 dely);urr=(psi[i+1,j+1]-psi[i+1,j-1])/(2 dely);vbb=-(psi[i+1,j-1]-psi[i-1,j-1])/(2 delx);vc=-(psi[i+1,j]-psi[i-1,j])/(2 delx);vtt=-(psi[i+1,j+1]-psi[i-1,j+1])/(2 delx);ul=(ull+uc)/2;ur=(urr+uc)/2;vb=(vbb+vc)/2;vt=(vtt+vc)/2;a=deltime(1/re/delx^2-(ur-Abs[ur])/2/delx);b=deltime(1/re/delx^2+(ul+Abs[ul])/2/delx);c=deltime(1/re/dely^2-(vt-Abs[vt])/2/dely);d=deltime(1/re/dely^2+(vb+Abs[vb])/2/dely);e=1-deltime(2/re/delx^2+2/re/dely^2+(ur+Abs[ur]-ul+Abs[ul])/2/delx+(vt+Abs[vt]-vb+Abs[vb])/2/dely);a omega[i+1,j]+b omega[i-1,j]+c omega[i,j+1]+d omega[i,j-1]+ e omega[i,j]+source[i,j] deltime]• The inttable module is used to set up an interpolation table from a calculated data array.inttable[dat_]:=Module[{},Flatten[Table[{N[i],N[j],dat[i,j]},{i,0,nx},{j,0,ny}],1]]3• The parameters are now defined. A 40×20 mesh is used, and ReL= 104, P r = 1.0. The startinglength is 0.1, and the height H is taken to be around twice the estimated boundary layer thickness atthe