Is it possible to use wall functions in CFD simulation of flow around a cylinder or an airfoil when the objective is to find aerodynamic forces, vortex shedding behind and also separation point or we have to be careful about Y+ selection?
It depends on the flow regime (Reynolds and Mach Numbers), one of the best way is to use (RANS) with suitable turbulence model related to application and boundary conditions
Any wall-modelled BCs. (in both LES and RANS) are subject on a strong assumption on the boundary layer. In other words, you are somehow already prescribing the drag on the body. If you want to get the correct drag prediction you need to use a very fine grid close to the walls (at least 3-4 nodes within y+
If you are simulating the performance of an airfoil, you have to satisfy Y-plus distribution to be less than 1 over the airfoil. For this purpose, you can use a fine mesh with sufficient height of the first cell. The aerodynamics coefficients of the airfoil are so sensitive to y-plus value. In addition, the SST k-omega model has been proven to perform well in simulation of the attached flow of the airfoil.