I tried it a bit different with a somewhat disturbing result

the undefined inner integral

j0 = \[Integral](Sqrt[196 - y^2 - x^2]) \[DifferentialD]x

insert the the upper bound yields a division by zero

j12 = j0 /. x -> Sqrt[14^2 - y^2] // Simplify

So try a limit, which seems to work

j12 = Limit[j0 /. x -> u, u -> Sqrt[14^2 - y^2]]

The lower limit is no problem

j11 = j0 /. x -> a + y/e

so the resulting inner integral should be

j2 = j12 - j11

Specify that for a = 1 and e = 1

j211 = j2 /. a -> 1 /. e -> 1

The upper bound for the outer integral is (and directly specified for a = e = 1 )

uL = (-a e + e Sqrt[196 + 196 e^2 - a^2 e^2])/(1 + e^2);
uL11 = uL /. {a -> 1, e -> 1}

But the integration gives something different from f [ 1, 1 ] given above:

In[11]:= NIntegrate[j211, {y, 0, uL11}]

Out[11]= -1843.8

Where is the flaw?

I think that is not true, especially when taking into account that Sqrt[ 196 -x^2- y^2 ] is not a complex function. The flaw is the strict confidence in the Limit Procedure, which gave

1/4 \[Pi] (-196 + y^2)

as value of the inner integral at the upper bound. When I had a close look at this limit

Limit[j0 /. x -> uu Sqrt[14^2 - y^2], uu -> 1]

it turned out that the negative value must be used. So writing

j2 = -j12 - j11 // FullSimplify

instead of the j2 given above and

jj2[a_, e_, y_] := Evaluate[j2]

and

ff[a_, e_] :=  NIntegrate[   jj2[a, e, y], {y,  0, ((-a)*e + e*Sqrt[196 + 196*e^2 - a^2*e^2])/(1 + e^2)}]

which when plotting gives the same picture as the other pics above

Plot3D[ff[a, e], {a, 0.01, 2}, {e, 0.01, 2}, ColorFunction -> "SolarColors"]

What is the learning? Mathematica is really a great software, but you can't always believe everything.

The flaw was that I didn't check if this limiting process gave a correct answer.

Adding assumptions seems to help. Below is one possibility. It took many minutes to evaluate though.

In[1]:= Integrate[Sqrt[196 - x^2 - 
       y^2], {y, 0, 
     ((-a)*e + e*Sqrt[196 + 
                196*e^2 - a^2*e^2])/
       (1 + e^2)}, {x, a + y/e, 
     Sqrt[14^2 - y^2]}, Assumptions -> {0 < a < 14, e > 0}]

(* Out[1]= ConditionalExpression[(1/(
 12 (1 + e^2)^(
  3/2)))(2 a^2 e Sqrt[(196 - a^2) (1 + e^2)] - 588 a e \[Pi] - 
   588 a e^3 \[Pi] + a^3 e^3 \[Pi] - 5488 Sqrt[1 + e^2] \[Pi] - 
   5488 e^2 Sqrt[1 + e^2] \[Pi] + 5488 I Sqrt[1 + e^2] Log[2] + 
   5488 I e^2 Sqrt[1 + e^2] Log[2] - 2744 I Sqrt[1 + e^2] Log[3] - 
   2744 I e^2 Sqrt[1 + e^2] Log[3] + 5488 I Sqrt[1 + e^2] Log[49] + 
   5488 I e^2 Sqrt[1 + e^2] Log[49] - 
   2744 I Sqrt[1 + e^2] Log[9604/3] - 
   2744 I e^2 Sqrt[1 + e^2] Log[9604/3] + 
   I a e (-588 + (-588 + a^2) e^2) Log[1 + e^2] + 
   I a e (-588 + (-588 + a^2) e^2) Log[(196 - (-196 + a^2) e^2)/(
     1 + e^2)] - 
   5488 I Sqrt[1 + e^2]
     Log[-14 (Sqrt[196 - a^2] + 14 I e) + I a^2 e + 
      a (14 I - Sqrt[196 - a^2] e)] - 
   5488 I e^2 Sqrt[1 + e^2]
     Log[-14 (Sqrt[196 - a^2] + 14 I e) + I a^2 e + 
      a (14 I - Sqrt[196 - a^2] e)] + 
   5488 I Sqrt[1 + e^2]
     Log[14 (Sqrt[196 - a^2] - 14 I e) + I a^2 e - 
      a (14 I + Sqrt[196 - a^2] e)] + 
   5488 I e^2 Sqrt[1 + e^2]
     Log[14 (Sqrt[196 - a^2] - 14 I e) + I a^2 e - 
      a (14 I + Sqrt[196 - a^2] e)] + 
   1176 I a e Log[-I a + Sqrt[-(-196 + a^2) (1 + e^2)]] + 
   1176 I a e^3 Log[-I a + Sqrt[-(-196 + a^2) (1 + e^2)]] - 
   2 I a^3 e^3 Log[-I a + Sqrt[-(-196 + a^2) (1 + e^2)]] - 
   5488 I Sqrt[1 + e^2]
     Log[14 + a e + 14 e^2 - e Sqrt[196 - (-196 + a^2) e^2]] - 
   5488 I e^2 Sqrt[1 + e^2]
     Log[14 + a e + 14 e^2 - e Sqrt[196 - (-196 + a^2) e^2]] + 
   5488 I Sqrt[1 + e^2]
     Log[14 - a e + 14 e^2 + e Sqrt[196 - (-196 + a^2) e^2]] + 
   5488 I e^2 Sqrt[1 + e^2]
     Log[14 - a e + 14 e^2 + e Sqrt[196 - (-196 + a^2) e^2]] - 
   5488 I Sqrt[1 + e^2]
     Log[-(-196 + a^2) e^3 + 14 Sqrt[196 - (-196 + a^2) e^2] + 
      14 e^2 Sqrt[196 - (-196 + a^2) e^2] + 
      e (196 - a Sqrt[196 - (-196 + a^2) e^2])] - 
   5488 I e^2 Sqrt[1 + e^2]
     Log[-(-196 + a^2) e^3 + 14 Sqrt[196 - (-196 + a^2) e^2] + 
      14 e^2 Sqrt[196 - (-196 + a^2) e^2] + 
      e (196 - a Sqrt[196 - (-196 + a^2) e^2])] + 
   5488 I Sqrt[1 + e^2]
     Log[(-196 + a^2) e^3 + 14 Sqrt[196 - (-196 + a^2) e^2] + 
      14 e^2 Sqrt[196 - (-196 + a^2) e^2] + 
      e (-196 + a Sqrt[196 - (-196 + a^2) e^2])] + 
   5488 I e^2 Sqrt[1 + e^2]
     Log[(-196 + a^2) e^3 + 14 Sqrt[196 - (-196 + a^2) e^2] + 
      14 e^2 Sqrt[196 - (-196 + a^2) e^2] + 
      e (-196 + a Sqrt[196 - (-196 + a^2) e^2])]), 
 196 e <= a (14 + a e)] *)

Among this double integral, "a" and "e" are assumed to be unknown parameters. I want to calculate this integral when "a" or "e" is set as a constant. And if possible, I want to draw the 3D plot of this integral with unknown "a" and "e". Thanks very much.

I speed up computation at factor 10.

 n = 1/4;(*Decrease the value to smooth a plot,the calculation time will be longer *)
 ListPlot3D[Partition[Flatten[Table[{a, e, f[a, e]}, {a, -2.01, 2, n}, {e, -2.01, 2, n}]],
 3], ColorFunction -> "SolarColors", PlotLegends -> Automatic]

regards, Mariusz

Hi

Yours $r$ variable need a value. You forgot to add it.

NIntegrate needs to be all variable constans a numeric. I'm add Method -> "LocalAdaptive" possible Integral may be singular at some values.

r=1;(* you may change this value*)
f[a_, e_, r_] := 196 Pi a - 2 (NIntegrate[Sqrt[196 - x^2], {y, 0, r Tan[e]}, {x, y/Tan[e], r}, 
Method -> "LocalAdaptive"] - (e r^3)/3) + 2 Re[NIntegrate[Sqrt[196 - x^2], {y, 0, r Tan[e] - a}, {x, (a + y)/Tan[e], r}, 
Method -> "LocalAdaptive"] - NIntegrate[Sqrt[-x^2 - y^2 + 196], {y, 0, -a Cos[e]^2 + (Sqrt[2] Sin[e])/2 Sqrt[
392 - a^2 - (a^2) Cos[2 e]]}, {x, (a + y)/Tan[e], Sqrt[14^2 - y^2]}, Method -> "LocalAdaptive"]]

n = 1/3;(*Decrease the value to smooth a plot,the calculation timewill be longer*)
ListPlot3D[Partition[Flatten[Table[{a, e, f[a, e, r]}, {a, -2.01, 2, n}, {e, -2.01, 2, n}]], 3], ColorFunction -> "SolarColors", 
PlotLegends -> Automatic]

For more information read Documentation Center it's all there,everything you need.

Regards, Mariusz