Drag calculation FDM method (incompressible)
 

Hey,

1. What is the most appropriate method to calculate drag if using a FDM method (incompressible flow)?


2. So, if we have the force at the wall then how do you appropriately judge what area that force is acting on (i.e. for each node on the wall)?

1) you need to compute the surface integral of the stress at wall, you can use some trapezoidal rule for the integral and classic FD for the velocity gradient

2) the question is not clear to me ...force/area is a pressure, you have a force computed at each node?

I think you might also be able to compute the line integral around your control volume using the Reynolds Transport Theorem.

http://en.wikipedia.org/wiki/Reynolds_transport_theorem

OK, let me try to clarify.

In the FDM method we have no information on any area that the computed forces are acting on, right? We only know the node-wise distribution of forces per unit area. We do not have a continuous information of the boundaries, only node-wise information.

As an easy example, consider a stair-step inclined wall where the nodes naturally are distributed evenly. What information do we have about the real shape of the wall? In this easy example I would assume that all nodes contribute over an equally large area. --> But should that area be calculated from the known real shape of the wall or should it be calculated from say the area of the stair-step shape in the numerical grid?

A more difficult (perhaps) example is when the geometry is not well defined. Then it would seem appropriate to calculate the area each node is affecting on the wall.

More background:
I have computed 2d flow over a cylinder, complete with good results with regards to the pressure difference (before and after the cylinder). When it comes to drag and lift the results become good only if I use the known area (circumference) of the cylinder and distribute it evenly in all wall nodes. If I try to approximate the area however I fail to give a good answer.

but this problem is not due to FD ... it is a very general issue caused by the discrete world of computation!

when the geometry is defined only by a finite set of variable you can only approximate the continuous geometry by proper interpolation, for example splines...therefore the area depends on the type of reconstruction.

on the other hand, when you have a FD method of some order of accuracy, you can reconstruct the continuous field (for example the force) only by interpolation of degree congruent to order of the scheme.

For example if you use a second order FD then it makes sense to use linear interpolation for the variables.

However, for a well refined grid, I don't think that appears large discrepancy in the results

OK, but I think my problems come from the reconstruction of the surface.

If you look at the attached figure I have a schematic of how I deal with the reconstruction.

Green dots represent interior and black dots represent wall, I have no flag that states wheter a node is on the actual boundary between fluid and solid - instead I assume that this happens if a wall node is connected to an interior node along any cartesian axis (nodes 'a' and 'b' in the figure are not, so I never use Case 2 as reconstruction).

Node Area (length) calculation:

Case 2: All areas are equal to delta_x (=delta_y)

Case 1: If we look at the two nodes in the middle of the red line (the assumed shape of the surface) then I calculate the bottom one as delta_X*( sqrt(2)/2 + sqrt(2)/2 ) and the top one as delta_X*( sqrt(2)/2 + 1/2 )

The second option I consider (but have not tested) is to apply forces between each node using a mean of the nodal values, but I think this will give the same result as in Case 1 above.

Is Case 1 the correct way to do this (assuming the order is ok)?

but you are using structured cartesian grid for curvilinear geometries? that makes sense of your problem....

Yes, I am. This is of course a problem when it comes to boundary layers, but if we ignore that for now (low Re flow) I think it is possible to run cases with curved boundaries. The problem is not the fluid flow, but rather the post-processing of the forces. As I said, it is possible to get good results depending on which area reconstruction method is used (even on the relatively coarse meshes I have tested).

Could you please clarify what you mean by "that makes sense of your problem..."?

pressure distribution is highly sensitive to the geometry shape, independently from the Re number... dp/dn -> 0 in the boundary layer

OK. This is also part of the boundary condition where I have dp/dn=0.

By highly sensitive, in what manner do you mean? Do you mean it is impossible to get a good solution when using a non-boundary conforming mesh?

just think about this approximation: dp/dn = 0 in the BL means you have on the body the same pressure distribution of the external (inviscid-like) flow. But the pressure in that region depends on the velocity around the body (approximate with Bernouilli in the external region) and the velocity distribution depends on the geometry (potential-like flow). Therefore, if you in computation represents a curvilinear surface by step-wise grid you can have strong errors.
The only remedy is to use a ver very refined geometry.

Otherwise is better to adopt a different approach such as the immersed boundary method

The error comes from the inability to put computational nodes at the curved boundary, right?

The stair-step geometry does not really exist in the FDM, that defines nodes as in the previous figure, it only exist as a post-processing artifact. If you consider the two cases in the previous figure, they will produce identical flow over the obstacle, but the forces calculated - after the simulation - will be different.

I agree that a fine mesh is needed when comparing to an unstructured mesh that can place nodes directly at the wall. This might be acceptable in certain cases though.

After a couple of months focusing on other matters than CFD on my spare time I got the opportunity to sit down with my structured FDM code this weekend. I would like to revisit this old question and continue with:

If the code manages to produce the correct pressure difference for flow over a cylinder (pressure before cylinder and after) then how can the structured Cartesian grid approach be so wrong? The pressure behind the cylinder should be very dependent on the geometry. I still believe the drag calculation itself is the problem (it is not extremely far off but still not acceptable at Cd=5.69 for Re=20).

What do you think?