【CFD理论】对流项-02

知识回顾
- 边界问题
properties of discretisation schemes
- conservation 守恒性
- Boundedness 有界性
- - example
  - case2- wiggle-对流扩散
- Transportiveness 输运性
- Accuracy 精度
- - QUICK格式3阶

知识回顾

aE=D−F2a_E=D-\frac{F}{2}aE=D−2F
aW=D+F2a_W=D+\frac{F}{2}aW=D+2F

边界问题

Dirichlet
Newmann

properties of discretisation schemes

conservation 守恒性

用纯扩散性方程分析：
(Γϕ2−ϕ1δx−qA)+(Γϕ3−ϕ2δx−Γϕ2−ϕ1δx)+(Γϕ4−ϕ3δx−Γϕ3−ϕ2δx)+(qB−Γϕ4−ϕ3δx)=qB−qA(\Gamma\frac{\phi_2-\phi_1}{\delta x}-q_A)+(\Gamma\frac{\phi_3-\phi_2}{\delta x}-\Gamma\frac{\phi_2-\phi_1}{\delta x})+(\Gamma\frac{\phi_4-\phi_3}{\delta x}-\Gamma\frac{\phi_3-\phi_2}{\delta x})+(q_B-\Gamma\frac{\phi_4-\phi_3}{\delta x})=q_B-q_A(Γδxϕ2−ϕ1−qA)+(Γδxϕ3−ϕ2−Γδxϕ2−ϕ1)+(Γδxϕ4−ϕ3−Γδxϕ3−ϕ2)+(qB−Γδxϕ4−ϕ3)=qB−qA

Boundedness 有界性

对角占优 diagonally dominant
ϕ\phiϕ应该在neighboring nods的值范围内
离散方程的所有系数应具有相同的符号
不过不符合以上条件，就可能越界产生“wiggles”现象

example

一维稳态扩散方程
ddx(kdTdx)+S=0\frac{d}{dx}(k\frac{dT}{dx})+S=0dxd(kdxdT)+S=0
(kATR−TPdPR)−(kATP−TLdLP)+S‾V=0(kA\frac{T_R-T_P}{d_{PR}})-(kA\frac{T_P-T_L}{d_{LP}})+\overline S V=0(kAdPRTR−TP)−(kAdLPTP−TL)+SV=0

内部网格：
(kAdTdx)r−(kAdTdx)l+S‾V=0(kA\frac{dT}{dx})_r-(kA\frac{dT}{dx})_l+\overline S V=0(kAdxdT)r−(kAdxdT)l+SV=0
TP(DlAl+DrAr+0)=TL(DlAl)+TR(DrAr)+S‾V=0T_P(D_lA_l+D_rA_r+0)=T_L(D_lA_l)+T_R(D_rA_r)+\overline SV=0TP(DlAl+DrAr+0)=TL(DlAl)+TR(DrAr)+SV=0
边界网格：
(kAdTdx)r−(kAdTdx)l+S‾V=0(kA\frac{dT}{dx})_r-(kA\frac{dT}{dx})_l+\overline S V=0(kAdxdT)r−(kAdxdT)l+SV=0
left
(dTdx)l=TP−TAdLP/2(\frac{dT}{dx})_l=\frac{T_P-T_A}{d_{LP/2}}(dxdT)l=dLP/2TP−TA
(kATR−TPdPR)−(klTP−TAdLP/2)+V‾=0(kA\frac{T_R-T_P}{d_{PR}})-(kl\frac{T_P-T_A}{d_{LP}/2})+\overline V=0(kAdPRTR−TP)−(kldLP/2TP−TA)+V=0
TP(0+DA+2DA)=TL(0)+TR(DA)+TA(2DA)+S‾VT_P(0+DA+2DA)=T_L(0)+T_R(DA)+T_A(2DA)+\overline SVTP(0+DA+2DA)=TL(0)+TR(DA)+TA(2DA)+SV
right
TP(DA+0+2DA)=TL(DlAl)+TR(0)+TB(2DA)+S‾VT_P(DA+0+2DA)=T_L(D_lA_l)+T_R(0)+T_B(2DA)+\overline SVTP(DA+0+2DA)=TL(DlAl)+TR(0)+TB(2DA)+SV

分割网格
Lcell=LN=1mL_cell=\frac{L}{N}=1mLcell=NL=1m
分配材料属性
DA=10∗0.1/1=10[W/K]DA=10*0.1/1=10[W/K]DA=10∗0.1/1=10[W/K]
S‾V=1000∗0.1∗1=100[W]\overline S V=1000*0.1*1=100[W]SV=1000∗0.1∗1=100[W]
矩阵系数

A matrix:30.0  -10.0    0.0    0.0    0.0 -10.0   20.0  -10.0    0.0    0.0 0.0  -10.0   20.0  -10.0    0.0 0.0    0.0  -10.0   20.0  -10.0 0.0    0.0    0.0  -10.0   30.0 B vector:
2100.0100.0100.0100.0
4100.0

求解矩阵
T=A/BT=A/BT=A/B

case2- wiggle-对流扩散

close all; clear; clc
%--------------------------------------------------------------------------
% User Inputs
%--------------------------------------------------------------------------% Thermal Conductivity of the fluid (W/mK)
cond = 0.1;% Specific heat capacity of the bar (J/Kg K)
cp = 1;% Cross-sectional Area of the channel (m2)
area = 1;% Length of the channel (m)
barLength = 1;% Number of cells in the mesh
nCells = 5;% Temperature at the left hand side of the channel (deg C)
tempLeft = 1;% Temperature at the right hand side of the channel (deg C)
tempRight = 0;% Heat source per unit volume (W/m3)
heatSourcePerVol = 0;% Flow velocity (m/s)
flowVelocity = 2.5;% Fluid Density (kg/m3)
fluidDensity = 1.0;% Plot the data? ('true' / 'false')
plotOutput = 'true';% Print the set up data? (table of coefficients and matrices)
printSetup = 'true';% Print the solution output (the temperature vector)
printSolution = 'true';% =========================================================================
% Code Begins Here
% =========================================================================
%--------------------------------------------------------------------------
% 1. Print messages
%--------------------------------------------------------------------------
fprintf('================================================\n');
fprintf('\n');
fprintf('   solve1DConvectionDiffusionEquationUpwind.m\n')
fprintf('\n')
fprintf(' - Fluid Mechanics 101\n')
fprintf(' - Author: Dr. Aidan Wimshurst\n')
fprintf(' - Contact: FluidMechanics101@gmail.com\n')
fprintf('\n')
fprintf(' [Exercise 1: Chapter 4]\n')
fprintf('================================================\n')%--------------------------------------------------------------------------
% 2. Create the mesh of cells
%--------------------------------------------------------------------------fprintf(' Creating Mesh ...\n');
fprintf('------------------------------------------------\n');% Start by calculating the coordinates of the cell faces
xFaces = linspace(0, barLength, nCells+1);% Calculate the coordinates of the cell centroids
xCentroids = 0.5*(xFaces(2:end) + xFaces(1:end-1));% Calculate the length of each cell
cellLength = xFaces(2:end) - xFaces(1:end-1);% Calculate the distance between cell centroids
dCentroids = xCentroids(2:end) - xCentroids(1:end-1);% For the boundary cell on the left, the distance is double the distance
% from the cell centroid to the boundary face
dLeft = 2*(xCentroids(1) - xFaces(1));% For the boundary cell on the right, the distance is double the distance from
% the cell centroid to the boundary cell face
dRight = 2*(xFaces(end) - xCentroids(end));% Append these to the vector of distances
dCentroids = [dLeft, dCentroids, dRight];% Compute the cell volume
cellVolumes = cellLength*area;%--------------------------------------------------------------------------
% 2.0 Calculate the Matrix Coefficients
%--------------------------------------------------------------------------fprintf(' Calculating Matrix Coefficients\n');
fprintf('------------------------------------------------\n');% Diffusive flux per unit area
DA = area*cond./dCentroids;% Convective flux
velocityVector = flowVelocity*ones(1, nCells+1);
densityVector = fluidDensity*ones(1, nCells+1);
F = densityVector.*velocityVector*area*cp;% Evaluate the Peclet Number
Pe = F./DA;% Calculate the source term Sp
Sp = zeros(1, nCells);% Assign sources to the left and right hand boundaries
Sp(1) = -1.0*(2.0*DA(1) +F(1));
Sp(end) = -1.0*(2.0*DA(end) -F(end));% Calculate the source term Su
Su = heatSourcePerVol*cellVolumes;% Assign sources to the left and right hand boundaries
Su(1) = 3.5*tempLeft;
Su(end) = -3.5*tempRight;% Left Coefficient (aL)
aL = DA(1:end-1) + (F(1:end-1)/2);%Right Coefficient (aR)
aR = DA(2:end) + (-F(2:end)/2);% Set the first element of aL to zero. It is a boundary face
aL(1) = 0;% Set the last element of aR to zero. It is a boundary face
aR(end) = 0;% Create the central coefficients
aP = aL + aR  - Sp;%--------------------------------------------------------------------------
% 3.0 Create the matrices
%--------------------------------------------------------------------------fprintf(' Assembling Matrices\n');
fprintf('------------------------------------------------\n');% Start by creating an empty A matrix and an empty B matrix
Amatrix = zeros(nCells, nCells);
BVector = zeros(nCells, 1);% Populate the matrix one row at a time (i.e one cell at a time)
%
% NOTE: this method is deliberately inefficient for this problem
% but it is useful for learning purposes. We could populate the
% diagonals and the off-diagonals directly.for i = 1:nCells% Do the A matrix Coefficients% Left boundary Cellif (i == 1)Amatrix(i, i) = aP(i);Amatrix(i, i+1) = -1.0*aR(i);% Right Boundary Cellelseif(i == nCells)Amatrix(i, i-1) = -1.0*aL(i);Amatrix(i, i) = aP(i);% Interior CellselseAmatrix(i, i-1) = -1.0*aL(i);Amatrix(i, i) = aP(i);Amatrix(i, i+1) = -1.0*aR(i);end% Do the B matrix coefficientsBVector(i) = Su(i);
end%--------------------------------------------------------------------------
% 4.0 Print the setup
%--------------------------------------------------------------------------if strcmp(printSetup, 'true')fprintf(' Summary: Set Up\n')fprintf('------------------------------------------------\n')fprintf('aL: ')for i = 1:nCellsfprintf('%6.1f ', aL(i)); endfprintf('\n')fprintf('aR: ')for i = 1:nCellsfprintf('%6.1f ', aR(i)); endfprintf('\n')fprintf('aP: ')for i = 1:nCellsfprintf('%6.1f ', aP(i)); endfprintf('\n')fprintf('Sp: ')for i = 1:nCellsfprintf('%6.1f ', Sp(i)); endfprintf('\n')fprintf('Su: ')for i = 1:nCellsfprintf('%6.1f ', Su(i)); endfprintf('\n')fprintf('Pe: ')for i = 1:nCells+1fprintf('%6.1f ', Pe(i)); endfprintf('\n')fprintf('A matrix:\n')for i = 1:nCellsfor j = 1:nCellsfprintf('%6.1f ', Amatrix(i, j)); endfprintf('\n')endfprintf('B vector:\n')for i = 1:nCellsfprintf('%6.1f\n', BVector(i)); endfprintf('\n')fprintf('------------------------------------------------\n')
end%--------------------------------------------------------------------------
% 5.0 Solve the matrices
%--------------------------------------------------------------------------fprintf(' Solving ...\n')
fprintf('------------------------------------------------\n')% Use MATLAB's default linear algebra solver (AX = B)
Tvector = Amatrix \ BVector;fprintf(' Equations Solved.\n')
fprintf('------------------------------------------------\n')%--------------------------------------------------------------------------
% 6.0 Print the Solution
%--------------------------------------------------------------------------
if strcmp(printSolution, 'true')fprintf(' Solution: Temperature Vector\n');fprintf('------------------------------------------------\n');fprintf('T vector:\n')for i = 1:nCellsfprintf('%10.6f\n', Tvector(i)); endfprintf('\n')fprintf('------------------------------------------------\n')
end%--------------------------------------------------------------------------
% 7.0 Plot the solution
%--------------------------------------------------------------------------% Plot the data if desired
if strcmp(plotOutput,'true')fprintf(' Plotting ...\n')fprintf('------------------------------------------------\n')% Append the boundary temperature values to the vector, so we can% plot the complete solutionxPlotting = [xFaces(1), xCentroids, xFaces(end)];temperaturePlotting = [tempLeft; Tvector; tempRight];% Figure Size ParametersfigSizeXcm = 8.5;aspectRatio = 4.0/3.0;figSizeYcm = figSizeXcm / aspectRatio;figSizeXpixels = figSizeXcm * 37.79527559055118;figSizeYpixels = figSizeYcm * 37.79527559055118;% Figure font, text and line widthsfontSize = 10;fontChoice = 'Arial';lineWidth = 1.0;markerSize = 4;% Colours for the line plotscolour1 = 'black';colour2 = [0.0, 0.129, 0.2784];% Adjust anchor for where figure shows up on screen (px)x0 = 100;y0 = 100;fig1 = figure('Name', '1D Diffusion');box on;hold on;plot(xPlotting, temperaturePlotting, 'k-o', 'Color', colour2, ...'MarkerSize', markerSize, 'linewidth', lineWidth + 0.2);hold off;xlabel('x [m]', 'FontSize', fontSize)ylabel('T [^{\circ} C]', 'FontSize', fontSize)grid on;set(gca, 'fontsize', fontSize);set(gca, 'linewidth', lineWidth);set(gcf, 'position', [x0, y0, figSizeXpixels, figSizeYpixels]);legend({'CFD'}, 'NumColumns', 1, 'Location', 'Best')end%==========================================================================
% END OF FILE
%==========================================================================

结果：

================================================Creating Mesh ...
------------------------------------------------Calculating Matrix Coefficients
------------------------------------------------Assembling Matrices
------------------------------------------------Summary: Set Up
------------------------------------------------
aL:    0.0    1.8    1.8    1.8    1.8
aR:   -0.8   -0.7   -0.7   -0.8    0.0
aP:    2.8    1.0    1.0    1.0    0.3
Sp:   -3.5    0.0    0.0    0.0    1.5
Su:    3.5    0.0    0.0    0.0   -0.0
Pe:    5.0    5.0    5.0    5.0    5.0    5.0
A matrix:2.8    0.8    0.0    0.0    0.0 -1.8    1.0    0.7    0.0    0.0 0.0   -1.8    1.0    0.7    0.0 0.0    0.0   -1.8    1.0    0.8 0.0    0.0    0.0   -1.8    0.3
B vector:3.50.00.00.0-0.0------------------------------------------------Solving ...
------------------------------------------------Equations Solved.
------------------------------------------------Solution: Temperature Vector
------------------------------------------------
T vector:1.0356300.8693551.2573310.3520532.464370------------------------------------------------Plotting ...
------------------------------------------------

CD方法计算，Pe数大（5）造成越界

参考：

book-An Introduction to Computational Fluid Dynamics, The Finite Volume Method

================================================solve1DConvectionDiffusionEquationUpwind.m- Fluid Mechanics 101- Author: Dr. Aidan Wimshurst- Contact: FluidMechanics101@gmail.com[Exercise 1: Chapter 4]

Transportiveness 输运性

Accuracy 精度

QUICK格式3阶

ϕe=38ϕE+68ϕP−18ϕW\phi_e=\frac{3}{8}\phi_E+\frac{6}{8}\phi_P- \frac{1}{8} \phi_Wϕe=83ϕE+86ϕP−81ϕW

用了三个点，三阶格式

迎风格式，非常稳定的格式，不会越界，守恒，一阶格式，会带来smear现象。
CD，可能会导致越界，二阶
二阶迎风格式
QUICK，三阶格式