Simplified explosive lens approximation method

So i was considering creating my own nuke simulator because of some perceived inaccuracies on the current nuclear payload designer. Still failing to find info on things like interstage design and electric neutron initiators, but i can get a pretty decent idea on imploder design already - bunch of explosive lenses(and very precise detonators) form the implosion wavefront, whose only purpose is detonating a high explosive hollow spherical charge(that performs the actual implosion) evenly enough to enable the weapon to achieve critical mass. The ingame model ignores this high explosive layer and instead assumes the slow explosive is directly responsible for compressing the payload.

There's not much on the open literature regarding explosive lens design, and the closest thing i could find was this paper detailing simulating one. It isn't directly appliable to nuclear weapon design as it studies plane wave generators instead of your usual spherical/prolate/ovoid implosion wave, but it can be adapted, and also reveals you have a significant pressure gradient across the wavefront, which if thrown directly into your compression design would likely cause a design using the explosive lens alone to fizzle due to irregular compression.

Given how the pressure wavefront isn't shown to significantly deform on later stages, it gave me the idea of ignoring it altogether and representing the wavefront as a series of intersecting slow explosive explosions generated by a fast explosive liner moving faster across the lens. You can even reverse-simulate your intended wavefront over discrete time periods and obtain your fast explosive liner's shape by doing it! (spoiler: it approximates as a revoluted truncated exponential spiral). I've wrote a program to do this(in Lua, of all things - shame on me):

2-lens design:

7-lens design:

1-lens design pictured before using up all the lens(for all those nations that haven't developed exploding bridgewire or slapper detonators yet):



Here's the source code for anyone interested (using Paul Kuchenko's turtle library - the easiest way to get it running is probably by running ZeroBrane Studio):



require "turtle"

ImploderChargeSize=80 --a 1 pixel/meter scale is terrible for display purposes with normal nuke sizes.
ImploderInnerSize=50 --let's just assume time is going off several orders of magnitude faster instead

Timestep=0.00025 --assuming a 30cm radius imploder charge, the time scale is actually 100 times faster than this.
--MaxSimulationSteps=1000

--Comp B
FastExplosiveVelocity = 8050
--Baratol. Seems like a pretty shitty explosive overall - going slower massively improves lens performance.
SlowExplosiveVelocity = 4900


NumberOfLenses =2--note this is still in 2D and would likely have issues accurately representing designs with over 2 lenses as this doesn't try to take a cross-section of them


size(600,600)

 if FastExplosiveVelocity<=SlowExplosiveVelocity then
--   error("The fast explosive doesn't detonate faster than the slow explosive. This makes shaping an explosive lens impossible.")
 end
SlowExplosiveVelocity=SlowExplosiveVelocity*Timestep
FastExplosiveVelocity=FastExplosiveVelocity*Timestep

function rotate(x,y,a)
 --rotates 2d coordinates x,y by a
 local x_prime=0
 local y_prime=0
 x_prime = (math.cos(a)*x)+(math.sin(a)*y)
 y_prime = (-math.sin(a)*x)+(math.cos(a)*y)
 return x_prime,y_prime
end


function circle_circle_intersection(x0, y0, r0, x1, y1, r1)
   --returns the intersection points between 2 circles.
   --code ported from http://paulbourke.net/geometry/circlesphere/tvoght.c
 local a, dx, dy, d, h, rx, ry;
 local x2, y2;
 --dx and dy are the vertical and horizontal distances between the circle centers.
 dx = x1 - x0;
 dy = y1 - y0;
 --Determine the straight-line distance between the centers.
 d = math.sqrt((dy*dy) + (dx*dx))
 --d = hypot(dx,dy); // Suggested by Keith Briggs

 --Check for solvability
 if (d > (r0 + r1)) then
   --no solution. circles do not intersect
   error("no solution. circles do not intersect")
 end
 if (d < math.abs(r0 - r1)) then
   --no solution. one circle is contained in the other
   error("no solution. one circle is contained in the other")
 end

   --[[
     'point 2' is the point where the line through the circle
     intersection points crosses the line between the circle centers.
     ]]

 --Determine the distance from point 0 to point 2.
 a = ((r0*r0) - (r1*r1) + (d*d)) / (2*d)
 --Determine the coordinates of point 2.
 x2 = x0 + (dx * a/d);
 y2 = y0 + (dy * a/d);
 --Determine the distance from point 2 to either of the intersection points.
 h = math.sqrt((r0*r0) - (a*a));
 --Now determine the offsets of the intersection points from point 2.
 rx = -dy * (h/d)
 ry = dx * (h/d)
 --Determine the absolute intersection points
 xi = x2 + rx
 xi_prime = x2 - rx
 yi = y2 + ry
 yi_prime = y2 - ry
 return xi,yi,xi_prime,yi_prime
 end

--draws the compression explosive layer
oval(0,0,ImploderChargeSize,ImploderChargeSize,"#FFFF00")
oval(0,0,ImploderInnerSize,ImploderInnerSize,"#FFFFFF")


--sets the final fast explosive shockwave position
ContourPosX = math.cos(math.pi-(math.pi/NumberOfLenses))*ImploderChargeSize
ContourPosY = math.sin(math.pi-(math.pi/NumberOfLenses))*ImploderChargeSize
--print(math.deg(math.pi/NumberOfLenses))

ContourPosList={}
table.insert(ContourPosList,{ContourPosX,ContourPosY})

local n=ImploderChargeSize
pncl("#F0F0F0")
repeat
 n=n+SlowExplosiveVelocity
 --gives some time to enjoy the view
 --wait(Timestep*20)
 --draws your expanding wavefront
 oval(0,0,n,n)
 --finds out where the wavefront has moved since then, then records it
   local intersectionpoints= table.pack(circle_circle_intersection(0,0,n, ContourPosX,ContourPosY,FastExplosiveVelocity))
   --pncl("#000000")
   --line(ContourPosX,ContourPosY,intersectionpoints[1],intersectionpoints[2])
   ContourPosX = intersectionpoints[1]
   ContourPosY = intersectionpoints[2]
   table.insert(ContourPosList,{ContourPosX,ContourPosY})
   
until ContourPosY<=0

-- now let's wrap up loose ends
do --i don't know what variable names i've used anymore, better wrap this in a separate scope just in case
 local n = table.maxn(ContourPosList)
 local dx = ContourPosList[n][1]-ContourPosList[n-1][1]
 local dy = ContourPosList[n][2]-ContourPosList[n-1][2]
 --print(ContourPosList[n-1][1],ContourPosList[n-1][2])
 --print(ContourPosList[n][1],ContourPosList[n][2])
 --print(dx,dy)
 scalingratio = -ContourPosList[n-1][2]/dy
 --print(scalingratio)
 dx =dx*scalingratio
 dy =dy*scalingratio
 ContourPosList[n]={ContourPosList[n-1][1]+dx, ContourPosList[n-1][2]+dy}
end

--draws an outline assuming a spherical lens assembly. you could cut quite a bit of the fast eplosive in practice.
--figuring out how much you can cut and still retain a concave shape is an exercise left to the reader.
oval(0,0,ContourPosList[#ContourPosList][1],ContourPosList[#ContourPosList][1])

--and finally, draw the explosive lens contour.
pncl("#000000")
for n=2,table.maxn(ContourPosList) do
 local pos1x = ContourPosList[n-1][1]
 local pos1y = ContourPosList[n-1][2]
 local pos2x = ContourPosList[n][1]
 local pos2y = ContourPosList[n][2]
   for a=(2*math.pi)/NumberOfLenses, 2*math.pi, (2*math.pi)/NumberOfLenses do
     local rotatedpos1x,rotatedpos1y = rotate(pos1x,pos1y,a)
     local rotatedpos2x,rotatedpos2y = rotate(pos2x,pos2y,a)
     line(rotatedpos1x,rotatedpos1y,rotatedpos2x,rotatedpos2y)
     --mirrors the generated curve along the rotation axis
     rotatedpos1x,rotatedpos1y = rotate(pos1x,-pos1y,a)
     rotatedpos2x,rotatedpos2y = rotate(pos2x,-pos2y,a)
     line(rotatedpos1x,rotatedpos1y,rotatedpos2x,rotatedpos2y)

 end
end

--waits a bit
wait(3)
--now, time for some (likely extremely laggy due to limited drawing capabilities) simulation!
--store a snapshot so we don't have all wavefronts overlaying eachother
local snapshot = snap()
--turns off automatic display every function
updt(false)


--creates the shockwaves
pncl("#888888")
local time=1
repeat
 undo(snapshot)
 for n=table.maxn(ContourPosList),table.maxn(ContourPosList)-time,-3 do -- yeah, i know, magic number. whatever.
   local i = (table.maxn(ContourPosList)-time)-n
   local pos1x = ContourPosList[n][1]
   local pos1y = ContourPosList[n][2]
     for a=(2*math.pi)/NumberOfLenses, 2*math.pi, (2*math.pi)/NumberOfLenses do --cut this to a single lens if you have perf issues
       local rotatedpos1x,rotatedpos1y = rotate(pos1x,pos1y,a)
       oval(rotatedpos1x,rotatedpos1y,SlowExplosiveVelocity*i,SlowExplosiveVelocity*i)
       --mirrors the generated curve along the rotation axis, again
         rotatedpos1x,rotatedpos1y = rotate(pos1x,-pos1y,a)
       oval(rotatedpos1x,rotatedpos1y,SlowExplosiveVelocity*i,SlowExplosiveVelocity*i)
   end
 end
 updt() --force an update
 time=time+1
until time==table.maxn(ContourPosList)
wait()


It's interesting to see how smaller you can make your explosive lens by having a better fast explosive:slow explosive detonation velocity ratio or by using more lenses. The two-lens design used ingame doesn't seem that practical anymore due to lens mass - especially when taking a look at a W80 layout greenpeace unwittingly released some time ago.