Scribd Logo
Upload

Welcome to Scribd!

  • 39 views

    CFD Applications

    Uploaded by

    qwe123
    This document contains the solutions to 4 questions from a numerical methods homework assignment. Question 1 involves using the 4th order Runge-Kutta method to solve the Van der Pol equation for different time steps. Question 2 uses RK-4 to model the free fall of a steel sphere in air and calculates the drag coefficient. Question 3 is not solved. Question 4 applies the TDMA algorithm to solve a fin problem by defining a new temperature difference parameter. Diagonal dominance is also checked for the matrix coefficients.

    Copyright:

    © All Rights Reserved
    Download as PDF, TXT or read online from Scribd

    CFD Applications

    Uploaded by

    qwe123
    39 views18 pages
    This document contains the solutions to 4 questions from a numerical methods homework assignment. Question 1 involves using the 4th order Runge-Kutta method to solve the Van der Pol equation for different time steps. Question 2 uses RK-4 to model the free fall of a steel sphere in air and calculates the drag coefficient. Question 3 is not solved. Question 4 applies the TDMA algorithm to solve a fin problem by defining a new temperature difference parameter. Diagonal dominance is also checked for the matrix coefficients.

    Original Description:

    CFD Applications

    Original Title

    CFD applications

    Copyright

    © © All Rights Reserved

    Available Formats

    PDF, TXT or read online from Scribd

    Did you find this document useful?

    Is this content inappropriate?

    This document contains the solutions to 4 questions from a numerical methods homework assignment. Question 1 involves using the 4th order Runge-Kutta method to solve the Van der Pol equation for different time steps. Question 2 uses RK-4 to model the free fall of a steel sphere in air and calculates the drag coefficient. Question 3 is not solved. Question 4 applies the TDMA algorithm to solve a fin problem by defining a new temperature difference parameter. Diagonal dominance is also checked for the matrix coefficients.

    Copyright:

    © All Rights Reserved
    Download as PDF, TXT or read online from Scribd
    Download as pdf or txt
    39 views18 pages

    CFD Applications

    Uploaded by

    qwe123
    This document contains the solutions to 4 questions from a numerical methods homework assignment. Question 1 involves using the 4th order Runge-Kutta method to solve the Van der Pol equation for different time steps. Question 2 uses RK-4 to model the free fall of a steel sphere in air and calculates the drag coefficient. Question 3 is not solved. Question 4 applies the TDMA algorithm to solve a fin problem by defining a new temperature difference parameter. Diagonal dominance is also checked for the matrix coefficients.

    Copyright:

    © All Rights Reserved
    Download as PDF, TXT or read online from Scribd
    Download as pdf or txt
    of 18

    NUMERICAL METHODS IN FLUID FLOW AND

    HEAT TRANSFER

    MIA 502E

    SPRING 2017

    Homework 1

    Name and Surname: Alperen YILDIZELİ


    Student ID: 503161105
    Question1
    Fourth order Runge-Kutta method is used for equation and new parameter (y) is defined as the
    derivative of x with respect to time. Solution is done for 3 time steps seperately. Thus, code is
    repeated itself 3 times. Definitions are:

    RK-4 method is not suitable for time steps given, Divergent results for 100 of dynamic viscosity :

    Graphic 1: Divergent solutions for given step sizes


    In order to obtain meaningful results time steps are divided by 1000. This change gives better
    results:

    Graphic 2: Convergent solutions for reorganized step sizes


    Question2
    In order solve 2nd Question, Newton’s 2nd law of motion is applied on steel sphere, Then
    velocity and displacement of sphere solved numerically as a set of linear equations with 4th
    order Runge-Kutta method.

    Problem is explained below with free body diagram and force equlibrium equation for x
    direction:

    For determinig Cd coefficient, given colerations are used:

    For Re <<1:
    Cd=24/Re

    For 1<Re<<400:
    Cd=24/(Re^0.646)

    For higher Reynolds Numbers:

    Re=0.5
    Results are given below:

    For h=0.1:

    Graphic 3: Solutions for d=0.001m

    Graphic 4: Solutions for d=0.01m


    Graphic 5: Solutions for d=0.02m

    Graphic 6: Solutions for d=0.07m


    For h=0.05:

    Graphic 7: Solutions for d=0.001m

    Graphic 8: Solutions for d=0.01m


    Graphic 9: Solutions for d=0.02m

    Graphic 10: Solutions for d=0.07m


    Also drag coefficient have been calculated and plotted:

    Graphic 11: Cd vs Reynolds Number

    Question3 didnt solved.


    Question4
    In order to solve fin problem crossectional area and the surrounding area is defined and
    TDMA alogrithm applied to problem. In order to homogenize the equation, new parameter
    “Teta” is defined as temperature diffrence between fin and ambient temperature.

    Crosssectional Area and Surrounding Area:

    Before applying TDMA alogrithm. Diagonal domiance condition controlled if it is satisfied:


    Temperature distribution of fin:
    Codes:
    Question 1:
    clc
    clear

    %% Explanation of Code

    % This code aims to solve non-linear Van-der-Pool equation


    numericaly
    % As numerical method, RK-4 is used and new parameter (y) is
    defined
    % as the derivative of x with respect to time. Solution is done
    for 3 time
    % steps seperately. Thus, code is repeated itself 3 times.

    %% Problem Data
    m1 = 100; %Dynamic viscosity [kg/ms]
    m2 = 100; %Dynamic viscosity [kg/ms]
    m3 = 100; %Dynamic viscosity [kg/ms]

    h1 = 0.0001 ; %Time Step value [s]


    h2 = 0.001 ; %Time Step value [s]
    h3 = 0.001 ; %Time Step value [s]

    t1 =[0:h1:500] ; %Time domain [s]


    t2 =[0:h2:500] ; %Time domain [s]
    t3 =[0:h3:500] ; %Time domain [s]

    %% Initial Values And Arrays for x and y

    x1=zeros(1,length(t1));
    x1(1)=2 ;
    x2=zeros(1,length(t2));
    x2(1)=2 ;
    x3=zeros(1,length(t3));
    x3(1)=2 ;

    y1=zeros(1,length(t1));
    y1(1)=0 ;
    y2=zeros(1,length(t2));
    y2(1)=0 ;
    y3=zeros(1,length(t3));
    y3(1)=0 ;

    %% Numerical Solution

    %% For h=0.0001
    for i=1:1:(length(t1)-1)

    f1 = ( m1 * (1-(x1(i)^2)) * y1(i) ) - x1(i) ;


    f2 = ( m1 * (1-((x1(i)+(h1/2))^2)) * (y1(i)+(h1*f1/2)) ) -
    x1(i) ;
    f3 = ( m1 * (1-((x1(i)+(h1/2))^2)) * (y1(i)+(h1*f2/2)) ) -
    x1(i) ;
    f4 = ( m1 * (1-((x1(i)+(h1))^2)) * (y1(i)+(h1*f3)) ) -
    x1(i) ;

    y1(i+1) = y1(i) + ((h1/6) * (f1 + (2*f2) + (2*f3) + f4 )) ;

    g1 = y1(i) ;
    g2 = y1(i)+ (h1*g1/2) ;
    g3 = y1(i)+ (h1*g2/2) ;
    g4 = y1(i)+ (h1*g3) ;

    x1(i+1) = x1(i) + ((h1/6) * (g1 + (2*g2) + (2*g3) + g4 ))


    ;
    end

    %% For h=0.001
    for i=1:1:(length(t2)-1)

    f1 = ( m2 * (1-(x2(i)^2)) * y2(i) ) - x2(i) ;


    f2 = ( m2 * (1-((x2(i)+(h2/2))^2)) * (y2(i)+(h2*f1/2)) ) -
    x2(i) ;
    f3 = ( m2 * (1-((x2(i)+(h2/2))^2)) * (y2(i)+(h2*f2/2)) ) -
    x2(i) ;
    f4 = ( m2 * (1-((x2(i)+(h2))^2)) * (y2(i)+(h2*f3)) ) -
    x2(i) ;

    y2(i+1) = y2(i) + ((h2/6) * (f1 + (2*f2) + (2*f3) + f4 )) ;

    g1 = y2(i) ;
    g2 = y2(i)+ (h2*g1/2) ;
    g3 = y2(i)+ (h2*g2/2) ;
    g4 = y2(i)+ (h2*g3) ;

    x2(i+1) = x2(i) + ((h2/6) * (g1 + (2*g2) + (2*g3) + g4 ))


    ;
    end

    %% For h=0.005
    for i=1:1:(length(t3)-1)

    f1 = ( m3 * (1-(x3(i)^2)) * y3(i) ) - x3(i) ;


    f2 = ( m3 * (1-((x3(i)+(h3/2))^2)) * (y3(i)+(h3*f1/2)) ) -
    x3(i) ;
    f3 = ( m3 * (1-((x3(i)+(h3/2))^2)) * (y3(i)+(h3*f2/2)) ) -
    x3(i) ;
    f4 = ( m3 * (1-((x3(i)+(h3))^2)) * (y3(i)+(h3*f3)) ) -
    x3(i) ;
    y3(i+1) = y3(i) + ((h3/6) * (f1 + (2*f2) + (2*f3) + f4 )) ;

    g1 = y3(i) ;
    g2 = y3(i)+ (h3*g1/2) ;
    g3 = y3(i)+ (h3*g2/2) ;
    g4 = y3(i)+ (h3*g3) ;

    x3(i+1) = x3(i) + ((h3/6) * (g1 + (2*g2) + (2*g3) + g4 ))


    ;
    end

    subplot(3,2,1)
    plot(t1,x1)
    axis([0 500 -2.5 +2.5])
    title('h=0.0001s')
    xlabel('t [s]')
    ylabel('x [m]')

    subplot(3,2,2)
    plot(x1,y1)
    axis([-2.5 2.5 -150 +150])
    title('h=0.0001s')
    xlabel('x [m]')
    ylabel('y [m]')

    subplot(3,2,3)
    plot(t2,x2)
    axis([0 500 -2.5 +2.5])
    title('h=0.001s')
    xlabel('t [s]')
    ylabel('x [m]')

    subplot(3,2,4)
    plot(x2,y2)
    axis([-2.5 2.5 -150 +150])
    title('h=0.001s')
    xlabel('x [m]')
    ylabel('y [m]')

    subplot(3,2,5)
    plot(t3,x3)
    axis([0 500 -2.5 +2.5])
    title('h=0.005s')
    xlabel('t [s]')
    ylabel('x [m]')

    subplot(3,2,6)
    plot(x3,y3)
    axis([-2.5 2.5 -150 +150])
    title('h=0.005s')
    xlabel('x [m]')
    ylabel('y [m]')

    Question 2:
    clc
    clear

    %% Explanation of Code

    % This code aims to solve free falling of a steel sphere in air


    % As numerical method, RK-4 is used. Since velocity is
    derivative of
    % displacepent wrt time it can be called set of ODE.

    %% Problem Data
    rho_fluid = 1.29 ; % Density of fluid [kg/m^3]
    rho_steel = 8000 ; % Density of steel [kg/m^3]
    nu = 1.49*(10^-5) ; % Kinematic viscosity [m^2/s]
    g = 9.81 ; % Gravitational accelaration [kg/m^3]

    d=0.001; % Diameter of steel sphere [m]


    V=4*pi*(d^3)/3 ; % Volume of steel sphere [m^3]
    Area=pi*d^2/4; % Facing area (Crosssectional Area) [m^2]

    m_fluid = V*rho_fluid ; % Mass of fluid [kg]


    m_steel = V*rho_steel ; % Mass of steel sphere [kg]

    h = 0.05 ; %Time Step value [s]

    t =[0:h:10] ; %Time domain [s]

    %% Initial Values And Arrays for x and y

    x=zeros(1,length(t));
    v=zeros(1,length(t));

    %% Numerical Solution

    A=m_steel+(m_fluid/2);
    B=g*(m_steel-m_fluid);

    for i=1:1:(length(t)-1)

    Re(i)=rho_fluid*v(i)*d/nu;
    Cd2(i)=24/Re(i);
    Cd3(i)=24/(Re(i)^0.646);
    Cd4(i)=0.5;
    if Re < 1
    Cd(i)=24/Re(i);
    elseif Re < 400
    Cd(i)=24/(Re(i)^0.646);
    else
    Cd(i)=0.5 ;
    end

    C=0.5*rho_fluid*Area ; %(v(i)^2)

    f1 = (B- (C * ( v(i) ^ 2 ) ) )/A ;


    f2 = (B- (C * ( (v(i)+(h*f1/2)) ^ 2 ) ) )/A ;
    f3 = (B- (C * ( (v(i)+(h*f2/2)) ^ 2 ) ) )/A ;
    f4 = (B- (C * ( (v(i)+(h*f3)) ^ 2 ) ) )/A ;

    v(i+1) = v(i) + ((h/6) * (f1 + (2*f2) + (2*f3) + f4 )) ;

    g1 = v(i) ;
    g2 = v(i)+ (h*g1/2) ;
    g3 = v(i)+ (h*g2/2) ;
    g4 = v(i)+ (h*g3) ;

    x(i+1) = x(i) + ((h/6) * (g1 + (2*g2) + (2*g3) + g4 )) ;


    end

    figure(1)
    subplot(1,2,1)
    plot(t,x)

    title('Displacement of Sphere for diameter of d=0.001m')


    xlabel('t [s]')
    ylabel('x [s]')

    subplot(1,2,2)
    plot(t,v)
    title('Velocity of Sphere for diameter of d d=0.001m')
    xlabel('t [s]')
    ylabel('v [m]')

    for i=length(Re)

    end

    figure(2)
    hold on
    semilogy(Re,Cd, 'LineWidth',1)
    semilogy(Re,Cd2,'x')
    semilogy(Re,Cd3,'+')
    semilogy(Re,Cd4,'g')
    legend('Cd','24/Re','(24/(Re^0.646))','0.5')

    figure(3)
    hold on
    semilogy(Re,Cd, 'LineWidth',1)
    xlabel('Re')
    ylabel('Cd')
    legend('Cd')

    Question 4:
    clc
    clear

    %% Explanation of Code

    % This code aims to solve fin problem numericaly


    % As numerical method, TDMA is used and new parameter (teta) is
    defined
    % as the T-Tinf.

    %% Problem Data

    %% Initial Values And Arrays for x and y

    y=[0:h:l]; %Spatial domain


    m=length(y);

    T=zeros(1,m);
    T(1)=473;
    Teta=zeros(1,m);
    Teta(1)=T(1)-293;
    %% Numerical Solution

    %Calculation of Area values and their derivatives with respect


    to y
    for i=1:1:m
    Ac(i) = L*(tb - (y(i)*((tb-tt)/(l)))) ;
    As(i) = 2*(L + tt - (y(i)*((tt-tb)/l)));

    dAs(i) = (2*L)-(2*((tt-tb)/(l))*y(i)) ;
    dAc(i) = (-1*L*(tb-tt))/(l);
    end
    %Matrix Coefficients
    for i=1:1:m-2
    a(i) = -2 - (h^2)*(hc*dAs(i)/(k*Ac(i)));
    b(i) = 1 + (h*dAc(i)/(2*Ac(i)));
    c(i) = 1 - (h*dAc(i)/(2*Ac(i)));
    d(i)=0 ;
    Domiance_Control(i)=abs(a(i))-abs(b(i)+c(i)) ;
    end

    d(1)=-c(1)*(Teta(1));
    a(m-2)=a(m-2)+(b(m-2)/(1+(h*hc/k)));

    alfa(1)=a(1);
    gama(1)=d(1);

    for i=2:1:m-2
    alfa(i)= a(i)-b(i-1)*c(i)/alfa(i-1);
    gama(i)=d(i)-gama(i-1)*c(i)/alfa(i-1);
    end

    Teta(10)=gama(9)/alfa(9);
    Teta(11)=Teta(10)/(1+(h*hc/k));

    for i=m-2:-1:2
    Teta(i)=(gama(i-1)-b(i-1)*Teta(i+1))/alfa(i-1);
    end

    for i=1:1:m

    T(i)=Teta(i)+293;

    end

    figure(1)
    plot(Domiance_Control)
    ylabel('|a|-|b+c|')
    title('Diagonal Domiance Control')

    Refferences:

    Numerical Methods in Fluid Flow and Heat Transfer Lecrure Notes

    You might also like