Find equations/theory for SLAB subroutine

[EmilSoleymani]

Find equations/theory for SLAB subroutine in lines 218-266 of the FORTRAN code. More info here in issue #14.

[EmilSoleymani]

This comment and this show my analysis of the subroutine, and the Boussineq equations I suspect are used in the code. I was hoping you could confirm which formula exactly is used.

[smiths]

@EmilSoleymani, could you be more specific about which lines of code you mean? Maybe this is a good question for @dstolle?

[EmilSoleymani]

I updated the link

[EmilSoleymani]

@dstolle @smiths

I have translated the original FORTRAN code of this subroutine into Julia in the wip-slab-calc branch. The code can be found here.

vdisp/src/OutputFormat/CalculationBehaviour.jl

Lines 358 to 395 in a8543ed

	function getSurchargePressure(behaviour, P::Array{Float64}, PP::Array{Float64})
	# foundationDepth = behaviour.bottomPointIndex * behaviour.dx
	dxx = 0.0
	wt = PP[behaviour.bottomPointIndex]
	pressure1 = behaviour.appliedPressure - wt
	pressure = pressure1
	for i=behaviour.bottomPointIndex:behaviour.nodalPoints
	if dxx < 0.01
	P[i] += pressure
	dxx += behaviour.dx
	continue
	end
	# material = behaviour.soilLayerNumber[i-1]
	if behaviour.foundation == "LongStripFooting"
	db = dxx / behaviour.foundationWidth
	ps = (db < 2.5 && behaviour.center) ? -0.28*db : -0.157 - 0.22*db
	ps *= 10
	P[i] += pressure*ps
	else
	basePressure = (behaviour.center) ? pressure : pressure/4
	length = (behaviour.center) ? behaviour.foundationLength/2 : behaviour.foundationLength
	width = (behaviour.center) ? behaviour.foundationWidth/2 : behaviour.foundationWidth
	ve2 = (length^2 + width^2 + dxx^2)/(dxx^2)
	ve = sqrt(ve2)
	an = length*width/(dxx^2)
	enm = (2 * an * ve / (ve2 + an^2)) * (ve2+1)/ve2
	fnm = (2 * an * ve / (ve2 - an^2))
	ab = (fnm < 0) ? atan(fnm) + pi : atan(fnm)
	P[i] += basePressure * (enm+ab)/pi
	end
	dxx += behaviour.dx
	end
	return P
	end

I am having a hard time matching the code to formulas from the textbook. I am getting really close, but I can't show a 1:1 correspondence. It seems the SLAB subroutine is using the Boussinesq equation for the stress at a point under the corner of a rectangular area, as seen on p.141 of Settlment's Analysis textbook:

Chunks of the code seem very alike to the above screenshot, however I can't seem to prove that they match exactly. For reference, here is their code:

vdisp/Resources/fortran_src/vdisp.for

Lines 242 to 247 in 3432b31

	50 VE2=(BL**2.+BW**2.+DXX**2.)/(DXX**2.)
	VE=VE2**0.5
	AN=BL*BW/(DXX**2.)
	AN2=AN**2.
	ENM=(2.*AN*VE/(VE2+AN2))*(VE2+1.)/VE2
	FNM=2.*AN*VE/(VE2-AN2)

The variable FNM is what gets passed in as the argument to the ATAN function. Looking at the equation, this matches ab/zC. BL is a, BW is b, and DXX is z. Here are screenshots of me trying to equate FNM to ab/zC:

Would it be better to copy how they did their FORTRAN code in our Julia (which I have already done), or to calculate exactly the formula given in the screenshots from the textbook.

Furthermore, everything discussed above was what happens for Rectangular Slab foundations. For Long Strip Footings, there is a different calculation used which I cannot find in any resources:

vdisp/src/OutputFormat/CalculationBehaviour.jl

Lines 371 to 376 in a8543ed

	if behaviour.foundation == "LongStripFooting"
	db = dxx / behaviour.foundationWidth
	ps = (db < 2.5 && behaviour.center) ? -0.28*db : -0.157 - 0.22*db
	ps *= 10
	P[i] += pressure*ps
	else

I couldn't think of any better/more meaningful names for all these intermediate variables (like ps, ve, ve2, an, ...), so I am hoping for some suggestions as well.

Also, the original fortran code calculates some variables that are never used:

vdisp/Resources/fortran_src/vdisp.for

Line 228 in 3432b31

ANBX=FLOAT(NBX)*DX

vdisp/Resources/fortran_src/vdisp.for

Line 234 in 3432b31

MTYP=IE(I-1,1)

I had originally copied these in the code, only to realize they are never used. I commented them out incase there is some use for them.

There is also this weird code:

vdisp/Resources/fortran_src/vdisp.for

Lines 230 to 231 in 3432b31

	BPRE1=Q-WT
	BPRES=BPRE1

I have copied it just like this (BPRE1 is called pressure1, BPRES is called pressure). However, I don't see the need for BPRE1 since they never update BPRES and always just use its value.

Finally, I am pretty confident that the following code uses the 2:1 method for calculating stress at center of slab (divides length and with by 2 and pressure by 4):

vdisp/src/OutputFormat/CalculationBehaviour.jl

Lines 377 to 379 in a8543ed

	basePressure = (behaviour.center) ? pressure : pressure/4
	length = (behaviour.center) ? behaviour.foundationLength/2 : behaviour.foundationLength
	width = (behaviour.center) ? behaviour.foundationWidth/2 : behaviour.foundationWidth

[dstolle]

Gentlemen – I will take a look at this. Dieter From: Emil Soleymani ***@***.***> Sent: Monday, June 6, 2022 3:59 PM To: smiths/vdisp ***@***.***> Cc: Stolle, Dieter ***@***.***>; Mention ***@***.***> Subject: Re: [smiths/vdisp] Find equations/theory for SLAB subroutine (Issue #52) @dstolle<https://github.com/dstolle> @smiths<https://github.com/smiths> I have translated the original FORTRAN code of this subroutine into Julia in the wip-slab-calc branch. The code can be found here<https://github.com/smiths/vdisp/blob/wip-slab-calc/src/OutputFormat/CalculationBehaviour.jl#L358-L395>. https://github.com/smiths/vdisp/blob/a8543ed448aa3053f6261c0600265bd2dd47bf8b/src/OutputFormat/CalculationBehaviour.jl#L358-L395 I am having a hard time matching the code to formulas from the textbook. I am getting really close, but I can't show a 1:1 correspondence. It seems the SLAB subroutine is using the Boussinesq equation for the stress at a point under the corner of a rectangular area, as seen on p.141 of Settlment's Analysis textbook<https://www.publications.usace.army.mil/Portals/76/Publications/EngineerManuals/EM_1110-1-1904.pdf>: [image]<https://user-images.githubusercontent.com/97846411/172234386-9670a070-cd80-4509-8cf9-447c62dee000.png> Chunks of the code seem very alike to the above screenshot, however I can't seem to prove that they match exactly. For reference, here is their code: https://github.com/smiths/vdisp/blob/3432b311995e18e32d7de5c7d1a92ad66710cacc/Resources/fortran_src/vdisp.for#L242-L247 The variable FNM is what gets passed in as the argument to the ATAN function. Looking at the equation, this matches ab/zC. BL is a, BW is b, and DXX is z. Here are screenshots of me trying to equate FNM to ab/zC: [image]<https://user-images.githubusercontent.com/97846411/172235746-a329bc28-af6f-4841-af22-09d743730162.png> [image]<https://user-images.githubusercontent.com/97846411/172235782-dce2dce4-c599-4631-bd47-1af17697965b.png> Would it be better to copy how they did their FORTRAN code in our Julia (which I have already done), or to calculate exactly the formula given in the screenshots from the textbook. Furthermore, everything discussed above was what happens for Rectangular Slab foundations. For Long Strip Footings, there is a different calculation used which I cannot find in any resources: https://github.com/smiths/vdisp/blob/a8543ed448aa3053f6261c0600265bd2dd47bf8b/src/OutputFormat/CalculationBehaviour.jl#L371-L376 I couldn't think of any better/more meaningful names for all these intermediate variables (like ps, ve, ve2, an, ...), so I am hoping for some suggestions as well. Also, the original fortran code calculates some variables that are never used: https://github.com/smiths/vdisp/blob/3432b311995e18e32d7de5c7d1a92ad66710cacc/Resources/fortran_src/vdisp.for#L228 https://github.com/smiths/vdisp/blob/3432b311995e18e32d7de5c7d1a92ad66710cacc/Resources/fortran_src/vdisp.for#L234 I had originally copied these in the code, only to realize they are never used. I commented them out incase there is some use for them. There is also this weird code: https://github.com/smiths/vdisp/blob/3432b311995e18e32d7de5c7d1a92ad66710cacc/Resources/fortran_src/vdisp.for#L230-L231 I have copied it just like this (BPRE1 is called pressure1, BPRES is called pressure). However, I don't see the need for BPRE1 since they never update BPRES and always just use its value. Finally, I am pretty confident that the following code uses the 2:1 method for calculating stress at center of slab (divides length and with by 2 and pressure by 4): https://github.com/smiths/vdisp/blob/a8543ed448aa3053f6261c0600265bd2dd47bf8b/src/OutputFormat/CalculationBehaviour.jl#L377-L379 [image]<https://user-images.githubusercontent.com/97846411/172237980-43897b72-f3b2-48b2-8330-eb3509543aa8.png> — Reply to this email directly, view it on GitHub<#52 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAGIDIWMC4GAW2GSLI3S463VNZKADANCNFSM5WGE4BGA>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

[smiths]

@EmilSoleymani good summary of your question. In the short term you should definitely use the translation of the FORTRAN code. We are confident that the FORTRAN code is correct. As we get a better understanding of the mapping between the theory and code, we may decide to change to something more understandable. As long as you continue to add test cases as you write the code, continuous integration should give us confidence that any changed code is still correct.

@dstolle that would be great if you could help @EmilSoleymani with this issue. When you reply, please use the web-interface for GitHub, rather than replying via e-mail. You can find the web-interface for this issue here:

#52

You should be able to just click on the hyperlink.

[smiths]

@dstolle, this isn't the proper way to use the issue tracker, but did you get the e-mail I sent about the meeting on Wednesday (June 8)? Emil and I are meeting at 1:00 pm in my office if you are available to join us.