Plotting CFD data contour in matlab

65 views (last 30 days)
Hi,
I have imported CFD data into matlab and I can see x,y,z and U column in my data and I seperate them in a column matrix. I want to plot the contour of this data. I am cluless here. I have tried many approaches. To my understanding about plotting a contour:
step1: plot x and y according to grid size(n) required and draw a meshgrid. [X,Y]=meshgrid(x,y). this will generate X and Y of n by n.
Step 2: constuct "V" of n by n matrix of velocity by rearranging the data.
Step 3: Plot contour using pcolor(x,y,V) or contour(x,y,V)
However I am stuck on step2 how to convert V into n by n order. I used reshape() function which is giving error (Index in position 1 is invalid. Array indices must be positive integers or logical values.
U_reshaped= reshape(U(zeros(200^2,1), X(:), Y(:)), [200 200]); % I want 200 by 200 however size of my U is only 33000 by 1. I have seen this line used by one code which is working fine in his code in this format. but when I am trying on my data it is giving me this error.
Is there any easy way to visualise data in Matlab or any help with this command. Thanks
  4 Comments
Muhammad Atif
Muhammad Atif on 18 Aug 2021
Hi Wan and KSSV,
Thanks for your reply.
My model is in 3D but I only extracted data in a 2D plane. so one axis is constant (it should be z but somehow ansys shuffeled y and z but i know the plane was in x-y )
I have attached my data.
Nikhil Shinde
Nikhil Shinde on 19 Aug 2022
Edited: Nikhil Shinde on 19 Aug 2022
Please take a look at Delaunay and Trisurf functions in matlab. Extract your vertices data in a matrix. Following is the code that I used in my posprocessing assignment, You can tailor it according to your needs
vel_files = dir('*.csv'); %Only velocity data
velnum = length(vel_files);
vel_data = cell(1, velnum);
for k = 1:velnum-1
filename = vel_files(k).name;
vel_data{k} = readmatrix(filename);
end
Tvert = readtable('vertices.csv'); %vertices data
x = Tvert.Var1;
y = Tvert.Var2;
vertices= cat(2,x ,y);
DT = delaunay(vertices);% Delaunay Triangulation
VELALL=cell2mat(vel_data); trisurf(DT,vertices(:,1),vertices(:,2),VELALL(:,141),'EdgeColor','none'); axis tight; daspect([1 1 1]); shading interp; colormap jet; view(2);title('Ux t=100s')

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 18 Aug 2021
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714947/snapshot_vel-0001.txt', 'VariableNamingRule','preserve')
T1 = 28947×7 table
nodenumber x-coordinate y-coordinate z-coordinate x-velocity y-velocity z-velocity __________ ____________ ____________ ____________ __________ __________ __________ 1 -50 -3.35 1.4744e-13 26.822 0 0 2 -49.575 -3.35 -4.3055e-09 26.822 0 0 3 -49.017 -3.35 -6.2855e-10 26.822 0 0 4 -48.298 -3.35 1.1733e-09 26.822 0 0 5 -47.532 -3.35 -3.9007e-10 26.822 0 0 6 -47.347 -3.35 2.0703e-09 26.822 0 0 7 -47.068 -3.35 8.5003e-09 26.822 0 0 8 -46.248 -3.35 -8.7366e-09 26.822 0 0 9 -45.919 -3.35 3.4381e-10 26.822 0 0 10 -45.776 -3.35 -1.3321e-08 26.822 0 0 11 -44.764 -3.35 2.1362e-08 26.822 0 0 12 -44.228 -3.35 8.2932e-09 26.822 0 0 13 -43.675 -3.35 5.3661e-09 26.822 0 0 14 -42.967 -3.35 -7.7576e-09 26.822 0 0 15 -42.153 -3.35 1.3898e-08 26.822 0 0 16 -41.239 -3.35 -6.0162e-09 26.822 0 0
[Uy,ix1] = unique(T1.('y-coordinate'))
Uy = 26085×1
-3.3500 -3.3500 -3.3500 -3.3500 -3.3500 -3.3500 -3.3500 -3.3500 -3.3500 -3.3500
ix1 = 26085×1
49 54 238 50 178 182 52 90 48 250
Unfortunately, the data do not appear to be gridded (here, regular intervals in ‘ix1’), so it will be necessary to interpolate. There are two likely options, griddata and scatteredInterpolant. This approach uses griddata:
N = 50; % Controls Grid Resolution
xv = linspace(min(T1.('x-coordinate')), max((T1.('x-coordinate'))), N);
yv = linspace(min(T1.('y-coordinate')), max((T1.('y-coordinate'))), N);
[Xc,Yc] = ndgrid(xv, yv);
Xv = griddata(T1.('x-coordinate'), T1.('y-coordinate'), T1.('x-velocity'), Xc, Yc);
figure
contourf(Xc, Yc,Xv)
figure
surfc(Xc, Yc, Xv)
This illustrates the general approach.
If you want to include a third dimension for the contour, it will be necessary to use the isosurface function. I leave those details to you.
.
  2 Comments
Muhammad Atif
Muhammad Atif on 18 Aug 2021
Hi Start,
I am truly grateful for your comment. It was really helpful. I have used griddata() command and it is making some sense of the results that I am getting but it is not capturing the geometry boundary properly. I have tried different plotting fuctions such as contour, contourf, surfc and pcolor. Last two are showing better boundary display but the boundary is not as smooth as it was in ansys and also it looks like draged in y axis and squeshed in x axis.
I have attached CFD and Matlab pictures. Any help improving the boundary capture and scaling the body properly would be highly appritiating.
Below is the code I used to generate meshgrid and griddata and contours.
X=raw_data.data(:,2);
Y=raw_data.data(:,3);
U=raw_data.data(:,5);
V=raw_data.data(:,6);
W=raw_data.data(:,7);
X1=min(X):0.1:max(X); Y1=min(Y):0.1:max(Y);
[xmesh, ymesh] = meshgrid(X1,Y1);
[x1 y1 u1]=griddata(X,Y,U,xmesh,ymesh);
[x1 y1 v1]=griddata(X,Y,V,xmesh,ymesh);
[x1 y1 w1]=griddata(X,Y,W,xmesh,ymesh);
res_vel=sqrt(u1.^2+v1.^2+w1.^2);
contourf(x1,y1,res_vel)
shading interp;
Thanks once again.
Star Strider
Star Strider on 19 Aug 2021
I adapted your code to mine:
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714947/snapshot_vel-0001.txt', 'VariableNamingRule','preserve')
T1 = 28947×7 table
nodenumber x-coordinate y-coordinate z-coordinate x-velocity y-velocity z-velocity __________ ____________ ____________ ____________ __________ __________ __________ 1 -50 -3.35 1.4744e-13 26.822 0 0 2 -49.575 -3.35 -4.3055e-09 26.822 0 0 3 -49.017 -3.35 -6.2855e-10 26.822 0 0 4 -48.298 -3.35 1.1733e-09 26.822 0 0 5 -47.532 -3.35 -3.9007e-10 26.822 0 0 6 -47.347 -3.35 2.0703e-09 26.822 0 0 7 -47.068 -3.35 8.5003e-09 26.822 0 0 8 -46.248 -3.35 -8.7366e-09 26.822 0 0 9 -45.919 -3.35 3.4381e-10 26.822 0 0 10 -45.776 -3.35 -1.3321e-08 26.822 0 0 11 -44.764 -3.35 2.1362e-08 26.822 0 0 12 -44.228 -3.35 8.2932e-09 26.822 0 0 13 -43.675 -3.35 5.3661e-09 26.822 0 0 14 -42.967 -3.35 -7.7576e-09 26.822 0 0 15 -42.153 -3.35 1.3898e-08 26.822 0 0 16 -41.239 -3.35 -6.0162e-09 26.822 0 0
X=T1{:,2};
Y=T1{:,3};
U=T1{:,5};
V=T1{:,6};
W=T1{:,7};
X1=min(X):0.1:max(X);
Y1=min(Y):0.1:max(Y);
[xmesh, ymesh] = meshgrid(X1,Y1);
[x1 y1 u1]=griddata(X,Y,U,xmesh,ymesh);
[x1 y1 v1]=griddata(X,Y,V,xmesh,ymesh);
[x1 y1 w1]=griddata(X,Y,W,xmesh,ymesh);
res_vel=sqrt(u1.^2+v1.^2+w1.^2)
res_vel = 214×1601
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN 27.2200 27.2208 27.2217 27.2226 27.2228 27.2227 27.2227 27.2226 27.2225 27.2224 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2202 27.2211 27.2220 27.2228 27.2228 27.2227 27.2226 27.2226 27.2225 27.2224 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2205 27.2218 27.2227 27.2228 27.2228 27.2227 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2205 27.2218 27.2225 27.2228 27.2228 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2197 27.2204 27.2210 27.2217 27.2224 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2197 27.2204 27.2210 27.2217 27.2223 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2197 27.2204 27.2211 27.2217 27.2224 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2197 27.2204 27.2211 27.2217 27.2224 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206 NaN 27.2197 27.2204 27.2211 27.2218 27.2224 27.2228 27.2227 27.2226 27.2225 27.2225 27.2224 27.2223 27.2222 27.2221 27.2220 27.2219 27.2218 27.2217 27.2216 27.2215 27.2214 27.2213 27.2212 27.2211 27.2210 27.2209 27.2208 27.2207 27.2206
figure
[mc,hcf] = contourf(x1,y1,res_vel, 150);
Lvls = hcf.LevelList
Lvls = 1×151
0 0.3045 0.6090 0.9135 1.2180 1.5225 1.8270 2.1315 2.4360 2.7405 3.0450 3.3495 3.6539 3.9584 4.2629 4.5674 4.8719 5.1764 5.4809 5.7854 6.0899 6.3944 6.6989 7.0034 7.3079 7.6124 7.9169 8.2214 8.5259 8.8304
colormap(turbo)
shading interp
colorbar
figure
[mc,hcf] = contourf(x1,y1,res_vel, [25:0.1:28]);
Lvls = hcf.LevelList
Lvls = 1×31
25.0000 25.1000 25.2000 25.3000 25.4000 25.5000 25.6000 25.7000 25.8000 25.9000 26.0000 26.1000 26.2000 26.3000 26.4000 26.5000 26.6000 26.7000 26.8000 26.9000 27.0000 27.1000 27.2000 27.3000 27.4000 27.5000 27.6000 27.7000 27.8000 27.9000
colormap(turbo)
shading interp
colorbar
This bears no strong resemblance to the .png file, and I am not certain that is even possible.
Adding the ‘150’ to the first contourf argument tells it to create 150 contours, so with this, it added some detail.
Restricting the contours from 25 to 28 with increments of 0.1 creates the second contourf plot. This is likely as good as it is possible to get.
I encourage you to experiment with the options the contourf function offers to get something that more closely approximates the result you want.
.

Sign in to comment.

More Answers (0)

Categories

Find more on Contour Plots in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!