sandbox/easystab/stab2014/free_surface_gravity_total_pressure

    Free-surface-gravity: displaying total pressure

    In this code, we take the code that display the animation of the free surface waves free_surface_gravity_particles.m and we show total pressure instead of only the pressure disturbance (it means that we add hydrostatic pressure to the pressure field) to see what changes.

    In a fluid, hydrostatic pressure is determined by the formula below: \displaystyle p=\rho.g.y where:

    • p is the hydrostatic pressure (Pa)

    • \rho is the fluid density (kg/m^{3})

    • g is the gravational acceleration (m/s^{2})

    • y is the height of the point (m)

    So we will add the hydrostatic pressure p to our pressure field to obtain the total pressure in our case.

    Dependency:

    • chebdif.m for the Chebychev differentiation matrices

    clear all; clf; n=100; % number of gridpoints alpha=1; % wavenumber in x L=2; % Fluid height in y rho=1; % fluid density mu=0.0001; % fuid viscosity g=1; % gravity % differentiation matrices scale=-2/L; [y,DM] = chebdif(n,2); D=DM(:,:,1)scale;
    DD=DM(:,:,2)    scale^2;
    y=(y-1)/scale; I=eye(n); Z=zeros(n,n); % renaming the matrices dy=D; dyy=DD; dx=ialphaI; dxx=-alpha^2I; Delta=dxx+dyy; % System matrices A=[muDelta, Z, -dx, Z(:,1); … Z, muDelta, -dy, Z(:,1); … dx, dy, Z, Z(:,1); … Z(1,:),I(n,:),Z(1,:),0]; E=blkdiag(rhoI,rhoI,Z,1); % boundary conditions loc=[1,n,n+1,2*n]; C=[I(1,:),Z(1,:),Z(1,:),0; … Z(1,:),I(1,:),Z(1,:),0; … Z(1,:),Z(1,:),-I(n,:),rho*g; … dy(n,:),dx(n,:),Z(1,:),0]; E(loc,:)=0;
    A(loc,:)=C; % compute eigenmodes [U,S]=eig(A,E); s=diag(S); [t,o]=sort(-real(s)); s=s(o); U=U(:,o); rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % particle animation of the eigenmodes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % parameters nper=1; % number of periods of oscillation nt=30; % number of time steps per period nx=40; % number of points in x modesel=1; % which more do animate % select the eigenmode u=1:n; v=u+n; p=v+n; eta=3    n+1; q=U(:,modesel); lambda=s(modesel); % The time and x extent tvec=linspace(0,nper2pi/abs(lambda),npernt); dt=tvec(2)-tvec(1); Lx=2pi/alpha; x=linspace(0,Lx,nx); % scale mode amplitude q=0.05q/abs(q(eta)); % initialize tracer particles [px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); py=py.(1+2real(exp(lambdatvec(1))q(eta)exp(ialphapx))/L); px=px(:);py=py(:);
    figure(1) % time loop for ind=1:npernt % expand mode to physical space qq=2real(exp(lambdatvec(ind))qexp(ialphax)); % plot pressure surf(x,y,qq(p,:)-10+kron(ones(1,nx),rhogy),‘facealpha’,0.3); view(2); shading interp; hold on % plot free surface plot(x,L+qq(eta,:),‘k-’,x,0x+L,‘k–’); hold on %%%% plot the particles plot(mod(px,Lx),py,‘k.’); pu=interp1(y,q(u),py); pv=interp1(y,q(v),py); % For particles above L, use Taylor expansion for velocity pu(py>L)=q(u(n))+q(eta)dy(n,:)q(u); pv(py>L)=q(v(n))+q(eta)dy(n,:)q(v); % expand to physical space puu=2real(exp(lambdatvec(ind))pu.exp(ialphapx)); pvv=2real(exp(lambdatvec(ind))pv.exp(ialphapx)); % advect particles px=px+puudt; py=py+pvvdt; xlabel(‘x’);
    ylabel(‘y’); title(‘Free surface gravity with total pressure - Mode 1’); axis equal; axis([0,Lx,0,1.3*L]); grid off hold off drawnow end set(gcf,‘paperpositionmode’,‘auto’); print(‘-dpng’,‘-r100’,‘free_surface_gravity_total_pressure.png’); figure(2) % initialize tracer particles [px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); py=py.(1+2real(exp(lambdatvec(1))q(eta)exp(ialphapx))/L); px=px(:);py=py(:); for ind=1:npernt % expand mode to physical space qq=2real(exp(lambdatvec(ind))qexp(ialphax)); % plot pressure surf(x,y,qq(p,:)-10,‘facealpha’,0.3); view(2); shading interp; hold on % plot free surface plot(x,L+qq(eta,:),‘k-’,x,0x+L,‘k–’); hold on %%%% plot the particles plot(mod(px,Lx),py,‘k.’); pu=interp1(y,q(u),py); pv=interp1(y,q(v),py); % For particles above L, use Taylor expansion for velocity pu(py>L)=q(u(n))+q(eta)dy(n,:)q(u); pv(py>L)=q(v(n))+q(eta)dy(n,:)q(v); % expand to physical space puu=2real(exp(lambdatvec(ind))pu.exp(ialphapx)); pvv=2real(exp(lambdatvec(ind))pv.exp(ialphapx)); % advect particles px=px+puudt; py=py+pvv*dt; xlabel(‘x’);
    ylabel(‘y’); title(‘Free surface gravity with pressure disturbance - Mode 1’); axis equal; axis([0,Lx,0,1.3*L]); grid off hold off drawnow
    end set(gcf,‘paperpositionmode’,‘auto’); print(‘-dpng’,‘-r100’,‘free_surface_gravity_pressure_disturbance.png’);

    Figures

    Pressure disturbance

    Pressure disturbance

    Total pressure

    Total pressure

    So we see that because the value of hydrostatic pressure is much higher than the pressure disturbance so the total pressure is almost the same as the hydrostatic pressure, it only depends on the height of the points. %}