MATLAB Examples

2D Unsteady Navier-Stokes problem

Based on "Finite Element Methods for flow problems" of Jean Donea and Antonio Huerta

```% Andrea La Spina % https://it.linkedin.com/in/andrealaspina % https://independent.academia.edu/AndreaLaSpina % https://www.researchgate.net/profile/Andrea_La_Spina ```

Equation

```% u_t-v*div(grad(u))+(u.div)u+grad(p)=b in W % div(u)=0 in W % u=u_D on W_D ```

Initialization

```initialization ```

Input

```% Model parameters -------------------------------------------------------- v=1; % Kinematic viscosity b_fun=@(x,y) ... % Body force [(12-24*y)*x^4+(-24+48*y)*x^3+(-48*y+72*y^2-48*y^3+12)*x^2+... (-2+24*y-72*y^2+48*y^3)*x+(1-4*y+12*y^2-8*y^3),... (8-48*y+48*y^2)*x^3+(-12+72*y-72*y^2)*x^2+... (4-24*y+48*y^2-48*y^3+24*y^4)*x+(-12*y^2+24*y^3-12*y^4)]; % ------------------------------------------------------------------------- % Space ------------------------------------------------------------------- x_i=0; % Initial point (x) x_f=1; % Final point (x) y_i=0; % Initial point (y) y_f=1; % Final point (y) % ------------------------------------------------------------------------- % Time -------------------------------------------------------------------- t_i=0; % Initial time t_f=0.1; % Final time dt=t_f/10; % Time step animation_time=10; % Animation time theta=1; % Theta for time integration scheme % ( 0 = Forward Euler) % (1/2 = Crank-Nicolson) % (2/3 = Galerkin) % ( 1 = Backward Euler) % ------------------------------------------------------------------------- % FEM --------------------------------------------------------------------- n_el_x=10; % Number of finite elements (x) n_el_y=10; % Number of finite elements (y) n_gauss_side=3; % Number of Gauss points per side body_force_selection='gauss'; % Selection points of body force % ('centre','gauss','nodes_v') toll_p=1e-6; % Tollerance for pressure toll_v=1e-6; % Tollerance for velocity % ------------------------------------------------------------------------- % Boundary conditions ----------------------------------------------------- bound_cond_p=0; % Boundary conditions (pressure) bound_cond_v=0; % Boundary conditions (velocity) dof_constrained_string_p=... % DOFs constrained (pressure) ['[1]']; dof_constrained_string_v=... % DOFs constrained (velocity) ['[2:n_np_x_v-1,',... 'n_np_v-n_np_x_v+2:n_np_v-1,'... '1:n_np_x_v:n_np_v-n_np_x_v+1,'... 'n_np_x_v:n_np_x_v:n_np_v]']; % ------------------------------------------------------------------------- % Initial conditions ------------------------------------------------------ radius=0.2; % Radius of the hill x_0=1/2; % Centre of the hill (x) y_0=1/2; % Centre of the hill (y) u_0_fun=@(x,y) [0,0]; % Initial condition (velocity) % ------------------------------------------------------------------------- % Plots ------------------------------------------------------------------- plot_geometry='yes'; % Plot geometry ('yes','no') plot_shape_functions='yes'; % Plot shape functions ('yes','no') plot_test_functions='yes'; % Plot test functions ('yes','no') % ------------------------------------------------------------------------- ```

Derived parameters

```% Common parameters ------------------------------------------------------- n_sd=2; % Number of dimensional space L_x=x_f-x_i; % Domain length (x) L_y=y_f-y_i; % Domain length (y) n_el=n_el_x*n_el_y; % Number of finite elements n_gauss=n_gauss_side^2; % Number of Gauss points L_el_x=L_x/n_el_x; % Length of a finite element (x) L_el_y=L_y/n_el_y; % Length of a finite element (y) % ------------------------------------------------------------------------- % Pressure ---------------------------------------------------------------- n_np_x_p=n_el_x+1; % Number of nodal points (x) n_np_y_p=n_el_y+1; % Number of nodal points (y) n_np_p=n_np_x_p*n_np_y_p; % Number of nodal points dof_el_p=4; % Number of DOFs per element dof_p=n_np_p; % Total number of DOFs h_x_p=L_el_x; % Spatial step (x) h_y_p=L_el_y; % Spatial step (y) x_p=x_i:h_x_p:x_f; % Space vector (x) y_p=y_i:h_y_p:y_f; % Space vector (y) % ------------------------------------------------------------------------- % Velocity ---------------------------------------------------------------- n_np_x_v=2*n_el_x+1; % Number of nodal points (x) n_np_y_v=2*n_el_y+1; % Number of nodal points (y) n_np_v=n_np_x_v*n_np_y_v; % Number of nodal points dof_el_v=9; % Number of DOFs per element dof_v=n_np_v; % Total number of DOFs h_x_v=L_el_x/2; % Spatial step (x) h_y_v=L_el_y/2; % Spatial step (y) x_v=x_i:h_x_v:x_f; % Space vector (x) y_v=y_i:h_y_v:y_f; % Space vector (y) % ------------------------------------------------------------------------- % Body force -------------------------------------------------------------- n_np_x_b=n_el_x; % Number of nodal points (x) n_np_y_b=n_el_y; % Number of nodal points (y) n_np_b=n_np_x_b*n_np_y_b; % Number of nodal points h_x_b=L_el_x; % Spatial step (x) h_y_b=L_el_y; % Spatial step (y) x_b=x_i+L_el_x/2:h_x_b:x_f-L_el_x/2; % Space vector (x) y_b=y_i+L_el_y/2:h_y_b:y_f-L_el_y/2; % Space vector (y) % ------------------------------------------------------------------------- % Time -------------------------------------------------------------------- T=t_i:dt:t_f; % Time vector % ------------------------------------------------------------------------- ```

Gauss parameters

```gauss=[]; [gauss]=Gauss_parameters_2D(n_gauss_side,gauss); % Trasformation of coordinated for the Gauss integration points for n=1:n_gauss x_gauss(n)=L_el_x/2*(1+gauss(n).csi); y_gauss(n)=L_el_y/2*(1+gauss(n).eta); end % Jacobian of the transformation J_mat=[L_el_x/2 0 0 L_el_y/2]; J=det(J_mat); % Computation of shape and test functions (and derivatives) at Gauss points % for pressure [gauss]=shape_functions_Gauss_points_2D_p(gauss); [gauss]=test_functions_Gauss_points_2D_p(gauss); % Computation of shape and test functions (and derivatives) at Gauss points % for velocity [gauss]=shape_functions_Gauss_points_2D_v(gauss); [gauss]=test_functions_Gauss_points_2D_v(gauss); ```

Plot of geometry

```if strcmp(plot_geometry,'yes')==1 if n_el<=100 plot_geometry_2D(x_p,y_p,x_v,y_v,L_x,L_y,n_gauss,... L_el_x,L_el_y,x_gauss,y_gauss); end end ```

Plot of shape and test functions

```% Normalized domain csi_plot=-1:0.1:1; eta_plot=-1:0.1:1; if strcmp(plot_shape_functions,'yes')==1 % Shape functions for pressure N_plot_p=f_N_plot_2D_p(csi_plot,eta_plot); plot_shape_functions_2D(csi_plot,eta_plot,N_plot_p,dof_el_p); % Shape functions for velocity N_plot_v=f_N_plot_2D_v(csi_plot,eta_plot); plot_shape_functions_2D(csi_plot,eta_plot,N_plot_v,dof_el_v); end if strcmp(plot_test_functions,'yes')==1 % Test functions for pressure W_plot_p=f_W_plot_2D_p(csi_plot,eta_plot); plot_test_functions_2D(csi_plot,eta_plot,W_plot_p,dof_el_p); % Test functions for velocity W_plot_v=f_W_plot_2D_v(csi_plot,eta_plot); plot_test_functions_2D(csi_plot,eta_plot,W_plot_v,dof_el_v); end ```

Evaluate matrices and vectors

```% Afference matrix for pressure [A_p]=afference_matrix_2D_p(n_np_x_p,n_np_y_p,dof_el_p); % Afference matrix for velocity [A_v]=afference_matrix_2D_v(n_np_x_v,n_np_y_v,dof_el_v); % Element mass matrix for n=1:n_el el(n).M=element_mass_matrix_2D(dof_el_v,gauss,J); end % Strain rate - velocity matrix [gauss]=strain_rate_velocity_matrix_2D(dof_el_v,gauss,L_el_x,L_el_y); % Element viscosity matrix for n=1:n_el el(n).K=element_viscosity_matrix_2D(v,dof_el_v,gauss,J); end % Element gradient operator matrix for n=1:n_el el(n).G=element_gradient_operator_matrix_2D(dof_el_v,dof_el_p,gauss,... J,L_el_x,L_el_y); end % Element load vector for velocity for n=1:n_el [r,c]=row_column(n,n_el_x); x_c=x_i+L_el_x/2+L_el_x*(c-1); y_c=y_i+L_el_y*(n_el_y-1/2)-L_el_y*(r-1); if strcmp(body_force_selection,'centre')==1 b_c=feval(b_fun,x_c,y_c); el(n).f=element_load_vector_centre_2D(b_c,dof_el_v,gauss,J); elseif strcmp(body_force_selection,'gauss')==1 for k=1:n_gauss x_g(k)=x_c-L_el_x/2+x_gauss(k); y_g(k)=y_c-L_el_y/2+y_gauss(k); b_g(k,:)=feval(b_fun,x_g(k),y_g(k)); end el(n).f=element_load_vector_gauss_2D(b_g,dof_el_v,gauss,J); elseif strcmp(body_force_selection,'nodes_v')==1 x_n(1)=x_c-L_el_x/2; y_n(1)=y_c-L_el_y/2; x_n(2)=x_c+L_el_x/2; y_n(2)=y_c-L_el_y/2; x_n(3)=x_c+L_el_x/2; y_n(3)=y_c+L_el_y/2; x_n(4)=x_c-L_el_x/2; y_n(4)=y_c+L_el_y/2; x_n(5)=x_c; y_n(5)=y_c-L_el_y/2; x_n(6)=x_c+L_el_x/2; y_n(6)=y_c; x_n(7)=x_c; y_n(7)=y_c+L_el_y/2; x_n(8)=x_c-L_el_x/2; y_n(8)=y_c; x_n(9)=x_c; y_n(9)=y_c; for k=1:dof_el_v b_n(k,:)=feval(b_fun,x_n(k),y_n(k)); end el(n).f=element_load_vector_nodes_v_2D(b_n,dof_el_v,gauss,J); end end % Element load vector for pressure for n=1:n_el el(n).h=zeros(dof_el_p,1); end ```

Assemblate matrices and vectors

```% Assemblage of mass matrix [M]=assemble_mass_matrix(el,dof_v,n_el,dof_el_v,A_v); % Assemblage of viscosity matrix [K]=assemble_viscosity_matrix(el,dof_v,n_el,dof_el_v,A_v); % Assemblage of gradient operator matrix [G]=assemble_gradient_operator_matrix(el,dof_v,dof_p,n_el,dof_el_v,... dof_el_p,A_v,A_p); % Divergence operator matrix [G_T]=G'; % Zero matrix ZERO=zeros(dof_p,dof_p); % Zero_1 matrix ZERO_1=zeros(dof_v*n_sd,dof_p); % Zero_2 matrix ZERO_2=zeros(dof_p,dof_v*n_sd); % Assemblage of load vector for velocity [f]=assemble_load_vector_v(el,dof_v,n_el,dof_el_v,A_v); % Assemblage of load vector for pressure [h]=assemble_load_vector_p(el,dof_p,n_el,dof_el_p,A_p); ```

Conversion of the matrices from full to sparse

```M=sparse(M); K=sparse(K); G=sparse(G); G_T=sparse(G_T); f=sparse(f); h=sparse(h); ZERO=sparse(ZERO); ZERO_1=sparse(ZERO_1); ZERO_2=sparse(ZERO_2); ```

Boundary conditions

```% Definition of the constrained DOFs for the pressure dof_constrained_p=eval(dof_constrained_string_p); dof_free_p=dof_p-length(dof_constrained_p); dof_constrained_p=sort(dof_constrained_p); n_dof_constrained_p=length(dof_constrained_p); % Definition of the constrained DOFs for the velocity dof_constrained_v=eval(dof_constrained_string_v); dof_free_v=dof_v-length(dof_constrained_v); dof_constrained_v=sort(dof_constrained_v); n_dof_constrained_v=length(dof_constrained_v); dof_constrained_v_X_Y=sort([dof_constrained_v.*2-1,dof_constrained_v.*2]); dof_free_v_X_Y=dof_v*n_sd-length(dof_constrained_v_X_Y); n_dof_constrained_v_X_Y=length(dof_constrained_v_X_Y); % Evaluation of the derivative of the boundary conditions for the pressure bound_cond_p_der=0; % Evaluation of the derivative of the boundary conditions for the velocity bound_cond_v_der=0; % Evaluation of boundary conditions over time for pressure for m=1:length(T) time(m).t=T(m); for n=1:n_dof_constrained_p constrain_p(n)=bound_cond_p; constrain_p_der(n)=bound_cond_p_der; end time(m).p_p(:,1)=constrain_p'; time(m).p_der_p(:,1)=constrain_p_der'; end % Evaluation of boundary conditions over time for velocity for m=1:length(T) time(m).t=T(m); for n=1:n_dof_constrained_v_X_Y constrain_v(n)=bound_cond_v; constrain_v_der(n)=bound_cond_v_der; end time(m).u_p(:,1)=constrain_v'; time(m).u_der_p(:,1)=constrain_v_der'; end % Mass matrix [M_ff,M_fp,M_pf,M_pp]=constrain_matrix(M,... dof_constrained_v_X_Y,dof_constrained_v_X_Y); % Viscosity matrix [K_ff,K_fp,K_pf,K_pp]=constrain_matrix(K,... dof_constrained_v_X_Y,dof_constrained_v_X_Y); % Gradient operator matrix [G_ff,G_fp,G_pf,G_pp]=constrain_matrix(G,... dof_constrained_v_X_Y,dof_constrained_p); % Divergence operator matrix [G_T_ff,G_T_fp,G_T_pf,G_T_pp]=constrain_matrix(G_T,... dof_constrained_p,dof_constrained_v_X_Y); % Zero matrix [ZERO_ff,ZERO_fp,ZERO_pf,ZERO_pp]=constrain_matrix(ZERO,... dof_constrained_p,dof_constrained_p); % Zero_1 matrix [ZERO_1_ff,ZERO_1_fp,ZERO_1_pf,ZERO_1_pp]=constrain_matrix(ZERO_1,... dof_constrained_v_X_Y,dof_constrained_p); % Zero_2 matrix [ZERO_2_ff,ZERO_2_fp,ZERO_2_pf,ZERO_2_pp]=constrain_matrix(ZERO_2,... dof_constrained_p,dof_constrained_v_X_Y); % Load vector for velocity [f_f,f_p]=constrain_vector(f,dof_constrained_v_X_Y); % Load vector for pressure [h_f,h_p]=constrain_vector(h,dof_constrained_p); ```

Initial conditions

```% Initial conditions for the velocity time(1).iter(1).u=zeros(dof_v*n_sd,1,1); for n=1:n_np_v [r,c]=row_column(n,n_np_x_v); xp=x_i+L_el_x/2*(c-1); yp=y_i+L_el_y*n_el_y-L_el_y/2*(r-1); u_0_xp_yp=feval(u_0_fun,xp,yp); time(1).iter(1).u(n_sd*n-1,1)=u_0_xp_yp(1); time(1).iter(1).u(n_sd*n,1)=u_0_xp_yp(2); end time(1).iter(1).u_f=constrain_vector(time(1).iter(1).u,... dof_constrained_v_X_Y); % Pressure from the static problem ---------------------------------------- K_ff_inv=inv(K_ff); time(1).iter(1).p_f=(G_T_ff*K_ff_inv*G_ff)\... (... G_T_ff*K_ff_inv*(f_f-K_fp*time(1).u_p-G_fp*time(1).p_p)... +G_T_fp*time(1).u_p-h_f... ); [time(1).iter(1).p]=data_all_dof(time(1).iter(1).p_f,time(1).p_p,... dof_constrained_p); % ------------------------------------------------------------------------- time(1).iter(1).z=[time(1).iter(1).u_f;time(1).iter(1).p_f]; % Conversion of data from vector to matrix for j=1:n_np_y_v [time(1).iter(1).u_x_matrix(:,j)]=time(1).iter(1).u(... ((n_np_y_v-j)*n_np_x_v+1)*2-1:2:... ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-1); [time(1).iter(1).u_y_matrix(:,j)]=time(1).iter(1).u(... ((n_np_y_v-j)*n_np_x_v+1)*2-0:2:... ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-0); end time(1).iter(1).u_matrix=sqrt(time(1).iter(1).u_x_matrix.^2+... time(1).iter(1).u_y_matrix.^2); for j=1:n_np_y_p [time(1).iter(1).p_matrix(:,j)]=time(1).iter(1).p(... (n_np_y_p-j)*n_np_x_p+1:... (n_np_y_p-j)*n_np_x_p+n_np_x_p); end ```

```time(1).iter(1).D=K; [time(1).iter(1).D_ff,time(1).iter(1).D_fp,... time(1).iter(1).D_pf,time(1).iter(1).D_pp]=... constrain_matrix(time(1).iter(1).D,... dof_constrained_v_X_Y,dof_constrained_v_X_Y); disp('TIME INTEGRATION') disp('*******************************************************************') time(1).n_iter=1; tic for m=1:length(T)-1 % Initialization of variables time(m+1).iter(1).u=time(m).iter(end).u; time(m+1).iter(1).p=time(m).iter(end).p; % Initialization of unknowns [time(m+1).iter(1).u_f,time(m+1).iter(1).u_p]=... constrain_vector(time(m+1).iter(1).u,dof_constrained_v_X_Y); [time(m+1).iter(1).p_f,time(m+1).iter(1).p_p]=... constrain_vector(time(m+1).iter(1).p,dof_constrained_p); time(m+1).iter(1).z=[time(m+1).iter(1).u_f;time(m+1).iter(1).p_f]; % Initialization of Viscosity+Convection matrix time(m+1).iter(1).D=time(m).iter(time(m).n_iter).D; [time(m+1).iter(1).D_ff,time(m+1).iter(1).D_fp,... time(m+1).iter(1).D_pf,time(m+1).iter(1).D_pp]=... constrain_matrix(time(m+1).iter(1).D,... dof_constrained_v_X_Y,dof_constrained_v_X_Y); % Initialization of errors (>tollerances) time(m+1).iter(1).norm_err_p=toll_p+1; time(m+1).iter(1).norm_err_v=toll_v+1; % Convergence plots figure('Color',[1 1 1]) axes('FontSize',14) subplot(1,2,1,'YScale','log') hold on title(['Time step ',num2str(m),'/',num2str(length(T)-1)],'FontSize',14) xlabel('Iteration','FontSize',14) ylabel('Norm((p_k-p_k_-_1)/p_k)','FontSize',14) grid on grid minor subplot(1,2,2,'YScale','log') hold on title(['Time step ',num2str(m),'/',num2str(length(T)-1)],'FontSize',14) xlabel('Iteration','FontSize',14) ylabel('Norm((u_k-u_k_-_1)/u_k)','FontSize',14) grid on grid minor pause(eps) fprintf('\nTIME STEP %d/%d\n',m,length(T)-1) disp('-------------------------------------------------------------------') clear time_analysis time_inversion time_solution_p time_solution_v k=1; while time(m+1).iter(end).norm_err_p>toll_p ||... time(m+1).iter(end).norm_err_v>toll_v % Previous evaluation p_old=time(m+1).iter(k).p; u_old=time(m+1).iter(k).u; % Element convection matrix for n=1:n_el u_el=zeros(dof_el_v*n_sd,1); for a=1:dof_el_v for i=1:n_sd r=n_sd*(a-1)+i; u_el(r,1)=time(m+1).iter(k).u(A_v(n,r),1); end end el(n).C=element_convection_matrix_2D(u_el,dof_el_v,gauss,... J,L_el_x,L_el_y); end % Assemblage of convection matrix [C]=assemble_convection_matrix(el,dof_v,n_el,dof_el_v,A_v); % Conversion of the convection matrix from full to sparse C=sparse(C); % Viscosity+Convection matrix time(m+1).iter(k).D=K+C; % Constrain viscosity+convection matrix [time(m+1).iter(k).D_ff,time(m+1).iter(k).D_fp,... time(m+1).iter(k).D_pf,time(m+1).iter(k).D_pp]=... constrain_matrix(time(m+1).iter(k).D,... dof_constrained_v_X_Y,dof_constrained_v_X_Y); % ----------------------------------------------------------------- % Definition of matrices H=[M_ff,ZERO_1_ff;ZERO_2_ff,ZERO_ff]; Q2=[time(m+1).iter(k).D_ff,G_ff;G_T_ff,ZERO_ff]; Q1=[time(m).iter(time(m).n_iter).D_ff,G_ff;G_T_ff,ZERO_ff]; S2=[f_f;h_f]-... [M_fp,ZERO_1_fp;ZERO_2_fp,ZERO_fp]*... [time(m+1).u_der_p;time(m+1).p_der_p]-... [time(m+1).iter(k).D_fp,G_fp;G_T_fp,ZERO_fp]*... [time(m+1).u_p;time(m+1).p_p]; S1=[f_f;h_f]-... [M_fp,ZERO_1_fp;ZERO_2_fp,ZERO_fp]*... [time(m).u_der_p;time(m).p_der_p]-... [time(m).iter(time(m).n_iter).D_fp,G_fp;G_T_fp,ZERO_fp]*... [time(m).u_p;time(m).p_p]; % Evaluation of the new unknowns time(m+1).iter(k+1).z=(H+dt* theta *Q2)\... (... (H-dt*(1-theta)*Q1)*time(m).iter(time(m).n_iter).z+... dt* theta *S2+... dt*(1-theta)*S1... ); time(m+1).iter(k+1).u_f=... time(m+1).iter(k+1).z(1:dof_free_v_X_Y,1); time(m+1).iter(k+1).p_f=... time(m+1).iter(k+1).z(dof_free_v_X_Y+1:end,1); % ----------------------------------------------------------------- % Data for all dof [time(m+1).iter(k+1).u]=... data_all_dof(time(m+1).iter(k+1).u_f,time(m+1).u_p,... dof_constrained_v_X_Y); [time(m+1).iter(k+1).p]=... data_all_dof(time(m+1).iter(k+1).p_f,time(m+1).p_p,... dof_constrained_p); % Current evaluation u_new=time(m+1).iter(k+1).u; p_new=time(m+1).iter(k+1).p; % Evaluation of errors time(m+1).iter(k+1).norm_err_p=norm((p_new-p_old)/p_new); time(m+1).iter(k+1).norm_err_v=norm((u_new-u_old)/u_new); % Display results fprintf('\nIteration = %d',k-1) fprintf('\tNorm_p = %.2e',time(m+1).iter(k+1).norm_err_p) fprintf('\tNorm_v = %.2e',time(m+1).iter(k+1).norm_err_v) fprintf('\tElapsed time %.2f sec',toc) fprintf('\tReynolds = %.1f',... max(time(m+1).iter(k+1).u)*mean(L_x,L_y)/v) % Convergence plot update subplot(1,2,1) semilogy(k-1,time(m+1).iter(k+1).norm_err_p,'k+','LineWidth',2) plot([0,k+1],[toll_p,toll_p],'r','LineWidth',2) xlim([0,k+1]) ylim([min(toll_p/10,min(time(m+1).iter(k+1).norm_err_p)),... max(10,max(time(m+1).iter(k+1).norm_err_p)*10)]) subplot(1,2,2) semilogy(k-1,time(m+1).iter(k+1).norm_err_v,'k+','LineWidth',2) plot([0,k+1],[toll_v,toll_v],'r','LineWidth',2) xlim([0,k+1]) ylim([min(toll_v/10,min(time(m+1).iter(k+1).norm_err_v)),... max(10,max(time(m+1).iter(k+1).norm_err_v)*10)]) pause(eps) k=k+1; end disp(' ') disp('-------------------------------------------------------------------') time(m+1).n_iter=k-1; end disp(' ') disp('*******************************************************************') disp(' ') % Conversion of data from vector to matrix for m=1:length(T) for j=1:n_np_y_v [time(m).u_x_matrix(:,j)]=time(m).iter(1).u(... ((n_np_y_v-j)*n_np_x_v+1)*2-1:2:... ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-1); [time(m).u_y_matrix(:,j)]=time(m).iter(1).u(... ((n_np_y_v-j)*n_np_x_v+1)*2-0:2:... ((n_np_y_v-j)*n_np_x_v+n_np_x_v)*2-0); end time(m).u_matrix=sqrt(time(m).u_x_matrix.^2+time(m).u_y_matrix.^2); for j=1:n_np_y_p [time(m).p_matrix(:,j)]=time(m).iter(1).p(... (n_np_y_p-j)*n_np_x_p+1:... (n_np_y_p-j)*n_np_x_p+n_np_x_p); end end ```
```TIME INTEGRATION ******************************************************************* TIME STEP 1/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 1.60e+00 Norm_v = 1.32e+01 Elapsed time 1.90 sec Reynolds = 0.0 Iteration = 1 Norm_p = 1.43e-06 Norm_v = 1.11e-06 Elapsed time 3.54 sec Reynolds = 0.0 Iteration = 2 Norm_p = 1.50e-13 Norm_v = 1.15e-12 Elapsed time 5.15 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 2/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 7.36e-01 Norm_v = 4.80e+00 Elapsed time 7.02 sec Reynolds = 0.0 Iteration = 1 Norm_p = 7.74e-08 Norm_v = 3.72e-07 Elapsed time 8.65 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 3/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 4.32e-01 Norm_v = 2.47e+00 Elapsed time 10.51 sec Reynolds = 0.0 Iteration = 1 Norm_p = 8.00e-08 Norm_v = 2.08e-07 Elapsed time 12.14 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 4/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 2.76e-01 Norm_v = 1.45e+00 Elapsed time 13.99 sec Reynolds = 0.0 Iteration = 1 Norm_p = 4.32e-08 Norm_v = 1.29e-07 Elapsed time 15.62 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 5/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 1.83e-01 Norm_v = 9.07e-01 Elapsed time 17.47 sec Reynolds = 0.0 Iteration = 1 Norm_p = 2.49e-08 Norm_v = 8.39e-08 Elapsed time 19.11 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 6/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 1.23e-01 Norm_v = 5.87e-01 Elapsed time 21.01 sec Reynolds = 0.0 Iteration = 1 Norm_p = 1.56e-08 Norm_v = 5.58e-08 Elapsed time 22.63 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 7/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 8.32e-02 Norm_v = 3.87e-01 Elapsed time 24.52 sec Reynolds = 0.0 Iteration = 1 Norm_p = 1.02e-08 Norm_v = 3.75e-08 Elapsed time 26.16 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 8/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 5.66e-02 Norm_v = 2.59e-01 Elapsed time 28.05 sec Reynolds = 0.0 Iteration = 1 Norm_p = 6.85e-09 Norm_v = 2.55e-08 Elapsed time 29.70 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 9/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 3.86e-02 Norm_v = 1.75e-01 Elapsed time 31.64 sec Reynolds = 0.0 Iteration = 1 Norm_p = 4.65e-09 Norm_v = 1.73e-08 Elapsed time 33.28 sec Reynolds = 0.0 ------------------------------------------------------------------- TIME STEP 10/10 ------------------------------------------------------------------- Iteration = 0 Norm_p = 2.64e-02 Norm_v = 1.19e-01 Elapsed time 35.68 sec Reynolds = 0.0 Iteration = 1 Norm_p = 3.17e-09 Norm_v = 1.18e-08 Elapsed time 37.80 sec Reynolds = 0.0 ------------------------------------------------------------------- ******************************************************************* ```

Plots

```% Grid matrices [X_v,Y_v]=meshgrid(x_v,y_v); [X_p,Y_p]=meshgrid(x_p,y_p); [X_b,Y_b]=meshgrid(x_b,y_b); % Evaluation of the body force field for n=1:n_np_b [r,c]=row_column(n,n_np_x_b); xp=x_i+L_el_x*(c-1)+L_el_x/2; yp=y_i+L_el_y*n_el_y-L_el_y*(r-1)-L_el_y/2; b(n,:)=feval(b_fun,xp,yp); end b_x=b(:,1); b_y=b(:,2); for j=1:n_np_y_b [b_x_matrix(:,j)]=b_x((n_np_y_b-j)*n_np_x_b+1:... (n_np_y_b-j)*n_np_x_b+n_np_x_b); [b_y_matrix(:,j)]=b_y((n_np_y_b-j)*n_np_x_b+1:... (n_np_y_b-j)*n_np_x_b+n_np_x_b); end b_matrix=sqrt(b_x_matrix.^2+b_y_matrix.^2); % Limit values u_0_x_min=min(min(time(1).iter(1).u_x_matrix)); u_0_x_max=max(max(time(1).iter(1).u_x_matrix)); u_0_y_min=min(min(time(1).iter(1).u_y_matrix)); u_0_y_max=max(max(time(1).iter(1).u_y_matrix)); p_0_min=min(min(time(1).iter(1).p_matrix)); p_0_max=max(max(time(1).iter(1).p_matrix)); u_f_x_min=min(min(time(length(T)).u_x_matrix)); u_f_x_max=max(max(time(length(T)).u_x_matrix)); u_f_y_min=min(min(time(length(T)).u_y_matrix)); u_f_y_max=max(max(time(length(T)).u_y_matrix)); p_f_min=min(min(time(length(T)).p_matrix)); p_f_max=max(max(time(length(T)).p_matrix)); % Body force vector field figure('Color',[1 1 1]) axes('FontSize',14) quiver(X_b',Y_b',b_x_matrix,b_y_matrix,... 'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5) title('Body force vector field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Body force field figure('Color',[1 1 1]) axes('FontSize',14) contourf(X_b',Y_b',b_matrix,'LineWidth',1) title('Body force field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Initial condition of velocity in x figure('Color',[1 1 1]) axes('FontSize',14) surf(X_v',Y_v',time(1).iter(1).u_x_matrix) title('Initial condition of u_x','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('u_0_x(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) if u_0_x_min~=u_0_x_max zlim([u_0_x_min,u_0_x_max]) end % Initial condition of velocity in y figure('Color',[1 1 1]) axes('FontSize',14) surf(X_v',Y_v',time(1).iter(1).u_y_matrix) title('Initial condition of u_y','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('u_0_y(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) if u_0_y_min~=u_0_y_max zlim([u_0_y_min,u_0_y_max]) end % Initial condition of velocity vector field figure('Color',[1 1 1]) axes('FontSize',14) quiver(X_v',Y_v',time(1).iter(1).u_x_matrix,time(1).iter(1).u_y_matrix,... 'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5) title('Initial condition of velocity vector field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Initial condition of velocity field figure('Color',[1 1 1]) axes('FontSize',14) contourf(X_v',Y_v',time(1).iter(1).u_matrix,'LineWidth',1) title('Initial condition of velocity field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Initial condition of pressure figure('Color',[1 1 1]) axes('FontSize',14) surf(X_p',Y_p',time(1).iter(1).p_matrix) title('Initial condition of p','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('p(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) if p_0_min~=p_0_max zlim([p_0_min,p_0_max]) end % Solution of velocity in x figure('Color',[1 1 1]) axes('FontSize',14) surf(X_v',Y_v',time(length(T)).u_x_matrix) title('Solution of u_x','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('u_x(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) zlim([u_f_x_min,u_f_x_max]) % Solution of velocity in y figure('Color',[1 1 1]) axes('FontSize',14) surf(X_v',Y_v',time(length(T)).u_y_matrix) title('Solution of u_y','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('u_y(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) zlim([u_f_y_min,u_f_y_max]) % Solution of velocity vector field figure('Color',[1 1 1]) axes('FontSize',14) quiver(X_v',Y_v',time(length(T)).u_x_matrix,time(length(T)).u_y_matrix,... 'LineWidth',1,'Color',[1 0 0],'AutoScaleFactor',1.5) title('Solution of velocity vector field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Solution of velocity field figure('Color',[1 1 1]) axes('FontSize',14) contourf(X_v',Y_v',time(length(T)).u_matrix,'LineWidth',1) title('Solution of velocity field','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Solution of pressure figure('Color',[1 1 1]) axes('FontSize',14) surf(X_p',Y_p',time(length(T)).p_matrix) title('Solution of p','FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('p(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) zlim([p_f_min,p_f_max]) ```
```Warning: Contour not rendered for constant ZData ```

Display in command window

```disp('MODEL PARAMETERS') disp('-------------------------------------------------------------------') fprintf('Length (x)\t\t\t=\t%.1f\n',L_x) fprintf('Length (y)\t\t\t=\t%.1f\n',L_y) fprintf('Body force\t\t\t=\t%s\n',char(b_fun)) fprintf('Kinematic viscosity\t\t=\t%.2f\n',v) fprintf('DOFs constrained (velocity)\t=\t%s\n',dof_constrained_string_v) fprintf('DOFs constrained (pressure)\t=\t%s\n',dof_constrained_string_p) fprintf('Boundary conditions (velocity)\t=\t%.1f\n',bound_cond_v) fprintf('Boundary conditions (pressure)\t=\t%.1f\n',bound_cond_p) fprintf('Initial condition\t\t=\t%s\n',char(u_0_fun)) fprintf('Initial time\t\t\t=\t%.2f\n',t_i) fprintf('Final time\t\t\t=\t%.2f\n',t_f) fprintf('Time step\t\t\t=\t%.2e\n',dt) fprintf('Theta\t\t\t\t=\t%.2f\n',theta) disp('-------------------------------------------------------------------') disp(' ') disp('FEM PARAMETERS') disp('-------------------------------------------------------------------') fprintf('Number of finite elements (x)\t=\t%d\n',n_el_x) fprintf('Number of finite elements (y)\t=\t%d\n',n_el_y) fprintf('Length of a finite element (x)\t=\t%.2f\n',L_el_x) fprintf('Length of a finite element (y)\t=\t%.2f\n',L_el_y) fprintf('Gauss integration points\t=\t%d\n',n_gauss) fprintf('Number of nodes (velocity)\t=\t%d\n',n_np_v) fprintf('Number of nodes (pressure)\t=\t%d\n',n_np_p) fprintf('DOFs per element (velocity)\t=\t%d\n',dof_el_v) fprintf('DOFs per element (pressure)\t=\t%d\n',dof_el_p) fprintf('Total number of DOF (velocity)\t=\t%d\n',dof_v) fprintf('Total number of DOF (pressure)\t=\t%d\n',dof_p) disp('-------------------------------------------------------------------') ```
```MODEL PARAMETERS ------------------------------------------------------------------- Length (x) = 1.0 Length (y) = 1.0 Body force = @(x,y)[(12-24*y)*x^4+(-24+48*y)*x^3+(-48*y+72*y^2-48*y^3+12)*x^2+(-2+24*y-72*y^2+48*y^3)*x+(1-4*y+12*y^2-8*y^3),(8-48*y+48*y^2)*x^3+(-12+72*y-72*y^2)*x^2+(4-24*y+48*y^2-48*y^3+24*y^4)*x+(-12*y^2+24*y^3-12*y^4)] Kinematic viscosity = 1.00 DOFs constrained (velocity) = [2:n_np_x_v-1,n_np_v-n_np_x_v+2:n_np_v-1,1:n_np_x_v:n_np_v-n_np_x_v+1,n_np_x_v:n_np_x_v:n_np_v] DOFs constrained (pressure) = [1] Boundary conditions (velocity) = 0.0 Boundary conditions (pressure) = 0.0 Initial condition = @(x,y)[0,0] Initial time = 0.00 Final time = 0.10 Time step = 1.00e-02 Theta = 1.00 ------------------------------------------------------------------- FEM PARAMETERS ------------------------------------------------------------------- Number of finite elements (x) = 10 Number of finite elements (y) = 10 Length of a finite element (x) = 0.10 Length of a finite element (y) = 0.10 Gauss integration points = 9 Number of nodes (velocity) = 441 Number of nodes (pressure) = 121 DOFs per element (velocity) = 9 DOFs per element (pressure) = 4 Total number of DOF (velocity) = 441 Total number of DOF (pressure) = 121 ------------------------------------------------------------------- ```

Animation

```% Screen dimensions scrsz=get(0,'ScreenSize'); bar=64; % Simulation parameters fact_ampl_sim=1; fact_vel_sim=(t_f-t_i)/animation_time; interval=0.0001; i=1; j=1; tic figure('Color',[1 1 1],'Position',[0 0 scrsz(3) (scrsz(4)-bar)]) while j<=length(T) % Plot quiver(X_v',Y_v',time(j).u_x_matrix,time(j).u_y_matrix,... 'LineWidth',1,'Color',[0 0 0],'AutoScaleFactor',1.5) title(['Solution - t = ',num2str(round(T(j)*100)/100),' sec'],... 'FontSize',14) xlabel('x','FontSize',14) ylabel('y','FontSize',14) zlabel('u(x,y)','FontSize',14) grid on grid minor xlim([x_i,x_f]) ylim([y_i,y_f]) % Time calibration i=i+1; time_real=toc; time_simulation=T(j)/fact_vel_sim; j=ceil(time_real/dt*fact_vel_sim+eps); pause(interval); end relative_error=abs(time_simulation-time_real)/abs(time_real); ```