Finite Differences using fsolve for non linear equations

8 次查看（过去 30 天）
Ellson Chow
Ellson Chow 2020-8-21
评论： Abdallah Daher 2022-3-15
I need help writing a code that solves the boundary layer problem for fluid mechanics. I am using fsolve to solve non-linear equations (continuity and momentum equations) for the boundary layer, but I am stuck getting no solutions for my code. Any help is greatly appreciated! Here I have attempted to convert my function F (equations) and its variables, u and v, into a single vector (f and x respectively) in an attempt to solve it as I assume that fsolve takes in vector inputs. Total 2*ny*nx equations and hence unknowns to be solved, 1 for each cell.
My code:
clear all, close all, clc
nx = 10;
ny = 10;
x0 = ones(2*ny*nx,1);
k = fsolve(@blayer, x0);
function f = blayer(x)
nx = 10;
ny = 10;
L = 1;
H = 1;
dx = L/(nx-1);
dy = H/(ny-1);
P = 1;
pgrad = P/L;
v = zeros(ny,nx);
u = zeros(ny,nx);
u(:,1) = 1; % Incoming
v(:,1) = 0; % Incoming
u(1,:) = 0; % Bottom
v(1,:) = 0; % Bottom
u(ny,:) = 0; % Top
v(ny,:) = 0; % Top
F = zeros(2*ny, nx);
X = [u;v];
Xt = X';
x = Xt(:);
j = 0;
while j <= ny
j = j+1;
if j == 1
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
end
end
end
if j > 1 && j < ny
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
end
end
end
if j == ny
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
end
end
end
end
Ft = F';
f = Ft(:);
end
  1 个评论
Ellson Chow
Ellson Chow 2020-8-21
I have changed my initial guess to u = 1 and v = 0;
x0 = zeros(2*ny*nx,1);
x0(1:ny*nx) = 1;
but still does not work.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2020-8-21
Do you need fsolve at all? Doesn't your blayer function do the solving itself? Couldn't you simply have
nx = 10;
ny = 10;
L = 1;
H = 1;
dx = L/(nx-1);
dy = H/(ny-1);
P = 1;
pgrad = P/L;
v = zeros(ny,nx);
u = zeros(ny,nx);
u(:,1) = 1; % Incoming
v(:,1) = 0; % Incoming
u(1,:) = 0; % Bottom
v(1,:) = 0; % Bottom
u(ny,:) = 0; % Top
v(ny,:) = 0; % Top
F = zeros(2*ny, nx);
% X = [u;v];
% Xt = X';
% x = Xt(:);
j = 0;
while j <= ny
j = j+1;
if j == 1
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
end
end
end
if j > 1 && j < ny
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
end
end
end
if j == ny
for i = 1:nx
if i == 1
F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i > 1 && i < nx
F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
elseif i == nx
F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j,i)-v(j-1,i))/dy;
F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
end
end
end
end
surf(F)
or am I missing something?
  3 个评论
显示 1更早的评论
Alan Stevens
Alan Stevens 2020-8-21
So you want to find the values of u and v that make F zero everywhere? In that case you need to pass u and v to blayer, not x (the values of which are fixed). i.e. you want fsolve to find values of u and v that set all the Fs to zero.
Abdallah Daher
Abdallah Daher 2022-3-15
Were you able it out? I am currently trying to write an explicit finite difference scheme for a flat plate boundary layer

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Mechanics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by