Your approach is absolutely acceptable so long as you realize its limits. The reason you are not seeing a method like this in a text book is because it is a mix of numerical methods and conventional empirical methods that is not particularly common.
The reason it is not common to approach a problem this way is because the natural convection heat transfer coefficient (htc) is a bulk property. It is based on the overall geometry and is going to give you a prediction of the heat transfer over the entire surface of the object. Thus the htc is based on the overall system geometry and properties and not the local geometry and properties right around each node. My hope is this image will make that clear:
Here we have a hot plate in some cold fluid. We have broken the hot plate into 6 nodes to estimate the surface temperature, but we can't calculate the htc based on only each nodes temperature because that is not the way natural convection works.
The most robust way to solve this is to solve the natural convection numerically as well. This is what a commercial CFD package like comsol would do. 108
However, this would require that you model the fluid numerically which as you say above would get very complicated.
The other approach you could take is to continue using the bulk htc but recalculate it each time you iterate through your internal numerics. Then you would come up with a bulk htc and you could calculate the heat loss from each node via Q=dA*htc*(T_i-T_bulk). When you are running your solver you just need to recalculate the htc based on a new average surface temperature each iteration.
There are other compromises you could make here to increase the accuracy of your model, but the overarching theme is to know the limits of each compromise you make. Using a bulk htc is an ok move if your fluid has a high Cp and thus will not have a big dT as it flows over the cube. It is also an ok move if the geometry is a cube or a sphere where the water is not flowing over the part for a long time. This would be a bad approach if your fluid had a low Cp and thus had a big dT as it flowed over your part or if the part had a very long flow path as it does not capture the change in duty caused by the change in dT between the first and last node the water flows past.
All in all, I think the model is simple enough you should just write it and see what the results are. You can then start messing with the temperature you use for the bulk temperature to try to capture the water's temperature better. 108