function problem3_bvp(k_input) %solves the problem 3 of hw 8 (6.C.5 in beers book) % clear all; % clc; nz = 20; %number of grid points in the z direction ny = 25; %number of grid points in the y direction k = k_input; %rate constant in 1/M/s K = 1e3; %equilibrium constant Ha = 1e3; %henrys constant for A mol/m^3/atm pa = 1e-4; %partial pressure of A mu = 1e-3; %viscosity of water rho = 1e3; %density of water Da = 1e-5; %diffusivity of A Db = 1e-5; %diffusivity of B Dc = 1e-5; %diffusivity of C b = 0.1; %width of channel L = 50; %length of channel theta = 80; %incline of the channel cb_init = 1;%concentration of b in the reservoir %initialize the grid points on y direction for i=1:20 y_i(i) = (i-1)*0.01/19; end for i=21:ny y_i(i) = 0.01 + (b-0.01)/(ny-20)*(i-20); end %intitialize the grid points on z direction for i=1:10 z_i(i) = (i-1)*5/(10-1); end for i=11:20 z_i(i) = 5+ (i-10)*(L-5)/(nz-10); end %populate the params params.nz = nz;params.k = k; params.K = K; params.Ha = Ha; params.pa = pa; params.mu = mu; params.rho = rho; params.Da = Da; params.Db = Db; params.Dc = Dc; params.b = b; params.L = L; params.theta = theta; params.cb_init = cb_init; params.ny = ny; params.z_i = z_i; params.y_i = y_i; %give the initial guesses of the profiles of each species ca0 = pa*Ha*ones(ny*(nz),1); cb0 = cb_init*ones(ny*(nz),1); cc0 = zeros(ny*(nz),1); %make a initial guess vector with all concentrations one below the other c_guess = [ca0;cb0;cc0]; %solve the system of equations options = optimset('Display','iter'); c = fsolve(@solve_c_nonuniformgrid,c_guess,options,params); %get concentrations of individual species from the ouput concentration %vector ca = c(1:ny*(nz)); cb = c(ny*(nz)+1:2*ny*(nz)); cc = c(2*ny*(nz)+1:3*ny*(nz)); %make concentration matrices for each species so that they can be plotted ca_matrix = array_to_matrix(ca,ny,nz); cb_matrix = array_to_matrix(cb,ny,nz); cc_matrix = array_to_matrix(cc,ny,nz); %plot the concentration of each species pcolor(y_i',z_i',ca_matrix);xlabel('y-axis (cm)'); ylabel('z-axis (cm)');title('Concentration of A');shading interp;colorbar;%zlabel('Concentration of A (mol/cm3)'); figure; pcolor(y_i',z_i',cb_matrix);xlabel('y-axis (cm)'); ylabel('z-axis (cm)');title('Concentration of B');shading interp;colorbar;%zlabel('Concentration of B (mol/cm3)'); figure; pcolor(y_i',z_i',cc_matrix);xlabel('y-axis (cm)'); ylabel('z-axis (cm)');title('Concentration of AB');shading interp;colorbar;%zlabel('Concentration of AB (mol/cm3)'); %calculate the flux of A flux = 0; for i=2:nz fi = (-3*ca_matrix(i,1)+4*ca_matrix(i,2)-ca_matrix(i,3))/(y_i(2)-y_i(1))/2; fi_1 = (-3*ca_matrix(i-1,1)+4*ca_matrix(i-1,2)-ca_matrix(i-1,3))/(y_i(2)-y_i(1))/2; flux = flux + -Da/1000*(fi+fi_1)/2*(z_i(i)-z_i(i-1))/L; end output = strcat('The total flux of A into the channel is:\t', num2str(flux),'\t mol/cm^2/s\n'); fprintf(output); return;