Need some guidance to solve this problem for velocity and pressure variables

Hello all,

In my problem, I need to solve the following equations to calculate the variables velocity (from Eq .(1 and 2)) and pressure (from Eq .(3)).

∇.G=S (1)

∇×G-G×∇(ln(F1+F2*G))=0 (2)

∇^2 (P^2)=-(F1+F2*G)S - G∙∇(F1+F2*G) (3)

where, G, G = mass velocity vector and its magnitude, P = Pressure and S = source term.

Here is my code below and queries to follow:

1) While defining continuity Eq. (1), I'm facing an error,

% ERROR: Array indices must be positive integers or logical values.

2) Without Eq. (1), the velocity and pressure values obtained are not good as expected.

I'm not sure about my approach as how to find the appropriate velocity and pressure.

Someone kindly examine the code and judge my approach and assist me further.

Thank you

clear; clc; close all;

L=0.25; W =5;
Nx=10; Ny=200;

x = linspace(0,L,Nx);
y = linspace(0,W,Ny);
[X, Y] = meshgrid(x,y);
dx = x/Nx; dy = y/Ny;

rhoo = 0.8;
por = 0.3; mu =0.8E-5; dia = 0.003;
P = 101325.*ones(Ny,Nx); 
V = zeros(Ny,Nx);
Gvec = zeros(Ny,Nx);

TimeStep =1;
Vx = gradient(X, TimeStep);
Vy = gradient(Y, TimeStep);
Vmag = sqrt(Vx.^2 + Vy.^2);

% Boundary conditions
V(:,1) = 0;
V(Ny,:) = 0;
V(1,:) = V(Ny-1,:);
V(:,Nx) = V(:,Nx-1);
V(Ny,1) = 0;
V(1,1)  = 0;

S = [0   0   0.1458   0.2058   0.2161   0.2171    0   0   0   0];

f1 = (300*((1-por).^2)*P*mu)/(rhoo.*(dia.^2)*(por.^3));
f2 = (3.5.*(1-por).*P)/(rhoo.*dia.*(por.^3));
F1 = f1.*ones(Ny,Nx);
F2 = f2.*ones(Ny,Nx);
Gvec =rhoo.*V;
Gmag =rhoo.*Vmag;
%divergence(Gvec) = 0        % Eq. (1)
% ERROR: Array indices must be positive integers or logical values.

d1 = log (F1 + (F2.*Gmag));
d2 = gradient(d1);

A = curl(Gvec,X);
B = curl(d2,Gvec);
V = A - B              % Eq. (2)

V = 200×10

0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139 0.0139

d3 = F1 + (F2.*Gmag);
d4 = gradient(d3);
P  = gradient(gradient((P.^2))) + divergence(d3,Gvec) + d3.*S   % Eq. (3)

P = 200×10

1.0e+08 * 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000 0 0 1.0178 1.4367 1.5086 1.5155 0 0 0.0000 0.0000

[Torsten]

Your code above just assumes values for G and P - you must solve for them as they should satisfy your differential equations.

I think it's best if you use a commercial solver, e.g. ANSYS or COMSOL for the equations of fluid dynamics.

OPENFOAM also has prepared code for fluid dynamics problems.

Although I only recognize the continuity equation - pressure and momentum equations look quite strange.

[Kumaresh Kumaresh]

May I know why you are guiding to use black box tools like ANSYS or COSMOL ?

Is that impossible to solve by MATLAB ? Kindly elaborate my mistakes. I would like to learn and solve.

https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.690200519

Hope this link helps to understand about the pressure momentum equation.

Thank you for your response

[Torsten]

May I know why you are guiding to use black box tools like ANSYS or COSMOL ?

The laminar Navier Stokes equations need a special treatment because of the velocity-pressure coupling. They are hard to solve for a beginner. Or how much time do you plan to invest in reading before you start programming ?

Is that impossible to solve by MATLAB ?

If it can be done in OPENFOAM, ANSYS and COMSOL, it can also be done in MATLAB. But you will have to program it from scratch if you don't find code somewhere else. So again the question: how much time do you plan to invest ? I would estimate half a year if you have to do all by yourself.

[Kumaresh Kumaresh]

Before working on this problem in MATLAB, I tried implementing it in OpenFOAM but failed to achieve the expected outcome. As you mentioned, the above equations are coupled sets of equations, which need special treatment. I don't find code anywhere or any default solvers related to it. This is a work by Ergun known by his name as the Ergun equation to examine the flow distribution in a porous medium.

Yes, I'm a beginner in MATLAB and this is part of my thesis work. I don't have too much time to invest in it. However, I will do my best in all possible ways to make it possible.

Kindly guide me.

Thank you for your valuable comments.