Volume of a bead with square hole- Program approach?

[Albert Chan]

Here is a simplied version of Werner's code, assumed diameter = 1

I like dimensionless version because you can get a sense of results.
Example, with diameter of 1, sphere volume = Pi/6 ≈ 0.523598775598 ≈ ½
You can sense how big the hole is ...

Also, A, B is limited to inside a unit circle.

Code:

00 { 83-Byte Prgm }
01▸LBL "RB1"
02 MVAR "A"
03 MVAR "B"
04 MVAR "V"
05 1
06 RCL "A"
07 X↑2
08 -
09 RCL "B"
10 X↑2
11 -
12 SQRT
13 LSTO "K"
14 RCL "A"
15 RCL× "B"
16 STO× ST Y
17 RCL÷ "K"
18 ATAN
19 -
20 3
21 ÷
22 RCL "A"
23 RCL "B"
24 XEQ 10
25 RCL- "V"
26 RCL "B"
27 RCL "A"
28▸LBL 10
29 RCL÷ "K"
30 ATAN
31 X<>Y
32 STO× ST Y
33 X↑2
34 3
35 -
36 6
37 ÷
38 ×
39 -
40 END

Example, to get sphere - hole, diameter = 20, hole = 2x2

Solve with A=0.1, B=0.1, we have hole = 0.00996658845699, about 2% of sphere

Pi/6 - V = 0.513632187141

Scale up by D^3 => sphere - hole = 0.513632187141 * 20^3 = 4109.05749713

Update:
Perhaps we should compare volume against maximum possible volume
Solve with sphere - 4 caps, or simply solve with A = B = √(0.5), we get

V(max) ≈ 0.402001836518

So, above hole ≈ 2½ % of maximum possible volume.

[Albert Chan]

(06-11-2020 06:53 PM)Albert Chan Wrote:  Solve with A=0.1, B=0.1, we have hole = 0.00996658845699, about 2% of sphere

For small square hole, we can estimate with simple formula
For unit diameter (d=1), with square hole, b=a

XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> hv2(b = 0.1)             // just to confirm square hole formula
→ 0.00996658845699

XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)

Since taylor series all have even powers, let dimensionless x = b^2

XCas> pade(x - x^2/3 - 7*x^3/90, x, 3,3)
→ (17*x^2-30*x)/(7*x-30)

XCas> r(x) := x*(17*x-30)/(7*x-30)
XCas> r(0.1^2)                     // rough estimate for small 0.1×0.1 hole
→ 0.00996658870698

[Albert Chan]

(06-16-2020 12:35 PM)Albert Chan Wrote:  For unit diameter (d=1), with square hole, b=a

XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)

hv2 ≈ b^2 * √(1 - (2b^2)/3) = b^2 * √(1 - (1-k^2)/3)

Remove b=a requirement: hv2 ≈ a*b * √(1 - (a^2+b^2)/3)
Remove d=1 requirement: hv2 ≈ a*b * √(d^2 - (a^2+b^2)/3)

It is not as accurate hv2 pade approximation, but likely good enough
Previous example, d = 1, a = b = 0.1

lua> d, a, b = 1, 0.1, 0.1
lua> a*b * sqrt(d*d - (a*a+b*b)/3)
0.009966610925150703

Hole volume with 5 good digits, over-estimated.

Last term (estimated average height) is equivalent to RMS(d, d, sqrt(d*d-a*a-b*b))
Repalce RMS with arithemetic mean, we get the under-estimated hole volume.

lua> a*b * (2*d + sqrt(d*d-a*a-b*b)) / 3
0.00996649831220389

This is not as good as RMS version, does not reduce computation, thus not recommended.

[DM48]

To be clear, my son's integral solves for the volume left behind in the sphere after the hole has been removed. I have listed three volumes below. I wasn't sure if I was clear on this.

d = 12
a = 2
b = 2

Volume of sphere solid = 904.779

Volume of sphere with hole removed = 857.0638

Volume of hole removed = 47.715

Thoughts?

[Albert Chan]

(11-05-2023 04:28 PM)DM48 Wrote:  d = 12
a = 2
b = 2
...

Volume of sphere solid = 904.779

Volume of sphere with hole removed = 857.0638

Volume of hole removed = 47.715

sphere volume = pi*d^3/6, looks good --> we check removed hole only.

This is a good test of my recent estimate, lo < (average height) < hi

lua> d, a, b = 12, 2, 2
lua> lo = (2*d + sqrt(d*d-a*a-b*b))/3
lua> hi = sqrt(d*d - (a*a+b*b)/3)
lua> a*b*lo, a*b*hi
47.5492050529208      47.55347866700536

Your hole volume is outside this range. sphere (with hole) volume formula is incorrect.
Perhaps you can show steps to derive the formula?

Here is Derive6 calculated hole volume

#1: r := 6
#2: a := 2
#3: b := 2
#4: 8·∫(∫(√(r·r - x·x - y·y), x, 0, a/2), y, 0, b/2)

Shift-Enter (or click [✔≈])      → 47.55262983

Comment:

8 = 2^3, we could scale-up all 3 dimensions, to remove factor 8

#5 d := 2*r
#6: ∫(∫(√(d·d - x·x - y·y), x, 0, a), y, 0, b)      → 47.55262983

[DM48]

This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.

[rprosperi]

(11-05-2023 08:01 PM)DM48 Wrote:  This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.

I'd say you can bring down the size of that image to about 10-15% of it's current size, it's silly the way it is now.

[DM48]

(11-05-2023 08:27 PM)rprosperi Wrote:  

(11-05-2023 08:01 PM)DM48 Wrote:  This gave us a reason to finally get somewhat serious with 3D AutoCAD. 857.22 is correct. My son's formula is off a hair. He is going to work on the derivation of this formula for anyone to scrutinize possibly explain why it is off.

I'd say you can bring down the size of that image to about 10-15% of it's current size, it's silly the way it is now.

I have tried to resize it with img=##x## but it does not display.

[rprosperi]

(11-05-2023 09:52 PM)DM48 Wrote:  I have tried to resize it with img=##x## but it does not display.

Not sure what's going on, as Joe said maybe there are limits. I wasn't suggesting you make it unusable, it just looked a bit silly the way it was. Feel free to restore it if that's the only way it's useful to see the details.

[Albert Chan]

(11-05-2023 04:16 PM)Albert Chan Wrote:  

(06-16-2020 12:35 PM)Albert Chan Wrote:  For unit diameter (d=1), with square hole, b=a

XCas> k := sqrt(1-2*b*b)
XCas> hv2 := b*b*k/3 - atan(b*b/k)/3 + b*(1-b*b/3)*atan(b/k)
XCas> taylor(hv2,b,10)
→ b^2 - b^4/3 + -7*b^6/90 + -3*b^8/70 + -83*b^10/2520 + b^12*order_size(b)

hv2 ≈ b^2 * √(1 - (2b^2)/3) = b^2 * √(1 - (1-k^2)/3)

Remove b=a requirement: hv2 ≈ a*b * √(1 - (a^2+b^2)/3)
Remove d=1 requirement: hv2 ≈ a*b * √(d^2 - (a^2+b^2)/3)

To show hole volume is over-estimated, we check sign of errors

lua> taylor(hv2 - b^2 * √(1 - (2b^2)/3), b, 10)

- b^6/45 - 23*b^8/945 - 143*b^10/5670 + b^12*order_size(b)

Negative under-estimated error → RMS formula over-estimated hole volume.

Quote:Last term (estimated average height) is equivalent to RMS(d, d, sqrt(d*d-a*a-b*b))
Repalce RMS with arithemetic mean, we get the under-estimated hole volume.

height of rect cap (both of them), h = 1-k      // 1 is diameter, k is corner height

mean([1,1,k]) = 1 - h/3 = (1-h) + 2/3*h

MEAN formula under-estimated hole volume, if rect cap mean height > 2/3*h

XCas> ratio := (hv2/(b*b)-k) / (1-k)
XCas> taylor(ratio, b, 8, polynom)

2/3 + 4/45*b^2 + 5/63*b^4 + 23/252*b^6 + 1997/16632*b^8

Small hole, h and b^2 about the same size: k = 1-h = sqrt(1-2*b^2) ≈ 1-b^2
Ratio converge faster (smaller coefficients), if we use variable h.

XCas> taylor(ratio(b=sqrt((1-(1-h)^2)/2)), h, 6, polynom)

2/3 + 4/45*h + 11/315*h^2 + 1/84*h^3 + 25/8316*h^4



Compare mean height ratio with spherical cap.

V = pi*h/6 * (3*a^2 + h^2)
A = pi*a^2

ratio = V/A/h = 1/2 + (h/a)^2/6 > 1/2

[Albert Chan]

Code:

function hv2(d,a,b) -- sphere, rectanglar hole removed volume
    local d2, a2, b2 = d*d, a*a, b*b
    local k = sqrt(d2 - a2 - b2)
    d = a*b*k - d*d2*atan((a*b)/(d*k))
    return (2*d + a*(3*d2-a2)*atan(b/k) + b*(3*d2-b2)*atan(a/k)) / 6
end

-- fast estimate: m1 < hv2 < m2
function hv2_m2(d,a,b) return a*b * sqrt(d*d - (a*a+b*b)/3) end
function hv2_m1(d,a,b) return a*b * (2*d + sqrt(d*d-a*a-b*b))/3 end

Although estimation formula assumed square hole, inequality work with rectangular hole too.
For estimate, hv2_m2(d,a,b) is recommended: cheaper to compute, better accuracy.

[DM48]

Link to derivation of formula.


\( \int_{0}^{a} \int_{0}^{\sqrt{d^2-y^2-b^2}}\sqrt{d^2-y^2-x^2}-b\text{ }dx\text{ }dy\text{ }+\int_{a}^{d}\int_{0}^{\sqrt{d^2-y^2}}\sqrt{d^2-y^2-x^2}\text{ }dx\text{ }dy \)


In[842]:= Clear[a, b, q, r, x]

d = 12
a = 2
b = 2
q = NIntegrate[
Integrate[
Sqrt[d^2 - y^2 - x^2] - b, {x, 0, Sqrt[d^2 - y^2 - b^2]}], {y, 0,
a}]
r = NIntegrate[
Integrate[Sqrt[d^2 - y^2 - x^2], {x, 0, Sqrt[d^2 - y^2]}], {y, a, d}]
q + r


Out[848]= 857.2260543948448

a=4
b=2
d=12

Out[849]=811.0424359298954

[Albert Chan]

Hi, DM48

Nice to know there is another way to get volume:

sphere - hole = (x-strip - hole) + (2 spherical caps)

Since we have direct formula for spherical cap, 2nd integral can be skipped
2 caps = 2 * pi/3*h^2*(3r-h) = pi/12*(2h)^2*(3d-2h) = pi/12*(d-a)^2*(2d+a)

10 DEF FNH(X) @ N=N+1 @ FNH=SQRT(H-X*X) @ END DEF
20 DEF FNZ(Y) @ H=D*D-Y*Y @ FNZ=INTEGRAL(0,FNH(B),P,FNH(IVAR)-B) @ END DEF
30 INPUT "D,A,B = ";D,A,B @ N=0 @ P=1E-5
40 I1=INTEGRAL(0,A,P,FNZ(IVAR)) @ I2=PI/12*(D-A)^2*(2*D+A)
50 DISP I1;"+";I2;"=";I1+I2,N @ GOTO 30

>run
D,A,B = 12,2,2
 176.547653773 + 680.678408277 = 857.22606205       960
D,A,B = 12,2,10
 17.2338339713 + 680.678408277 = 697.912242248      992
D,A,B = 12,10,2
 662.307552043 + 35.6047167408 = 697.912268784      1888

Numerically, thinner, less curvy strip (A ≤ B) converge better.

For comparison, below patch directly get (-hole) + sphere

>20 DEF FNX(Y) @ H=D*D-Y*Y @ FNX=INTEGRAL(0,B,P,FNH(IVAR)) @ END DEF
>40 I1=-INTEGRAL(0,A,P,FNX(IVAR)) @ I2=PI/6*D*D*D

>run
D,A,B = 12,2,2
-47.552634005  + 904.778684234 = 857.226050229      225
D,A,B = 12,2,10
-206.866458199 + 904.778684234 = 697.912226035      465
D,A,B = 12,10,2
-206.866453454 + 904.778684234 = 697.91223078       465

[DM48]

The derivation of the formula with simplified equation.

Derivation PDF Direct Link

[Albert Chan]

Hi, DM48

Double integral reduced to single. Nice!
I simplified a bit, using Arc SOHCAHTOA method

10 DEF FNZ(Y) @ N=N+1 @ Y=D*D-Y*Y @ FNZ=Y*ACOS(B/SQR(Y)) @ END DEF
20 INPUT "D,A,B = ";D,A,B @ N=0 @ P=1E-5
30 I1=INTEGRAL(0,A,P,FNZ(IVAR))/2 @ C=SQR(D*D-A*A-B*B)
40 I2=PI/12*(D-A)^2*(2*D+A) - B/4*(A*C+(D*D-B*B)*ATAN(A/C))
50 DISP I1;"+";I2;"=";I1+I2,N @ GOTO 20

>run
D,A,B = 12,2,2
 200.09879005  +  657.127264322 = 857.226054372     15
D,A,B = 12,2,10
 82.547121343  +  615.365121066 = 697.912242409     31
D,A,B = 12,10,2
 764.410058495 + -66.4978178412 = 697.912240654    31

Below is direct formula for I1 (c is corner height, see code)

>a*(3*d^2-a^2)*pi/12 - a*b*c/12
 858.633041958
>res - b*(3*d^2+b^2)*atan(a/c)/12
 785.468340817
>res - a*(3*d^2-a^2)*atan(b/c)/6
 615.99486321
>res + d^3*atan(a*b/(c*d))/3
 764.410059992