Hi, I am solving a convection diffusion equation. There are four boundary conditions. I know the right and bottom boundary have no flux condition. Also, the left side looks irregular which denotes the ice water boundary. This part consists of more than 60 segments with each being a boundary line segment. Now, by guessing the temperature distribution at the top boundary, I was able to solve the equation and thus getting the boundary condition everywhere inside the domain. Now, I want to compare the flux at the left(irregular) boundary and the top boundary. How should I calculate the flux? Would the following code be right? Note that the domain is cylindrical, so I need to integrate using cylindrical coordinates. The code is following: //Create an array from 4
int[int] leftlabels(61); // Array to store values from 4
for (int i = 0; i < 61; i++) {
leftlabels[i] = i + 4; // Fill the array with values from 4
}
uvel = Uvel; // velocity in r direction
wvel = Wvel;
solve thermic(T, v)
= int2d(Mesh)(
+ k*(
dx(T) * dx(v)
+ dy(T) * dy(v)
+ (1/x) * dx(T) * v
)
- (+ uvel * dx(T) * v
+ wvel * dy(T) * v)
)
// Neumann Conditions
+ int1d(Mesh, bottomside) (0v)
+ int1d(Mesh, rightside) (0v)
// Dirchlet Conditions
+ on(watersurface, T=Tsurface) //x in func usurface depends on the watersurface coordinates
+ on(leftlabels, T=Tleft);
real totalfluxleft = 0.;
real totalfluxtop = 0.;
// Integrate heat flux over each boundary to compute total heat transfer (This is the air water boundary)
totalfluxtop = int1d(Mesh, 3)(-k * (dx(T) * N.x + dy(T) * N.y) * x); // Label 3 for the top boundary (* deltaR is deleted)
// Loop through left boundary sub-boundaries
for (int label = 4; label <= 64; label++) {
real flux = int1d(Mesh, label)(-k * (dx(T) * N.x + dy(T) * N.y) * x);
totalfluxleft += flux;
}
totalfluxleft = totalfluxleft * 2 * pi;
totalfluxtop = totalfluxtop * 2 * pi;
// Output the total flux values
cout << "Total Heat Flux on Left Boundary (Ice-Water): " << totalfluxleft << endl;
cout << "Total Heat Flux on Top Boundary (Water-Air): " << totalfluxtop << endl;
Also, the domain looks like the following graph:
In other words, I would like to know how does boundary flux in freeFEM work? How do we obtain it, is there a short cut?