Drug delivery in circulatory system

生物醫學研究中心的一組科學家被要求研究藥物在循環系統中的轉運。在初始階段，他們設計了一個簡單的場景來校準他們的初始設計。他們雇用你來協助他們的研究通過對一個二維板進行數值研究。你知道，要研究這種輸運，你必須解下列對流-擴散方程

其中為藥物的濃度，為速度，為藥物在血液中的擴散系數常數。血流是單向的(在上式中只存在分量)，其速度由下式模擬

結束時,你的工作,你需要提供一個壓縮文件夾包含所有您已經實現的功能除了主要腳本運行整個仿真的接口(你必須加載變量直接在主腳本,不允許用戶動態修復它們)。所有的編碼都需要用適當的注釋來解釋。此外，您必須提供一個PDF文件，其中包含以下結果和討論。

利用MATLAB的contourf命令繪制以下時間點的藥物沿平板分布:t=0,0.5,1,1.5,2,2.5,3s。討論結果和主要特征的觀察物理行為的濃度。

MATLAB求解代碼

clc close all clear %% initalisation lx=4*pi; y_max=pi; y_m=0.5*y_max; D=0.1; dt=0.1; nx=51; ny=21; tmax=3; %% intitalisation of x,y,t dx=lx/(nx-1); x=0:dx:lx; dy=y_max/(ny-1); y=0:dy:y_max; nt=tmax/dt+1; t=0:dt:tmax;%% Processing Velocity u=zeros(nx,ny,nt); for i=1:nxfor j=1:nyfor k=1:ntu(i,j,k)=3*(1-((x(i)-2*pi)^2/(2*pi^2))*cos(pi*t(k))*sin(x(i)))*1-((y(j)-y_m)^2/y_m^2);endend end %% concentration at inital condition c=zeros(nx,ny,nt); for i=1:4nxfor j=1:nyc(i,j,1)=1-(y(j)-y_m)^2/(y_m^2);end end %% Values for AP, AN, AS, AW, AE from discretisationas=-D*(dt)/(dy^2); ap= 1+(2*D*dt/dy^2)+(2*D*dt/dx^2); an=-D*(dt)/dy^2;b=reshape(c(:,2:ny-1,1),[nx*(ny-2),1]); f(:,1)=b; %stores the value of concentration A = zeros(nx*(ny-2),nx*(ny-2)); %creates a specific size of matrix Afor k=2:nt for i=1:nxfor j=2:ny-1ae(i,j,k)=u(i,j,k)*dt/(2*dx)-D*dt/(dx^2); %depends on the U velocity vectoraw(i,j,k)=-u(i,j,k)*dt/(2*dx)-D*dt/(dx^2); %''endend%%Types%T1 centrefor i=2:(nx-1)for j=3:ny-2pointer=(j-2)*nx+i; A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer+nx)=an;A(pointer,pointer-nx)=as; endend%T2 leftfor j=3:ny-2i=1;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer+nx-2)=aw(i,j,k);A(pointer,pointer+nx)=an;A(pointer,pointer-nx)=as; end%T3 rightfor j=3:ny-2i=nx;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer-nx+2)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer+nx)=an;A(pointer,pointer-nx)=as; end%T4 bottomfor i=2:nx-1j=2;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer+nx)=an;end%T5 topfor i=2:nx-1j=ny-1;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer-nx)=as; end%T6 bottom lefti=1;j=2;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer+nx-2)=aw(i,j,k);A(pointer,pointer+nx)=an;%T7 bottom righti=nx;j=2;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer-nx+2)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer+nx)=an;%T8 top lefti=1;j=ny-1;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer+1)=ae(i,j,k);A(pointer,pointer+nx-2)=aw(i,j,k);A(pointer,pointer-nx)=as;%T9 top righti=nx;j=ny-1;pointer=(j-2)*nx+i;A(pointer,pointer)=ap;A(pointer,pointer-nx+2)=ae(i,j,k);A(pointer,pointer-1)=aw(i,j,k);A(pointer,pointer-nx)=as;%solvingZ=A\b; f(:,k)=Z; b=Z; end%% Resultsc(:,2:ny-1,:)=reshape(f,[nx,ny-2,nt]);%reshapes the concentration values ffor h=1:5:31 figurecontourf(c(:,:,h)')xlabel('x (m)')ylabel('y (m)')time=(h-1)*dt;title([ 'Concentration over time' , num2str(time),'s']) end

總結

以上是生活随笔為你收集整理的在循环系统中给药的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。