Based on my understanding of what OP wants, no, your code does not work with arbitrary functions. Let's say that we want to swap fC3 with a new function that takes two arguments. The way your code is written, hF1 needs to know that, and so it also needs to change. Or let's say we wanted to replace fC3 with a new function that takes 5 arguments. Now your code is broken completely, because hF1 explicitly expects all of the intermediate functions to take at most 4 arguments. Your code has fairly tightly coupled the "shape" of the hF1 function (or any alternative we might want to use) to the "shape" of the intermediate fA1 (et al) functions. Now, maybe that's suffient for the OPs case, but I got the impression that something more general was wanted.

Thank you so much. I will work on your solution and update this post.

not really sure that you are actually wanting to find the functions. The code looks for a vector of 4 random numbers that when evaluated by a function stay within a range.

The function you are evaluating is

fA1[w1_, w2_, w3_, w4_] := 10 + w1 + w2 + w3 + w4
fB2[r1_, r2_, r3_, r4_] := 1 + r1 + r2 + r3 + r4
fC3[t1_] := 5*t1
hF1[fA1_, fB2_, fC3_] := (2*fA1)*(3*fB2)*(4*fC3)

FullSimplify[
 hF1[fA1[d1 + d2, d3, d1, d2], fB2[d1 + d2, d3, d1, d2], fC3[2*d1]]]

which evaluates to

240 d1 (1 + 2 d1 + 2 d2 + d3) (2 (5 + d1 + d2) + d3)

so:

f[{d1_, d2_, d3_, d4_}] := 
 240 d1 (1 + 2 d1 + 2 d2 + d3) (2 (5 + d1 + d2) + d3)

the input parameters are

n = 100
{min, max} = {1, 10^6};
ranges = {{1, 2}, {1, 2}, {1, 2}, {1, 2}};

so this would solve the problem in one line.

Column[{#, f[#]} & /@ Select[Transpose[RandomReal[#, n] & /@ ranges], min < f[#] < max &]]

or step by step

(*gerenrate the random vectors where each row of ds = {d1, d2, d3, d4}*)
ds = Transpose[RandomReal[#, n] & /@ ranges]
(*select the values for which each row of ds remains in the limits*)
sel = Select[ds, min < f[#] < max &]
(*show the output show the output, where it recomputes f again, which is not needed*)
Column[{#, f[#]} & /@ sel]

Let's start from the inside and work our way out.

LowerLimit <= hF1[fA1, fB2, fC3] <= UpperLimit

What you actually passed in for argument fA1 was

fA1[D1 + D2, D3, D1, D2]

When you start evaluating the Module, D1/2/3/4 have no values assigned to them. But you're also passing in these symbols when you call GT5. The argument names are x1/2/3/4. So the lines like

x1 = RandomReal[{p1min, p1max}]

actually get interpreted as

D1 = RandomReal[{p1min, p1max}]

That's why fA1[D1+D2,D3,D1,D2] actually resolves to some actual value.

Taking a tangent, this explains why you get those messages the second time around. The second time around, D1 has a numeric value. Thus, the argument x1 gets a numeric value. And so, the line

x1 = RandomReal[{p1min, p1max}]

gets interpreted as

1.6817646433788693 = RandomReal[{p1min, p1max}]

(or whatever, it's random). You can't assign a new value to an actual number.

The main problem you're having is you're duplicating information. The first argument you pass to GT5 is fA1[D1+D2,D3,D1,D2], but then your fourth argument is D1. You're basically passing in D1 several times. Maybe you were trying to parameterize GT5 to allow arbitrary functions being passed in, but let's defer that and just assume that you always want to compute with fA1, fB2, and fC3. This allows us to just apply those functions in the body of GT5. Specifically, the line we started looking at could be written like this:

LowerLimit <= hF1[fA1[x1 + x2, x3, x1, x2], fB2[x1 + x2, x3, x1, x2], fC3[2*x1]] <= UpperLimit

This means we can simplify the signature of GT5 like this:

GT5[x1_, x1min_, x1max_, x2_, x2min_, x2max_, x3_, x3min_, x3max_, x4_, x4min_, x4max_, NN_] := ...

And we can go further, because x1/2/3/4 were just dummy arguments. All we really need are the 4 pairs of min/max. We can use a local variable for each--we don't need to pass that in. Furthermore, the variables p1/2/3/4 aren't doing anything--let's eliminate them. So, at this point, we have something like this:

GT5[x1min_, x1max_, x2min_, x2max_, x3min_, x3max_, x4min_, x4max_, NN_] :=
  Module[
    {x1, x2, x3, x4,
     n, selectedData, LowerLimit, UpperLimit},
    LowerLimit = 1;
    UpperLimit = 10000000;
    selectedData = {};
    n = 1;
    Do[
      x1 = RandomReal[{x1min, x1max}];
      x2 = RandomReal[{x2min, x2max}];
      x3 = RandomReal[{x3min, x3max}];
      x4 = RandomReal[{x4min, x4max}]; 
      If[
        LowerLimit <= hF1[fA1[x1+x2,x3,x1,x2], fB2[x1+x2,x3,x1,x2], fC3[2*x1]] <= UpperLimit,
        selectedData = AppendTo[selectedData, {x1,x2,x3,x4,hF1[fA1[x1+x2,x3,x1,x2], fB2[x1+x2,x3,x1,x2], fC3[2*x1]]}], 
        0],
      {n, NN}];
    selectedData]

There is still more work that should be done. We're calculating some things twice, we still have some superfluous local variables, and we probably want to make the upper/lower limits arguments instead of magic local values. But hopefully the above is enough to get you started.

When you define this:

hF1[fA1_,fB2_,fC3_]:= (2*fA1)*(3*fB2)*(4*fC3)

Are you using the argument names fA1 and so forth as just a mnemonic device, or are you expecting some actual computational link to the functions fA1 etc that you previous defined?