!************************************************************************ !* Educationally-Designed Unstructured 2D (EDU2D) Code !* !* --- EDU2D-AirfoilSpline !* !* This program reads the OM6 wing section data (discrete points), !* and constructs cubic splines for the wing section. The section data !* contains only the upper surface of the airfoil. Spline will be !* constructed for the upper surface only. !* !------------------------------------------------------------------------ !* !* Input: !* !* 'om6_wing_section_sharp.dat': a set of discrete points, which defines !* the upper surface of an airfoil. !* !* !* Output: !* !* 'fort.100' : input data (in a slightly different format) !* 'fort.200' : 500 sampled points over the spline (smooth curve). !* 'edu2d_airfoil_input_tec.dat' : Input dat (for Tecplot) !* 'edu2d_airfoil_spline_tec.dat': 500 sampled points (for Tecplot) !* !------------------------------------------------------------------------ !* !* Note: Spline gives a picewise cubic polynomial representation of the !* input data. I used it a long time ago for grid adaptation for !* a Joukowsky airfoil, to move boundary nodes along the airfoil. !* !* Note: Matlab file 'edu2d_airfoil_plot.m' load the data from fort.100 !* and fort.200, and plot the airfoil. The matlab file will plot !* a complete airfoil by copying the upper surface to the lower !* surface. So, it shows a symmetric airfoil. !* !* !* This is Version 1 (February 08, 2017). !* !* This F90 code is written and made available for an educational purpose. !* This file may be updated in future. !* !* Katate Masatsuka, http://www.cfdbooks.com !************************************************************************ program edu2d_airfoil_spline implicit none integer , parameter :: dp = selected_real_kind(p=15) real(dp), parameter :: zero = 0.0_dp real(dp), dimension( :), allocatable :: t, x, y, xb, yb real(dp), dimension(:,:), allocatable :: ax, ay real(dp) :: tp, sf, ts integer :: n, nb, i, os !************************************************************************ !* Read a set of discrete data: x(i) and y(i), i = 1, n, which represent !* the upper surface of an OM6 airfoil section with a TE closure. !* !* Note: !* !* The data file "om6_wing_section_sharp.dat" is based on Table 3 in !* AIAA Journal, Vol. 54, No. 9, September 2016, !* "Reynolds-Averaged Navier£üStokes Simulations on NACA0012 and ONERA-M6 !* Wing with the ONERA elsA Solver" by Mayeur, A. Dumont, D. Destarac, !* and V. Gleize. !* !* Note: TE closure is also as described in the above paper. !* !************************************************************************ write(*,*) write(*,*) " Reading the data file: om6_wing_section_sharp.dat" write(*,*) " - Discrete points of OM6 wing section....." write(*,*) open(unit=1, file='om6_wing_section_sharp.dat', status="unknown", iostat=os) !Read the number of points: read(1,*) n allocate( ax(4,n), ay(4,n) ) allocate( t(n), x(n), y(n)) !Read x(i) do i = 1, n read(1,*) x(i) end do !Read y(i) do i = 1, n read(1,*) y(i) end do close(1) write(*,*) write(*,*) " --- Number of points read = ", n write(*,*) " --- This is the upper surface of a symmetric airfoil. " write(*,*) !************************************************************************ !* Construct the parametrization of (x(t),y(t)): t=[0,1]. !* The subroutine creates the parameter array t(1:n). !* !* Input: x(1:n), y(1:n), n; Output: t(1:n) !************************************************************************ write(*,*) write(*,*) " Parametrizing the data: x(t), y(t): " write(*,*) call param(x,y,t,n) write(*,*) " t(1) = ", t(1) write(*,*) " t(n) = ", t(n) write(*,*) !************************************************************************ !* Construct cubic splines: !* !* -> Construct spline matrices ax(1:4,1:n-1) and ay(1:4,1:n-1), !* for ( x(1:n),t(1:n) ) and ( y(1:n),t(1:n) ), respectively. !* !* Note: Zero slopes at LE, and zero second derivatives at TE. !* !* Note: ax(:,n) and ay(:,n) are used only for a temporary computation, !* and not represent a cubic polynomials. A compelte set of cubic !* poynomials are represented by ax(:,1:n-1) and ay(:,1:n-1). !************************************************************************ write(*,*) write(*,*) " Construct a cubic spline matrix for x(t)... " write(*,*) call construct_cubic_spline(x, t,n, 1,0, zero,zero, ax) write(*,*) write(*,*) " Construct a cubic spline matrix for y(t)... " write(*,*) call construct_cubic_spline(y, t,n, 1,0, zero,zero, ay) !************************************************************************ !* Just for the sake of visualization, sample a lot of points in the !* spline (to visualize the spline curve). !* !* Sample points are saved in xb and yb, i = 1, nb. !* These points will show how the spline looks like. !* !************************************************************************ !Sample 500 points in the spline. nb = 500 allocate( xb(nb), yb(nb) ) write(*,*) write(*,*) " Sampling points over the spline curve... nb = ", nb write(*,*) xb(1) = x(1) yb(1) = y(1) do i = 2, nb-1 tp = real(i-1,dp)/real(nb-1,dp) sf = 8.0_dp; ts = ( tanh(sf*(tp-0.5_dp)) - tanh(sf*(-0.5_dp)) ) & /( tanh(sf*( 0.5_dp)) - tanh(sf*(-0.5_dp)) ) xb(i) = cubic_spline_geometry(ts, t,ax,n) yb(i) = cubic_spline_geometry(ts, t,ay,n) end do xb(nb) = x(n) yb(nb) = y(n) !************************************************************************ ! Write out the original discerte data, and the interpolated data. !************************************************************************ !----------------------------------------------------------------- ! Data points (input discrete data) write(*,*) write(*,*) " Writing out the input discrete points in the file, fort.100 ..." write(*,*) " --- To be plotted by Matlab (using edu2d_airfoil_plot.m)" write(*,*) " --- Number of points = ", n write(*,*) do i = 1, n write(100,*) x(i) ,y(i) end do !----------------------------------------------------------------- ! Data points (input discrete data) in a Tecplot file. write(*,*) write(*,*) " Writing out the input discrete points in Tecplot file:" write(*,*) " --- edu2d_airfoil_input_tec.dat" write(*,*) " --- Number of points = ", n write(*,*) open(unit=10, file="edu2d_airfoil_input_tec.dat", status="unknown", iostat=os) write(10,*) 'TITLE = "2D airfoil input points (upper surface only)"' write(10,*) 'VARIABLES = "x","y"' write(10,*) 'zone i =', n, ' f=point' do i = 1, n write(10,'(2es25.15)') x(i) ,y(i) end do close(10) !----------------------------------------------------------------- ! Interpolated points (representing the spline; this is smooth). write(*,*) write(*,*) " Writing sampled points of spline in the file, fort.200 ..." write(*,*) " --- To be plotted by Matlab (using edu2d_airfoil_plot.m)" write(*,*) " --- Number of points = ", nb write(*,*) do i = 1, nb write(200,*) xb(i) ,yb(i) end do !----------------------------------------------------------------- ! Interpolated points (representing the spline; this is smooth) in a Tecplot file. write(*,*) write(*,*) " Writing sampled points of spline in Tecplot file:" write(*,*) " --- edu2d_airfoil_spline_tec.dat" write(*,*) " --- Number of points = ", nb write(*,*) open(unit=11, file="edu2d_airfoil_spline_tec.dat", status="unknown", iostat=os) write(11,*) 'TITLE = "2D airfoil with spline (upper surface only)"' write(11,*) 'VARIABLES = "x","y"' write(11,*) 'zone i =', nb, ' f=point' do i = 1, nb write(11,'(2es25.15)') xb(i) ,yb(i) end do close(11) !----------------------------------------------------------------- write(*,*) write(*,*) " Successfully finished the cubic spline..." write(*,*) " Open the tecplot files to view the airfoil" write(*,*) " or use edu2d_airfoil_plot.m to plot the results in Matlab." write(*,*) stop !************************************************************************ contains !************************************************************************ !* This program constructs a parametrization of a curve in (x,y) plane !* using arc-length as a parameter. !* !* Input : x, y : coordinates of the defining points. !* !* Output: t : parameter values of x(t) and y(t). !* !* !* Note: (x(i),y(i)) = ( x(t(i)), y(t(i)) ), i = 1, n. !* !* Hiroaki Nishikawa, October 29, 1998 !************************************************************************ subroutine param(x, y, t, n) implicit none integer , parameter :: dp = selected_real_kind(p=15) !input integer , intent(in) :: n real(dp), dimension(n), intent(in) :: x, y !Output real(dp), dimension(n), intent(out) :: t !Local variables integer :: i !************************************************************************ !* Compute the arc-length parameter. !************************************************************************ t(1) = zero do i = 2, n t(i) = t(i-1) + sqrt( (x(i)-x(i-1))**2 + (y(i)-y(i-1))**2 ) end do !************************************************************************ !* Normalize the parameter such that t = [0,1]. !************************************************************************ do i = 2, n t(i) = t(i)/t(n) end do return end subroutine param !************************************************************************ !************************************************************************ !************************************************************************ !* This program constructs a cubic spline of a function f(t), t=[0,1]. !* !* Input: n: number of points !* df_0, df_n: derivatives at the end points (t=0,1) !* in case of clamped splines. !* !* Output: a: 4x(n-1) spline matrix. !* !* !* Note: A cubic polynomial over [i-1, i] is defined by a(:,i-1). !* !* o------o------o----o----o---------o---o------o !* i=1 i-1 i i=n !* a(:,i-1) !* !* Hiroaki nishikawa, October 29, 1998 !************************************************************************ subroutine construct_cubic_spline(f,t,n,cnd0,cnd1,df_0,df_n,a) implicit none integer , parameter :: dp = selected_real_kind(p=15) !Input integer , intent(in) :: n integer , intent(in) :: cnd0, cnd1 real(dp), dimension(n), intent(in) :: f, t real(dp), intent(in) :: df_0, df_n !Output real(dp), dimension(4,n), intent(out) :: a !Local variables real(dp), dimension(n) :: h, q, u, z, l integer :: i, j !************************************************************************ !* a(1,:) are just the function values. !************************************************************************ do i = 1, n a(1,i) = f(i) end do !************************************************************************ !* Construct the right hand side. !************************************************************************ do i = 1, n-1 h(i) = t(i+1)-t(i) end do q = zero !Clamped spline (specify the derivative at the left end): if (cnd0 == 1) q(1) = 3.0d0*(a(1,1)-a(1,1))/h(1)-3.0d0*df_0 !Clamped spline (specify the derivative at the right end): if (cnd1 == 1) q(n) = 3.0d0*df_n-3.0d0*(a(1,n)-a(1,n-1))/h(n-1) do i = 2, n-1 q(i) = 3.0d0*(a(1,i+1)-a(1,i))/h(i)-3.0d0*(a(1,i)-a(1,i-1))/h(i-1) end do !************************************************************************ !* Solve the linear system for a(3,:). !************************************************************************ !f''=0 (natural spline) if (cnd0 == 0) then l(1) = 1.0d0 u(1) = 0.0d0 z(1) = 0.0d0 !f'=df_n (clamped spline) elseif (cnd0 == 1) then l(1) = 2.0d0*h(1) u(1) = 0.5d0 z(1) = q(1)/l(1) endif do i = 2, n-1 l(i) = 2.0d0*(t(i+1)-t(i-1))-h(i-1)*u(i-1) u(i) = h(i)/l(i) z(i) = (q(i)-h(i-1)*z(i-1))/l(i) end do !f''=0 (natural spline) if (cnd0 == 0) then l(n) = 1.0d0 z(n) = 0.0d0 a(3,n) = 0.0d0 !f'=df_n (clamped spline) elseif (cnd1 == 1) then l(n) = h(n-1)*(2.0-u(n-1)) z(n) = (q(n)-h(n-1)*z(n-1))/l(n) a(3,n) = z(n) endif !--------------------------------------------------------- ! Construct the coefficients of cubic polynomials: ! a(:,i-1) over [i-1, i], i = 1, n-1. do j = n-1, 1, -1 !--------------------------------------------------------- ! Back substitution to compute a(3,:) a(3,j) = z(j) - u(j)*a(3,j+1) !--------------------------------------------------------- ! Compute bj and dj from cj: ! ! cubic_j = aj + bj*t + cj*t^2 + dj*t^3 ! = a(1,j) + a(2,j)*t + a(3,j)*t^2 + a(4,j)*t^3 a(2,j)=(a(1,j+1)-a(1,j))/h(j) - h(j)*(a(3,j+1)+2.0d0*a(3,j))/3.0d0 a(4,j)=(a(3,j+1)-a(3,j))/(3.0d0*h(j)) end do !--------------------------------------------------------- !--------------------------------------------------------- return end subroutine construct_cubic_spline !*********************************************************************** !*********************************************************************** !************************************************************************ !* This program computes the coords of a point of a curve for a given !* parameter value. !* !* Input: tp: parameter at which (x,y) are sought !* a: spline matrix !* t: defining points in the parameter space !* !* Output: fp: the coordinate value at t=tp. !* !* !* Hiroaki Nishikawa, October 29, 1998 !************************************************************************ function cubic_spline_geometry(tp,t,a,n) result(fp) implicit none integer , parameter :: dp = selected_real_kind(p=15) !input integer , intent(in) :: n real(dp), dimension(n), intent(in) :: t real(dp), dimension(4,n), intent(in) :: a real(dp), intent(in) :: tp !Output real(dp) :: fp !Local variables integer :: i !************************************************************************ !* Find out in which interval tp lies. !* !* Example: !* !* o------o------o----o----o-----*---o---o------o !* i=1 2 3 4 5 tp 6 7 i=n=8 !* !* In this case, it finds i = 6, and uses the polynomial whose !* coefficients are defined by a(:,5). !* !************************************************************************ i = 1 do while (t(i) < tp) i = i + 1 if (i > n) then write(*,*) " Requested paramete is out of range... " write(*,*) " tp = ", tp, " > 1.0... Stop" stop endif end do if (i == 1) then i = 2 endif if ( tp < t(1) ) then write(*,*) " Requested paramete is out of range... " write(*,*) " t(1) = ", t(1) write(*,*) " tp = ", tp, " < t(1)... Stop" stop endif !************************************************************************ !* t is in between i and (i-1). !* then compute x by using the coefficeints for the cubic polynomial !* defined over [i-1, i], which are given by a(:,i-1). !************************************************************************ fp = a(1,i-1) & + a(2,i-1)*(tp-t(i-1)) & + a(3,i-1)*(tp-t(i-1))*(tp-t(i-1)) & + a(4,i-1)*(tp-t(i-1))*(tp-t(i-1))*(tp-t(i-1)) return end function cubic_spline_geometry !*********************************************************************** !*********************************************************************** end program edu2d_airfoil_spline