suzanne.soy/FreeCAD
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Improved free surface evolving method in order to support long simulations.

Browse Source
This commit is contained in:
Jose Luis Cercós Pita 2012-08-10 16:57:58 +02:00 committed by Yorik van Havre
parent 060d63c6c9
commit b349147ea2
4 changed files with 38 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
src/Mod/Ship/simRun/Sim/fsEvolution.py
View File

@ -45,33 +45,19 @@ class simFSEvolution:
@param t Actual time (without adding dt).
"""
self.fs = fs
# Allocate memory
nx = self.fs['Nx']
ny = self.fs['Ny']
nF = nx*ny
# Evaluate potential gradients
grad = self.evaluateGradient()
# Integrate variables
# In order to improve results in really long simulations free surface
# will performed considering external waves and second order effects
# in two different ways. First external waves at time t will be
# substracted, then second order waves will be computed, and finally
# external waves at t+dt will be added.
for i in range(0,nx):
for j in range(0,ny):
# Get value at pos using characteristics method
gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j])
gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal
# Velocity potential
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \
dt*self.fs['accPot'][i,j] + \
0.5*dt*dt*grav*self.fs['pos'][i,j][2]
# Acceleration potential. This is really hard to simulate
# accurately due to numerical diffusion of the function, so
# external waves, and diffracted waves will be computed
# in two different ways:
# * External waves will be considered analitically,
# substracting waves at t, and adding waves at t+dt
# * Second order waves will be computed substracting external
# waves to free surface height, and then imposing boundary
# condition.
pos = np.copy(self.fs['pos'][i,j])
# Substract external waves at time t.
for w in waves['data']:
A = w[0]
T = w[1]
@ -81,16 +67,39 @@ class simFSEvolution:
k = 2.0*np.pi/wl
frec = 2.0*np.pi/T
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
# Substract external waves height in order to know second
# order waves free surface amplitude.
amp = A*np.sin(k*l - frec*(t+dt) + phase)
amp = A*np.sin(k*l - frec*t + phase)
pos[2] = pos[2] - amp
# Compute analitic external waves acceleration potential
amp0 = grav*A*np.cos(k*l - frec*t + phase)
amp1 = grav*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp0 + amp1
# Now impose free surface boundary condition
# self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + grav*pos[2]
amp = - grav/frec*A*np.sin(k*l - frec*t + phase)
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] - amp
amp = grav*A*np.cos(k*l - frec*t + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] - amp
# Now compute second order waves using position copy,
# where external waves are excluded, in order impose
# free surface boundary condition relative to second
# order phenomena.
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + \
dt*self.fs['accPot'][i,j]
# self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + \
# grav*pos[2]
# Restore external waves to velocity and acceleration
# potentials.
for w in waves['data']:
A = w[0]
T = w[1]
phase = w[2]
heading = np.pi*w[3]/180.0
wl = 0.5 * grav / np.pi * T*T
k = 2.0*np.pi/wl
frec = 2.0*np.pi/T
l = pos[0]*np.cos(heading) + pos[1]*np.sin(heading)
amp = - grav/frec*A*np.sin(k*l - frec*(t+dt) + phase)
self.fs['velPot'][i,j] = self.fs['velPot'][i,j] + amp
amp = grav*A*np.cos(k*l - frec*(t+dt) + phase)
self.fs['accPot'][i,j] = self.fs['accPot'][i,j] + amp
# Update free surface point position
gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j])
gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal)
self.fs['pos'][i,j][2] = self.fs['pos'][i,j][2] + dt*gradVal
# Impose values at beach (far free surface)
for i in range(0,nx):
for j in [0,ny-1]:
@ -107,15 +116,12 @@ class simFSEvolution:
ny = self.fs['Ny']
nF = nx*ny
grad = np.ndarray((nF,3), dtype=np.float32)
FF = open('gradient', 'w')
for i in range(0,nx):
for j in range(0,ny):
pos = self.fs['pos'][i,j]
grad[i*ny+j] = self.gradientphi(pos)
gradVal = np.dot(np.abs(grad[i*ny+j]),grad[i*ny+j])
gradVal = np.copysign(np.sqrt(np.abs(gradVal)), gradVal)
FF.write('%g\t%g\n' % (pos[1], gradVal))
FF.close()
return grad
def gradientphi(self, pos):

BIN
src/Mod/Ship/simRun/Sim/fsEvolution.pyc
View File

Binary file not shown.

8
src/Mod/Ship/simRun/Simulation.py
View File

@ -100,7 +100,6 @@ class FreeCADShipSimulation(threading.Thread):
ny = FS['Ny']
msg = Translator.translate("\t[Sim]: Iterating...\n")
FreeCAD.Console.PrintMessage(msg)
count = 0
while self.active and self.t < self.endTime:
msg = Translator.translate("\t\t[Sim]: Generating linear system matrix...\n")
FreeCAD.Console.PrintMessage(msg)
@ -113,13 +112,6 @@ class FreeCADShipSimulation(threading.Thread):
fsEvol.execute(FS, waves, dt, self.t)
self.t = self.t + dt
FreeCAD.Console.PrintMessage('t = %g s\n' % (self.t))
count = count+1
FF = open('%d' % (count), 'w')
i=1
for j in range(0,ny):
FF.write("%g\t%g\t%g\t%g\n" % (FS['pos'][i,j,1], FS['pos'][i,j,2],
FS['velPot'][i,j], FS['accPot'][i,j]))
FF.close()
# Set thread as stopped (and prepare it to restarting)
self.active = False
threading.Event().set()

BIN
src/Mod/Ship/simRun/Simulation.pyc
View File

Binary file not shown.